PROGRAM fit
	USE analyze_timeseries
	USE analysis_functions
	USE make_model
	USE fit_library
	USE marquardt_library

	IMPLICIT NONE

	INTERFACE
		SUBROUTINE modelWrapper(X,A,Y,DYDA,NA)
			USE fit_library
			USE make_model

			REAL, DIMENSION(NA) :: A
			REAL, DIMENSION(13) :: X, Y
			REAL, DIMENSION(2) :: fRange, qyRange, zRange
			REAL, DIMENSION(13, NA) :: DYDA
			REAL :: fStep, qyStep, zStep
			INTEGER :: numRows, NA
		END SUBROUTINE modelWrapper
	END INTERFACE

	REAL, DIMENSION(:,:), POINTER :: model

	REAL, DIMENSION(:,:), POINTER :: dataFFT
	INTEGER :: numRowsFFT, numRowsModel, stringIdx, i
	REAL :: percentDiffTol, percentOutTol, sampleRate
	REAL :: alpha, AR, elongation, nu, theta, v, average, stdDev
	REAL, DIMENSION(2) :: chiSquared
	CHARACTER(LEN=64) :: inputFile, fftOutputFile, modelOutputFile, logString

	INTEGER :: MA, MFIT, NCA, NDATA
	REAL :: CHISQ, ALAMDA
	REAL, DIMENSION(:), ALLOCATABLE :: A, SIG, X, Y
	REAL, DIMENSION(:,:), ALLOCATABLE :: AMATRIX, COVAR
	INTEGER, DIMENSION(:), ALLOCATABLE :: LISTA

	CALL getTimeseriesArguments(inputFile, sampleRate, percentDiffTol, percentOutTol)
	CALL getModelArguments(alpha, AR, elongation, nu, theta, v)

	dataFFT => analyzeTimeseries(inputFile, sampleRate, percentDiffTol, percentOutTol, numRowsFFT)

	stringIdx = INDEX(inputFile,'T')
	
	WRITE(fftOutputFile, '(A4,A9,A5)') 'fft_', inputFile(stringIdx:stringIdx + 8), '.data' !27:35 for 3C48, 41:49 for 3C273
	CALL writeData(fftOutputFile, dataFFT, numRowsFFT)

	MA = 6
	MFIT = 4
	NCA = MA
	NDATA = 13
	CHISQ = 0.0
	ALAMDA = -0.1
	
	ALLOCATE(A(MA))
	ALLOCATE(AMATRIX(NCA,NCA))
	ALLOCATE(COVAR(NCA,NCA))
	ALLOCATE(LISTA(MA))
	ALLOCATE(SIG(NDATA))
	ALLOCATE(X(NDATA))
	ALLOCATE(Y(NDATA))

	IF (theta .EQ. 0.0) THEN
		theta = 0.1
	ENDIF

	A = (/alpha, AR, elongation, nu, theta, v/)
	LISTA = (/1,2,5,6/)

	X = dataFFT(3:15,1)
	Y = dataFFT(3:15,2)

	!average = SUM(Y) / NDATA
	!stdDev = SQRT(SUM((Y - average)**2)/NDATA)
	
	SIG = 0.1 * Y

	DO i=1,NDATA
		SIG(i) = SIG(i) / FLOAT(i)
	ENDDO

	chiSquared = (/1,2/)

	DO WHILE ( ABS(ALAMDA) .GT. 1.0E-04 .AND. ABS(ALAMDA) .LT. 1.0E08)
		CALL MRQMIN(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,COVAR,AMATRIX,NCA,CHISQ,modelWrapper,ALAMDA)
		PRINT*, "CHI, LAMBDA: ", CHISQ, ALAMDA
		PRINT*, "A: ", A

		chiSquared(2) = chiSquared(1)
		chiSquared(1) = CHISQ
	ENDDO

	WRITE(logString, '(A9,A1,F8.4,A1,F8.4,A1,F8.4,A1,F8.4,A1,F8.4,A1,F9.4)') inputFile(stringIdx:stringIdx + 8)," ", A(1)," ", A(2)," ", A(3)," ", A(4)/1000000," ", A(5)," ", A(6)/1000 
	CALL writeLog("polyFitLog.txt",TRIM(logString))

	model => makeModel(A(1),A(2),A(3),A(4),A(5),A(6),numRowsModel)
	WRITE(modelOutputFile, '(A6,A9,A5)') 'model_', inputFile(stringIdx:stringIdx + 8), '.data'			
	CALL writeData(modelOutputFile, model, numRowsModel)
