function gprime=gravity(phi,h)
% function gprime=gravity(phi,h)
% gprime is the acceleration of gravity in ft/sec*sec
% and is normal to the surface of the ellipsoid.
% Note: gprime=g-omegae X (omegae X p)) (Stevens and Lewis, p. 25
% phi is the geodetic latitude in radians
% h is the altitude above the ellipsoid in ft
% Ref: WGS-84 (NIMA TR8350.2, pages 4-1, 4-2)
% This approximation is good when h is small.

s2phi=(sin(phi))^2;
gammae=9.7803253359; %m/s^2 gravity at the equator on the ellipsoid
%gammap=9.8321849378; %m/s^2 gravity at the pole on the ellipsoid
k=.00193185265241;    % non-dim, % theoretical (normal) gravity formula constant
a=6378137.0;         % m, semi-major axis
%b=6356752.3142;      % m, semi-minor axis
%k=b*gammap/(a*gammae)-1 % theoretical (normal) gravity formula constant
e2=0.00669437999;    % e squared, e=first eccentricity of the ellipsoid
gamma=gammae*(1+k*s2phi)/sqrt(1-e2*s2phi);
%g=gamma/.3048 % convert g to ft/(sec*sec)
%f=1/298.257223563;
%m=.00344978650684;
%A=-2*(1+f+m)/a
A=-3.15704287057988e-07;
%B=4*f/a
B=2.10268965044023e-09;
%C=3/(a*a);
C=7.374516772942e-14;
H=h*.3048; % H (m)
gammah=gamma*(1 + (A + B*s2phi)*H + C*H*H);
gprime=gammah/.3048; % convert g to ft/(sec*sec)
