%% plasipt11.m is designed to be called by plasipt1.m. It calculates
%% the oxygen and hydrogen densities for points above 350km. 

function [oxden,hyden,heden]=plasipt11(lonmat,latmat,dismat,timmat,day,month,year,lshcp,nop,iricf,gms,hes,sms)

% lonmat = geographic longitude of points of interest (rads)
% latmat = geographic latitude of points of interest (rads)
% dismat = radial distance of points of interest (km)
% timmat = univerisal time (fractional hours) (scalar)
% day = day of points of interest (scalar)
% month = month of points of interest (scalar)
% year = year of points of interest (scalar)
% lshcp = maximum L-shell for calculation of profiles (scalar)
% nop = number of points per field line (scalar)
% iricf = IRI2001 FoF2 data set to be used (0 = CCIR, other = URSI) (scalar)
% gms = include IRI2001 geomagnetic storm model (0 = no, 1 = yes) (scalar)
% hes = He+ simulated (0 = no, 1 = yes) (scalar)
% sms = switch to determine which version of MSIS to use (0 = MSIS90, 1 = NRLMSIS2000)(scalar)
%
% oxden = O+ density (O+/m^3)
% hyden = H+ density (H+/m^3)
% heden = He+ density (He+/m^3)

re=6371;                                          %Earth radi
dno=daynum(day,month,year);                       %day number
[lonmatm,latmatm]=geo2mag3(lonmat,latmat);        %magnetic co-ordinates
sv=ones(size(lonmat));        
lshmat=lshpar(-sv,latmatm,dismat);                %L-shell of points of interest
timm=timut2ml2(timmat*sv,lonmatm,dno*sv);         %Magnetic local time of points of interest

%Result matrices
oxden=zeros(size(lonmat));
hyden=oxden;
heden=oxden;

%% Load in plasimod data set for required day
[hytot,hetot,hyeqp,heeqp,oxeqp]=plasires3(timmat,dno,year,timm,lshmat);
if hes==0
 hetot=nan*ones(size(hytot));
 heeqp=nan*ones(size(hyeqp));
end

%% Calculate the required densities for inner plasmasphere.
poip=lshmat<=min(lshcp,7);                   %points inner plasmasphere
if any(poip)
 [oxden(poip),hyden(poip),heden(poip)]=plasipt12(hyeqp(poip),heeqp(poip),oxeqp(poip),lshmat(poip),lonmatm(poip),latmatm(poip),timmat,day,month,year,nop,iricf,gms,hes,sms);
end

%% Calculate the required densities for outer plasmasphere and polar regions.
if any(~poip)
 [oxden(~poip),hyden(~poip),heden(~poip)]=plasipt15(hyeqp(~poip),heeqp(~poip),lshmat(~poip),latmatm(~poip),dismat(~poip),lonmatm(~poip),timmat,year,month,day,lshcp,nop,iricf,gms,hes,sms);
end 
