	program formula
c
	implicit none

	real*8 DSIND
	real*8 DCOSD

	real*8	A(3),freq,freq1,c,wavelen,arg,Ymod,Powr(24,40),sum
	real*8 freqsv(24),au,toyofreq,z00,p,R,v,pi4
	real*8 Fd,Fs,Ft,Vx,qSquared,qx,qxSquared,qy,q,z,z0,elongation
	integer*4 i,j,k,m
c
	data Powr /960*0.D0 /	
	data elongation /23.3D0/	! 78.6762 for 7 Feb 2013
	c	= 2.99792458D+8		! Speed of light, m/sec
	au	= 1.4959787066D+11	! Astronomical unit, meters
	pi4     = 12.56637061D0		! 4 pi
	toyofreq= 327.E+6
	wavelen = c/toyofreq
	print 100, wavelen
100	format(' Wavelength of array (m) ',f10.3)
C
C 	Julio's formula
C
c`	A(1) is alpha	!1 to 3.5
c	A(2) is theta	!.3 arc-sec  or 0.1
c	A(3) is speed V
c	AR assume for now to be 1.0
c	Here also we assume Sun is exactly 1 AU from Earth
c
c	A(1)=2.0D0
	A(1)=3.0D0			!Alpha
	A(2)=0.1D0			!source size in arc-sec
c	A(2)=A(2)*4.84813681D-6		!convert to radians, factor of 2 error...
	A(2)=A(2)*9.696273622D-6	!convert to radians
	p =au*DSIND(elongation)	
	z0=au*DCOSD(elongation)		!p and z0 in meters
	do i=1,40			!speed loop 100 --> 2050 km/sec, in steps of 50 km/sec
	v=(1.D+5)+(5.D+4)*dfloat(i-1)	!v in m/sec
	A(3)=v		!/au		!A(3) in m/sec
	do k=1,24			!frequency loop, 0.071 to 14.125 Hz
	freq=0.05D0+0.1D0*dfloat(k-3)
	freq1=10**(freq-1.0D0)
	freqsv(k)=freq1
c	
c	print *,k,freq,freq1,wavelen
c
	sum=0.D0
	do j=1,500			!qy loop, 1E-6 --> 500E-6, steps of 1E-6
	qy=1.D-6*dfloat(j)		!qy in 1/m
c	qy=qy/au
c
c	Here we flag a possible error in Julio's formula for "z" (substituting "-" in place of "+"),
c	if "z" is the distance along the line of sight between closest approach to the Sun "p",
c	and the right triangle having Earth-Sun line (au=1) as its base and "R" as its height.
c	On the other hand, integrating over z, with the indicated z range "-2 to 0", this formula 
c	looks OK: here the step size is 0.005. 
c	Also "qy" *is* integrated over, the out-of-the-plane dimension.
c
	do m=1,400			!z loop, from -1 to +1 au
	z=0.005D0*dfloat(m-201)
	z=au*z				!z in meters
	z00=z0+z
	R=DSQRT(z**2+p**2)		!R in meters
	Vx=A(3)*p/R			!Vx in m/sec
	qx=freq1/Vx			!qx in 1/m
	qxSquared=qx*qx
	qSquared=qxSquared+qy*qy
	q=DSQRT(qSquared)		!This statement could be omitted by using arg=-A(1)/2.D0 below
!	print *,m,qSquared,z,wavelen
	Fd=(DSIN(qSquared*z00*wavelen/pi4))**2	!As in Bernie's program
c	Fd=(DSIN(qSquared*z*wavelen/pi4))**2	!As in Julio's thesis
c
c	Here Bernie's treatment differs from Julio's formula in that "z0 here" is used instead of "z00 here".
c	This change also appears in "fit_fft2.f".
c	
	Fs=DEXP(-((qSquared*z*A(2)/2.36D0)**2))
	arg=-A(1)			!/2.D0: here we don't divide by two since we're using q directly
	R=R/au				!Hipshot: maybe R, only just here, should be in au???
	Ft=(R**(-3))*(q**arg)		!This is what Bernie thinks ...
c	Ft=(q**arg)			!No R at all ...
	Ymod=Fd*Fs*Ft/Vx
c	print *,m,Fd,Fs,Ft,Vx,R
	sum=sum+Ymod
	enddo	!z loop, m
	enddo	!qy loop, j
	Powr(k,i)=sum
	enddo	!frequency loop, k
	enddo	!speed loop, i
c
c	renormalize to an arbitrary factor of 10000 lower, needed for alpha = 3, but not for = 2.
c
	do k=1,24
	do i=1,40
		Powr(k,i)=Powr(k,i)/10000.0D0
	enddo
	enddo	
	!open (30,file='C:\Marquardt\JulioFormula.txt')
	open (30,file='JulioFormula.txt')
	do k=1,24
	write (30,30)k,freqsv(k),(Powr(k,i),i=1,40)
30	format(i3,f7.3,40f16.0)	
	enddo
	close (30)
C
	stop
	end

