function get_grav, z, lat, lon, angle

; Paul Withers, 2001.07.03
; Open University, Britain

; Called by:				get_acc.pro, get_acc_model.pro
; Calls:				N/A
; Declares common blocks:		N/A
; Reads from:				N/A
; Writes to:				N/A
; Options:				N/A
; Purpose:				Return acceleration due to gravity in 
;					inertial cartesian frame given 
;					position in momentary spherical frame 
;					and time

common marsblk
; Properties of Mars and physical constants
; Defined in marsprops

; Include only C20 term beyond spherical symmetry

; Source of GM, Rref, and C20 is PDS dataset mors_1006
; Equation for GMM1 in Smith et al (1993)
; Equation for mors_1006 dataset.cat says it's in Tyler et al (1992)
; JGR v97, 7759-7779, MO RS instrument paper
; Should be the same, I will assume so!
; V = (gm)/r + (gm)/r * (rref/r)^2 [P2(cos(colat))] * c20 from Smith et al
; P2(x) = 0.5 (3*x^2 -1)

gm = marspropsstruct.gm	
;Grav constant*Mass of Mars in m^3/s^2
rref = marspropsstruct.rref
; Reference radius of Mars in m, 
; often mean equatorial radius
c20 = marspropsstruct.c20
; C20 dimensionless coefficient
; Normalization conventions are hideous, see Kaula book (1966)
; and Smith et al to understand. Also some Tolson/Withers emails.

r = z 
; Gives you option to use z as altitude, not radial distance
; r = z+rref

colat = double(90.) - lat
colat = colat*double(!pi/180.)
lon = lon * double(!pi/180.)
; Convert degrees to radians in momentary spherical frame

g_r = double(-1.)*gm/(r^2) - double(3.)*gm*(rref^2)/(r^4)*double(0.5)*(double(3.)*(cos(colat)^2)-double(1.))*c20
g_colat = gm*(rref^2)/(r^4) * double(0.5) * (double(-6.)*cos(colat)*sin(colat))*c20
g_lon = double(0.)
; Acceleration due to gravity in momentary spherical frame

g_intx = g_r*sin(colat)*cos(lon) + g_colat*cos(colat)*cos(lon) - g_lon*sin(lon)
g_inty = g_r*sin(colat)*sin(lon) + g_colat*cos(colat)*sin(lon) + g_lon*cos(lon)
g_intz = g_r*cos(colat)          - g_colat*sin(colat)
; This transformation from Bradbury, Theoretical Mechanics, p64
; It's a standard one
; Transform from momentary spherical frame to momentary cartesian frame

gx1 = g_intx * cos(angle) - g_inty * sin(angle)
gy1 = g_intx * sin(angle) + g_inty * cos(angle)
gz1 = g_intz
; Transform from momentary cartesian frame to inertial cartesian frame

return, [gx1,gy1,gz1]
end