MODULE make_model
CONTAINS
	FUNCTION makeModel(alpha, AR, elongation, nu, theta, v, numRowsOutput)
		USE integrate_library

		IMPLICIT NONE

		REAL :: alpha, AR, elongation, nu, theta, v, fStep, qyStep, zStep
		REAL, DIMENSION(2) :: fRange, qyRange, zRange
		REAL, DIMENSION(6) :: parameters
		REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: output
		REAL, DIMENSION(:,:), POINTER :: makeModel
		INTEGER :: fMaxIndex
		INTEGER, INTENT(OUT) :: numRowsOutput
	
		!alpha = 3.0
		!AR = 1.0
		!elongation = 20.0
		!nu = 140.0E6
		!theta = 0.0
		!v = 500.0E3

		parameters = (/ alpha, AR, elongation, nu, theta, v /)

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

		fMaxIndex = getMaxIndex(fRange,fStep)

		IF(ALLOCATED(output)) THEN		
			DEALLOCATE(output)			
		ENDIF

		ALLOCATE(output(fMaxIndex, 2))

		output(:,1) = getFrequencyList(fRange, fStep)

		output(:,2) = integrate(integrand, parameters, fRange, fStep, qyRange, qyStep, zRange, zStep)	

		makeModel => output

		numRowsOutput = fMaxIndex
	END FUNCTION makeModel

	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 = ABS(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
END MODULE make_model