SUBROUTINE lfit(x,y,sig,ndat,a,ma,npc,funcs,ierr)
      implicit none

! ierr = 1, ok; =0 singular matrix

      INTEGER ma,ia(ma),npc,ndat,MMAX,ierr
      double precision chisq,a(ma),covar(npc,npc),
     &    sig(ndat),x(ndat),y(ndat)
      EXTERNAL funcs
      PARAMETER (MMAX=10000)
CU    USES covsrt,gaussj
      INTEGER i,j,k,l,m,mfit
      double precision sig2i,sum,wt,ym,afunc(MMAX),beta(MMAX)
      mfit=0
      ierr=1      
      do 11 j=1,ma
        ia(j)=1
        if(ia(j).ne.0) mfit=mfit+1
11    continue
      do 13 j=1,mfit
        do 12 k=1,mfit
          covar(j,k)=0.0d0
12      continue
        beta(j)=0.0d0
13    continue
      do 17 i=1,ndat
      
        call funcs(x(i),afunc,ma) 
        
        ym=y(i)
        if(mfit.lt.ma) then
          do 14 j=1,ma
            if(ia(j).eq.0) ym=ym-a(j)*afunc(j)
14        continue
        endif
        sig2i=1.0d0/sig(i)**2
        j=0
        do 16 l=1,ma
          if (ia(l).ne.0) then
            j=j+1
            wt=afunc(l)*sig2i
            k=0
            do 15 m=1,l
              if (ia(m).ne.0) then
                k=k+1
                covar(j,k)=covar(j,k)+wt*afunc(m)
              endif
15          continue
            beta(j)=beta(j)+ym*wt
          endif
16      continue
17    continue
      do 19 j=2,mfit
        do 18 k=1,j-1
          covar(k,j)=covar(j,k)
18      continue
19    continue
      call gaussj(covar,mfit,npc,beta,1,1,ierr)
      if(ierr.eq.0)then
!        singular matrix
         return
      endif
      
      j=0
      do 21 l=1,ma
        if(ia(l).ne.0) then
          j=j+1
          a(l)=beta(j)
        endif
21    continue
      chisq=0.0d0
      do 23 i=1,ndat
        call funcs(x(i),afunc,ma)
        sum=0.0d0
        do 22 j=1,ma
          sum=sum+a(j)*afunc(j)
22      continue
        chisq=chisq+((y(i)-sum)/sig(i))**2
23    continue
      call covsrt(covar,npc,ma,ia,mfit)
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.

      SUBROUTINE covsrt(covar,npc,ma,ia,mfit)
      implicit none
      INTEGER ma,mfit,npc,ia(ma)
      double precision covar(npc,npc)
      INTEGER i,j,k
      double precision swap
      do 12 i=mfit+1,ma
        do 11 j=1,i
          covar(i,j)=0.0d0
          covar(j,i)=0.0d0
11      continue
12    continue
      k=mfit
      do 15 j=ma,1,-1
        if(ia(j).ne.0)then
          do 13 i=1,ma
            swap=covar(i,k)
            covar(i,k)=covar(i,j)
            covar(i,j)=swap
13        continue
          do 14 i=1,ma
            swap=covar(k,i)
            covar(k,i)=covar(j,i)
            covar(j,i)=swap
14        continue
          k=k-1
        endif
15    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.

      SUBROUTINE gaussj(a,n,np,b,m,mp,ierr)
      implicit none
      INTEGER m,mp,n,np,NMAX,ierr
      double precision a(np,np),b(np,mp)
      PARAMETER (NMAX=10000)
      INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),
     &      ipiv(NMAX)
      double precision big,dum,pivinv
      ierr=1
      do 11 j=1,n
        ipiv(j)=0
11    continue
      do 22 i=1,n
        big=0.0d0
        do 13 j=1,n
          if(ipiv(j).ne.1)then
            do 12 k=1,n
              if (ipiv(k).eq.0) then
                if (dabs(a(j,k)).ge.big)then
                  big=dabs(a(j,k))
                  irow=j
                  icol=k
                endif
              else if (ipiv(k).gt.1) then
!                pause 'singular matrix in gaussj'
                ierr=0
                return
              endif
12          continue
          endif
13      continue
        ipiv(icol)=ipiv(icol)+1
        if (irow.ne.icol) then
          do 14 l=1,n
            dum=a(irow,l)
            a(irow,l)=a(icol,l)
            a(icol,l)=dum
14        continue
          do 15 l=1,m
            dum=b(irow,l)
            b(irow,l)=b(icol,l)
            b(icol,l)=dum
15        continue
        endif
        indxr(i)=irow
        indxc(i)=icol
        if (a(icol,icol).eq.0.0d0)then
!         pause 'singular matrix in gaussj'
          ierr=0
          return
        endif
        pivinv=1.0d0/a(icol,icol)
        a(icol,icol)=1.0d0
        do 16 l=1,n
          a(icol,l)=a(icol,l)*pivinv
16      continue
        do 17 l=1,m
          b(icol,l)=b(icol,l)*pivinv
17      continue
        do 21 ll=1,n
          if(ll.ne.icol)then
            dum=a(ll,icol)
            a(ll,icol)=0.0d0
            do 18 l=1,n
              a(ll,l)=a(ll,l)-a(icol,l)*dum
18          continue
            do 19 l=1,m
              b(ll,l)=b(ll,l)-b(icol,l)*dum
19          continue
          endif
21      continue
22    continue
      do 24 l=n,1,-1
        if(indxr(l).ne.indxc(l))then
          do 23 k=1,n
            dum=a(k,indxr(l))
            a(k,indxr(l))=a(k,indxc(l))
            a(k,indxc(l))=dum
23        continue
        endif
24    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.