      subroutine fit_g(ndat,x,y,sigy,nfit,a,covar,amat,chi2,q,subg,lpri)
c
c bernd berg, october 8, 2002.
c general fit with the levenberg-marquardt method.
c
c ndat:     number of data points.
c nfit:     number of fit paprameters.
c mfit:     maximum number of fit parameters.
c nax_iter: maximum number of steps (iterations).
c
c x:      arguments of data point, y(i)=y(x(i)), i=1,...,ndat.
c y:      y(i),    i=1,...,ndat, data points.
c sigy:   sigy(i), i=1,...,ndat, standard deviatons of the data points.
c a:      fit parameters;  last nfit-nfit values are held constant.
c subg:   user supplied fit subroutine  subg(xx,a,yfit,dyda,nfit).
c lpri:   if .true., information about iteration process is printed.
c eps:    (internal) stop parameter for the iteration.
c
c a:      array (vector) fit parameters.
c cova:   covariance matrix on output.
c amat:   curvature  matrix on output.
c chi2:   chi squared.
c q:      probability that observed deviations of data from best fit
c         are due to chance (assuming gaussian distribution of data)
c
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      parameter(iuo=6,mfit=20,max_iter=250,eps=ten**(-4))
      dimension x(ndat),y(ndat),sigy(ndat),a(nfit),covar(nfit,nfit),
     & amat(nfit,nfit), atry(mfit),b(mfit),dela(mfit),dyda(mfit)
      external subg
      if(nfit.gt.mfit) stop "fit_g: enlarge mfit."
      n0chi2=0
c
      xla=one/ten**3 ! initial value of lambda.
      call fit_gc(ndat,x,y,sigy,nfit,a,amat,b,dyda,chi2,subg)
      if(lpri) write(iuo,'(/," fit_g: chi2_in =",1g16.7,/)') chi2
      if(lpri) write(iuo,*) "  iter     chi2         dchi2/chi2"
c
      iter=0
1     iter=iter+1
      if(ndat.gt.nfit .and. chi2.lt.eps) then
        write(iuo,'(" check your error bars: chi2 =",g16.7)') chi2
        stop "fit_g: chi2 too small."
      end if
c
      do j=1,nfit
        dela(j)=b(j)
        do k=1,nfit
          covar(j,k)=amat(j,k)
        end do   
        covar(j,j)=amat(j,j)*(one+xla) ! marquard's definition.
      end do
c
      call mat_gau(covar,nfit,dela,1)
c
      do j=1,nfit
        atry(j)=a(j)+dela(j)
      end do
      call fit_gc(ndat,x,y,sigy,nfit,atry,covar,dela,dyda,chi2_try,subg)
      if(ndat.gt.nfit) dchi2=(chi2-chi2_try)/chi2
      if(ndat.eq.nfit) dchi2=chi2/eps
c 
      if(chi2_try.lt.chi2) then ! accept try
        xla=xla/ten ! and decrease lambda by the factor ten.
        chi2=chi2_try
        do j=1,nfit
          b(j)=dela(j) 
          a(j)=atry(j) 
          do k=1,nfit
            amat(j,k)=covar(j,k)
          end do
        end do
      else ! reject try 
        xla=ten*xla ! and enlarge lambda by the factor ten.
      endif
c
      if(lpri) write(iuo,'(i7,2g15.6)') iter,chi2,dchi2
      if(dchi2.lt.zero .or. dchi2.gt.eps) n0chi2=0
      if(dchi2.lt.zero .or. dchi2.gt.eps) goto 1
      if(dchi2.eq.zero) then
        n0chi2=n0chi2+1
        if(n0chi2.le.10) goto 1
        write(iuo,'(" fit_g warning: n0chi2.eq.10 reached.",/)')
      end if
c
      do j=1,nfit
        dela(j)=b(j)
        do k=1,nfit
          covar(j,k)=amat(j,k)
        end do   
      end do
      call mat_gau(covar,nfit,dela,1)
      q=-one
      ndgf=ndat-nfit
      if(ndgf.gt.0) q=one-gamma_p(half*ndgf,half*chi2)
c
      return
      end


      subroutine fit_gc(ndat,x,y,sigy,nfit,a,amat,b,dyda,chi2,subg)
c caculates the matrix amat and the vector b.
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      dimension x(ndat),y(ndat),sigy(ndat),amat(nfit,nfit),
     &          b(nfit),dyda(nfit)
      external subg
      do j=1,nfit
        b(j)=zero
        do k=1,j
          amat(j,k)=zero
        end do
      end do
      chi2=zero
      do i=1,ndat
        call subg(x(i),a,yfit,dyda,nfit)
        sigy2i=one/(sigy(i)*sigy(i))
        dely=y(i)-yfit
        do j=1,nfit
          das2i=dyda(j)*sigy2i
          b(j)=b(j)+das2i*dely
          do k=1,j
            amat(j,k)=amat(j,k)+das2i*dyda(k)
          end do
        end do
        chi2=chi2+dely*(dely*sigy2i)
      end do
      do j=2,nfit
      do k=1,j-1
        amat(k,j)=amat(j,k)
      end do
      end do
      return
      end