END PROGRAM fit

SUBROUTINE modelWrapper(X,A,Y,DYDA,NA)
	USE marquardt_library
	USE make_model
	USE integrate_library

	REAL, DIMENSION(NA) :: A, tempA
	REAL, DIMENSION(13) :: X, Y, left, right
	REAL, DIMENSION(2) :: fRange, qyRange, zRange
	REAL, DIMENSION(13, NA) :: DYDA
	REAL :: fStep, qyStep, zStep, stepSize
	INTEGER :: numRows, NA

	fRange = (/ 0.2929688,1.464844 /)
	qyRange = (/ 1.0E-6,500.0E-6 /)
	zRange = (/ -2.0,0.0 /)
	
	fStep = 0.09765625
	qyStep = 10.0E-6
	zStep = 0.05

	numRows = getMaxIndex(fRange,fStep)

	X = getFrequencyList(fRange, fStep)
	Y = integrate(integrand, A, fRange, fStep, qyRange, qyStep, zRange, zStep)

	DO j=1,6
		IF (j .NE. 3 .AND. j .NE. 4) THEN	
			tempA = A
			tempA(j) = tempA(j) * 0.95

			left = integrate(integrand, tempA, fRange, fStep, qyRange, qyStep, zRange, zStep)

			tempA = A
			tempA(j) = tempA(j) * 1.05

			right = integrate(integrand, tempA, fRange, fStep, qyRange, qyStep, zRange, zStep)

			tempA = A
			stepSize = tempA(j) * 1.05 - tempA(j) * 0.95

			DYDA(:,j) = numDerivitive(left, right, 13, stepSize)
		ELSE
			DYDA(:,j) = 0
		ENDIF
	ENDDO
END SUBROUTINE modelWrapper

FUNCTION integrand(parameters, f, qy, z)
!Elongation should be input in degrees
!Theta should be input in arc seconds
!Computations done in SI units
!z should be input as multipliers of an AU (ex: z = (actual distance)/AU)
	IMPLICIT NONE

	REAL, DIMENSION(6) :: parameters
	REAL :: alpha, AR, elongation, nu, theta, v, f, qy, z, z0
	REAL :: au, PI, c, lambda, Vx, p, beta
	REAL(8) :: Fd, Fs, Ft, R, qx, qSquared
	REAL(8) :: integrand

	c = 3E8
	PI = 3.141592653589
	au = 149597870700
	
	alpha = parameters(1)
	AR = parameters(2)
	elongation = parameters(3) * (PI/180)
	nu = parameters(4)
	theta = parameters(5) * 4.84813681E-6
	v = parameters(6)
		 
	lambda = c / nu
	
	z = au * z
	z0 = z + au * COS(elongation)
	p = au * TAN(elongation) 
	R = SQRT(z**2 + p**2)

	Vx = v * p / R
	qx = 2 * PI * f / Vx

	qSquared =  (qx**2 + qy**2)

	Fd = SIN(qSquared * z0 * lambda / (4.0*PI))**2
	!Fs = EXP(-qSquared * (z0*theta/2.35)**2)
	Fs = EXP(-qSquared * (au * COS(elongation)*theta/2.35)**2)
	Ft = R**(-4) * (qx**2 + qy**2 / AR**2)**(-alpha/2)

	!Fd = SIN( ( (( 2.0*PI*f*  (p**2+z**2.0)**(0.5)  )/(v*p))**2  + qy**2 )* ((z+Cos(elongation)*au)*lambda)/(4.0*PI)   )**2
	!Fs = EXP(   -(  ( (( 2*PI*f*(p**2+z**2)**(0.5) )/(v*p))**2 + qy**2)**(0.5) * au * Cos(elongation) * (theta/2.35)  )**2.0 )
	!Ft = ( (( 2*PI*f*(p**2+z**2)**(0.5) )/(v*p))**2 + qy**2/AR**2)**(-alpha/2.0)

	integrand = Fd * Fs * Ft / Vx

	RETURN
END FUNCTION integrand