pro gen_areoid

; Paul Withers, 2002.02.12
; LPL, Arizona

; Using Tolson's email of 2002.02.12, convert PDS altitudes into
; radial distances from centre of mass
; Then use these in detailed pressure integration...

n=181
m=361
; this is in areocentric lat/lon, PDS data are in areodetic...

lat = dblarr(n,m)
lon = lat
i=0
while i lt n do begin
j=0
while j lt m do begin
lat(i,j) = (2.*i - (n-1.))/(n-1.) ; -1 to 1
lon(i,j) = (2.*j - (m-1.))/(m-1.) ; -1 to 1
j=j+1
endwhile
i=i+1
endwhile

lat = lat*!pi/2.
lon = lon*!pi

;lat = double([-1., 0., 1.]*.45*!pi)
;lon = double([0.,0.,0.])

; coefficients from Tolson email
; pretty close to GMM-1 paper, Smith et al.
clm = dblarr(5,5)
clm(2,0) = -8759197.6
clm(2,1) = 132.
clm(2,2) = -843122.

clm(3,0) = -119340.7
clm(3,1) = 38656.8
clm(3,2) = -159257.9
clm(3,3) = 354141.3

clm(4,0) = 51505.2
clm(4,1) = 42392.1
clm(4,2) = -11164.7
clm(4,3) = 65141.9
clm(4,4) = 1130.

slm = clm*0.
slm(2,0) = 0.
slm(2,1) = 6.8
slm(2,2) = 496785.3

slm(3,0) = 0.
slm(3,1) = 252774.2
slm(3,2) = 84668.7
slm(3,3) = 251996.6

slm(4,0) = 0.
slm(4,1) = 37475.6
slm(4,2) = -89633.9
slm(4,3) = -2723.5
slm(4,4) = -128953.9

clm = clm*1e-10
slm = slm*1e-10

; values from Tolson email
; all seem reasonable
omega = double(7.0882e-5)	; rotation rate, rad s^-1
rref = double(3382.4) 		; reference radius, km, for clm, slm coeffs
gm = double(42828.2352)		; big-G * mass of Mars, km^3 s^-2
a = double(3393.4)		; IAU ellipsoid reference radius, km
f = double(0.0052083)		; IAU ellipsoid flattening, dimensionless

x = sin(lat)
y = (1-x^2)^0.5

; formulae from Tolson email and Arfken
p20 = 1.5*x^2 - 0.5
p21 = 3.*y*x
p22 = 3*y^2

p30 = 2.5*x^3 - 1.5*x
p31 = 7.5*y*(x^2-0.2)
p32 = 15.*y^2*x
p33 = 15.*y^3

p40 = (35*x^4 - 30*x^2+3)*0.125
p41 = 17.5*y*(x^3-3*x/7.)
p42 = 52.5*y^2*x^2-7.5*y^2
p43 = 105*y^3*x
p44 = 105*y^4   

; Normalization discussed in my notes of 2002.02.12
p20 = p20 * sqrt(5.)
p21 = p21 * sqrt(10./6.)
p22 = p22 * sqrt(10./24.)

p30 = p30 * sqrt(7.)
p31 = p31 * sqrt(14./12.)
p32 = p32 * sqrt(14./120.)
p33 = p33 * sqrt(14./720.)

p40 = p40 * sqrt(9.)
p41 = p41 * sqrt(18./20.)
p42 = p42 * sqrt(18./360.)
p43 = p43 * sqrt(18./5040.)
p44 = p44 * sqrt(18./40320.)

droverrref = $
 clm(2,0)*p20 + clm(3,0)*p30 + clm(4,0)*p40 + $
 clm(2,1)*p21*cos(lon) + clm(2,2)*p22*cos(2*lon) + $
 clm(3,1)*p31*cos(lon) + clm(3,2)*p32*cos(2*lon) + clm(3,3)*p33*cos(3*lon) + $
 clm(4,1)*p41*cos(lon) + clm(4,2)*p42*cos(2*lon) + clm(4,3)*p43*cos(3*lon) + $
  clm(4,4)*p44*cos(4*lon) + $
 slm(2,1)*p21*sin(lon) + slm(2,2)*p22*sin(2*lon) + $
 slm(3,1)*p31*sin(lon) + slm(3,2)*p32*sin(2*lon) + slm(3,3)*p33*sin(3*lon) + $
 slm(4,1)*p41*sin(lon) + slm(4,2)*p42*sin(2*lon) + slm(4,3)*p43*sin(3*lon) + $
;  slm(4,4)*p44 + $
  slm(4,4)*p44*sin(4*lon) + $
 0.5*omega^2*rref^3*(cos(lat))^2 / gm

; my first run had no sin(4*lon) on last term...
; found 2002.09.12    
equipot_radius = (droverrref + 1)*rref
ellipsoid_radius = a*(1-f)/sqrt(1-f*(2-f)*(cos(lat))^2)

deltah = equipot_radius - ellipsoid_radius
; PDS altitudes are measured from ellipsoid surface plus deltah
; so add PDS alts to (deltah + ellipsoid_radius) to get radius

openw, lun1, 'equipot_radius.dat', /get_lun
printf, lun1, equipot_radius, lat, lon
free_lun, lun1

stop
end