%% gravmod.m calculates the gravity parallel to dipole field
%% line at the point of interest.
%%

function gv=gravmod(dismat,latmatm)

% dismat = radial distance of point of interest (km)
% latmatm = magnetic latitude of point of interest (rads)
% gv = gravity at point of interest (m/s^2)

re=6371;                                %Earth radi (km)
latmatm=abs(latmatm);                   %check latmatm positive
gv=9.8*(re./dismat).^2;                 %gravity at point of interest (m/s^2)
ltc=latmatm~=0;                         %non-zero latitude terms
gv(ltc)=2*gv(ltc).*tan(latmatm(ltc))./sqrt(1+4*tan(latmatm(ltc)).^2);      %gravity component parallel to field line (m/s^2)
gv(~ltc)=0;                                                                %gravity at equator
