C+
C NAME:
C	lm_fit
C PURPOSE:
C	To make a Levenberg-Marquardt fit of a simple function.
C
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	run lm_fit
C INPUTS:
C	Main program
C OUTPUTS:
C	Fits parameters and errors of the function fit
C FUNCTIONS/SUBROUTINES:
C	mrqmin.f
C INCLUDE:
C COMMON BLOCKS:
C SIDE EFFECTS:
C RESTRICTIONS:
C PROCEDURE:
C	Latest 10/11/13 work
C	logical $CAM.
C MODIFICATION HISTORY:
C	October-2013, Bernard Jackson (UCSD)
C-

	program lm_fit
C
	parameter      (NdataM	=  50,		! Max # data points
     &			MAM	=  2,		! Max # of coefficients
     &			MfitM	=  2)		! Max # of coefficients fit
C
	real		X	(NdataM),	! Values of numbers input (spectral values)
     &			Xin	(NdataM),	! Input set of numbers (same as above but passed to mrqmin)
     &			Ydat	(NdataM),	! Data values at each value of Xin
     &			Y	(NdataM),	! The function values
     &			Sig	(NdataM),	! Individual standard deviations of the values fit
     &			A	(MAM),		! Coefficients to be fit
     &			Covar	(MfitM,MfitM),	! scratch array	
     &			Alpha	(MfitM,MfitM)	! scratch array
C
	integer		ListA	(MAM)
	external	Funcs			!  Put in to declare an external function by B. Jackson - needed here

	Ndata 	= 10				! # of actual data points
	MA 	= MAM				! # of coefficients
	Mfit	= MfitM				! # number of coeficients fit
	NCA	= Mfit				! a # greater or equal to Mfit 
C	Mfit = 1
C
C Set up the data points to be fit by Levenberg-Marquardt.
C
	a1 = 1.0				! This is the actual value of A(1)
	b1 = 2.0				! This is the actual value of A(2)

	do i=1,Ndata
	  X(i)    = float(i)
	  Xin(i) = X(i)
	  Ydat(i) = a1*X(i)**2 + b1
C	  Sig(i) = 1.0
	  Sig(i) = 1.0*X(i)**2
C	  Sig(i) = 0.5*X(i)**2
	end do
C	Ydat(2) = 4.0
	Ydat(4) = 16.0
	Ydat(8) = 68.0
C	Ydat(9) = 79.0
	Lista(1) = 1
	Lista(2) = 2
C
C End of set up of points to be fit.
C
C Guess the values of the coefficients.
C
	A(1)	= 1.3
	A(2)	= 1.8
	Alamda	= -1.0
C
C Call the Minimization
C
	write(*,'(A,12F6.2)') 'Xin     =',(Xin(I),I=1,12) 
	Chisqold = 15.0
	Chisq	 = 10.0
	do i=1,15
	  write(*,'(A,I5,2F15.4,F15.10)') 'Enter Chisqold, Chisq  =',i, Chisqold, Chisq, Chisq-Chisqold	
	  if((Chisq-Chisqold).lt.-0.001.or.(Chisq-Chisqold).gt.0.0) then
 	    print *, 'Outside mngmin'
	    Chisqold = Chisq	    
	    call mrqmin(Xin,Ydat,Sig,Ndata,A,MA,ListA,Mfit,Covar,Alpha,NCA,Chisq,Funcs,Alamda)
	    print *, '  '
	    print *, 'Chisqold, Chisq  =',Chisqold, Chisq
	    print *, 'Alamda =',Alamda
C	    write(*,'(A,12F6.2)') 'Ydat  =',(Ydat(J),J=1,12) 
C	    write(*,'(A,12F6.2)') 'Xin   =',(Xin (J),J=1,12) 
	  end if
	end do
	Alamda = 0.0
	print *, '  '
	print *, 'Outside last mrqmin'
	call mrqmin(Xin,Ydat,Sig,Ndata,A,MA,ListA,Mfit,Covar,Alpha,NCA,Chisq,Funcs,Alamda) 
	print *, 'Program sucessfully completed'
	write(*,'(A,2F10.5)') ' Final values of A     =', A(1),A(2)
	write(*,'(A,4F10.5)') ' Final values of Covar =', Covar(1,1),Covar(1,2),Covar(2,1),Covar(2,2)
	stop
	end
C************************************************************************************************
	subroutine Funcs(X,A,Xmod,Dxda,NA)
C
	real	A	(NA),
     &		Dxda	(NA)
	real    er	(NA),
     &		AA	(20)
C
	call formula(X,A,Xmod,NA)
C	Xmod = A(1)*X**2 + A(2)
C
C Determine A derivatives
C
	er(1) = 0.001
	er(2) = 0.001
C
	call derivatives(X,A,Dxda,er,NA)
C
	write(*,'(A,5F12.4)') 'Inside Funcs 1', Xmod, A(1), A(2), Dxda(1), Dxda(2)
C
C	Dxda(1) = X**2
C	Dxda(2) = 1.0
C end A derivative determination
C
C	write(*,'(A,5F12.4)') 'Inside Funcs 2', Xmod, A(1), A(2), Dxda(1), Dxda(2)
C
	return
	end
C
C**************************************************************************************************
C
	subroutine formula(X,A,Xmod,NA)
C
	real	A	(NA)
C
C The formula that fits the data
C
	  Xmod = A(1)*X**2 + A(2)
C
	return
	end
C
C**************************************************************************************************
C
	subroutine derivatives(X,A,Dxda,er,NA)
C
	real	A	(NA),
     &		Dxda	(NA)
	real    er	(NA),
     &		AA	(20)
C
C The derivatives
C
C	Dxda(1) = X**2
C	Dxda(2) = 1.0
C

	do i = 1, NA
	  do j = 1, 20
	    AA(j) = A(j)
          end do
	  AA(i) = A(i) + er(i)
	  call formula(X,AA,XmodAP,NA)
	  AA(i) = A(i) - er(i)
	  call formula(X,AA,XmodAM,NA)
	  Dxda(i) = (XmodAP-XmodAM)/(2.0*er(i))
	end do
C
	return
	end