function ainteg1(zz)
C
	common		/params/ elong,as,be,NA, 	! parameters used by the_model to calculate qyI
     &                  zs,ze,ccc,PI,au,au2,anu,alamb,freq,twopi,qxf,
     &			L1,L2,L3,L4,L5,L6,L7,L8,L9,A1,A2,A3,A4,A5,A6,A7,A8,A9,
     &			z,z0,R,p,Vx,  	 		! all these parameters are calculated here to integrate for qy
     &			qy,X,
     &			er1,er2,er3,er4,er5,er6,er7,er8,er9,er10
C
	external	ainteg2	
C		
	print *, 'inside ainteg1_1', as,be,qy,z,zz
		z = zz					! z comes into program in AU

	        call qsimp2(ainteg2,as,be,qyI)		! as and be into program as inverse au 

		qy = qyI

C
		if(L1.eq.1) V		= A1
		if(L1.eq.2) V		= A2
		if(L1.eq.3) V		= A3
		if(L1.eq.4) V		= A4
		if(L2.eq.1) alpha	= A1
		if(L2.eq.2) alpha	= A2
		if(L2.eq.3) alpha	= A3
		if(L2.eq.4) alpha	= A4
		if(L3.eq.1) theta	= A1
		if(L3.eq.2) theta	= A2
		if(L3.eq.3) theta	= A3
		if(L3.eq.4) theta	= A4
		if(L4.eq.1) ar		= A1
		if(L4.eq.2) ar		= A2
		if(L4.eq.3) ar		= A3
		if(L4.eq.4) ar		= A4
C
	print *, ' inside ainteg1_2', L1, L2, L3, L4, V,alpha,theta,AR,qy,z
C
		ztoau = z*au				! z in m
C
		z0    = z + COS(elong)			! z0 in AU			
		z0toau = z0*au				! z0 in m
C
		aau = au
C		aau = atoau
C
C		p = au * TAN(elong) 
		p = TAN(elong)				! p in AU, this is the heliocentric distance of the "p" point

		R2 = (z**2 + p**2)			! R2 in AU^2
		R  = SQRT(R2)

		Vx = V * p / R				! V in AU/sec - Vx in AU/sec
C
		qx = qxf / Vx				! qxf in 1/sec --> qx in m if V in m/sec (qy is integrated from 15 to 1500 m)
C
		qxSquared = qx**2			! qx^2 in (1/m)^2
		qySquared = qy**2			! qy^2 in (1/m)^2
		qSquared =  (qxSquared + qySquared)	! q^2 - (1/m)^2
C
		Fd = SIN(qSquared * z0toau * alamb / (4.0*PI))**2
		expval	= -qSquared * (au * COS(elong)*theta/2.35)**2
		if(expval.lt.-50.) expval = -50.0				! prevent underflow for the exponential function
		Fs = EXP(expval)
		Ft = R2**(-2) * ((qxSquared + qySquared / AR**2)*aau)**(-alpha/2.0) 
C
		ainteg1 = (Fd * Fs * Ft)*aau/Vx
C
	print *, ' At end of ainteg1_2', V,alpha,theta,AR, z, qx,qy,Fd,Fs,Ft,ainteg1

	return
	end		
C
C
C**********************
C
	function ainteg
	real		A	(NA)
C
	common		/params/ elong,as,be,NA, 	! parameters used by the_model to calculate qyI
     &                  zs,ze,ccc,PI,au,au2,anu,alamb,freq,twopi,qxf,
     &			L1,L2,L3,L4,L5,L6,L7,L8,L9,A1,A2,A3,A4,A5,A6,A7,A8,A9,
     &			z,z0,R,p,Vx,  	 		! some of these parameters are used to integrate qy
     &			qyI,X,				! this parameter determined here
     &			er1,er2,er3,er4,er5,er6,er7,er8,er9,er10
C
C		zs 	= -2.0
		ze	= 0.0
		c 	= 3.0E8
		PI 	= 3.141592653589
		au 	= 149597870700.
		anu 	= 140.0E6
	print *, 'Into ainteg1', X, A, NA, elong, as, be
		elongation = elong * (PI/180.0)
C
		do I = 1, NA
C		  if(lista(i).eq.1) V     = A(lista(i))
C		  if(lista(i).eq.2) alpha = A(Lista(i))
C		  if(lista(i).eq.3) theta = A(Lista(i))
C		  if(lista(i).eq.4) AR    = A(Lista(i))
		end do
C		if(Lista(1).eq.1) 			
		V 	= A(1)				! Velocity input in kilometers/sec
		alpha 	= A(2)				! integrated turbulence
		theta 	= A(3) *  4.84813681E-6		! Theta input in arc seconds
		AR	= A(4)				! Axial ratio
C		 
		alambda = c / anu
C	
		z = au * z

		z0 = z + au * COS(elongation)
		p = au * TAN(elongation) 
		R = SQRT(z**2 + p**2)
C
		Vx = v * p / R
C
		qx = 2.0 * PI * X / Vx
C
		qSquared =  (qx**2 + qy**2)

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

		ainteg1 = Fd * Fs * Ft / Vx
	print *, 'At end of ainteg1', X, ainteg1, Fd, Fs, Ft, Vx, p, v, R, z
	return
        end
C
C**************************************************************************************************
C