%% Modgall.m is the Gallagher et al. emprical plasmasphere model as
%% outlined in Advances in Space Research 1988 no. 8.  Note that the
%% model is based on satellite data from DE1 for 1981 and 1982 only.
%% It is a 'quiet period' model and does not take into account variation
%% with solar cycle etc.
%% The only inputs that it requires are L shell, altitude, magnetic 
%% longitude and univerisal time.
%%

function plsd=modgall(lsh,alt,timmat,lonmatm,dno)

% lsh = L-shell of the point of interest
% alt = altitude of the point of interest (km)
% timmat = univerisal time of point of interest (fractional hours)
% lonmatm = magnetic longitude of point of interest (rads)
% dno = day number
% plsd = plasma (H+) density (cm^-3)

mlt=timut2ml2(timmat,lonmatm,dno);                      %calculate local magnetic time
x=mlt;                                                  %x parameter required by model
x(x>12)=24-x(x>12);                                     %correct value for x in region 12 to 24 hours

%Constant parameters in model
a1=1.4;
a2=1.53;
a3=-0.036;
a4=30.76;
a5=159.9;
a7=6.27;

%Variable parameters in model
a6=-0.87+0.12*exp(-x.^2/9);
a8=0.7*cos(2*pi*(mlt-21)/24)+4.4;
a9=15.3*cos(2*pi*mlt/24)+19.7;

%Main three terms in model
fl=a2-exp(a3*(1-a4*exp(-alt/a5)));
gl=a6.*lsh+a7;
hl=(1+(lsh./a8).^(2*(a9-1))).^(-(a9./(a9-1)));

%Calculate plasma density
plsd=10.^(a1*fl.*gl.*hl);