%% plasisc11.m is designed to be called by plasisc1.m. It calculates
%% the electron density for points above 350km. 
%% The full model is used for the regions below min(7,lshcp), while for all
%% points above this are assumed to be 'polar' including points on
%% field lines in the bulge region.
%%
%% 16/11/2001 - updated to include He+.
%% 5/3/2002 - remove different files name extension for IRI CCIR or URSI parameters (store iricf in arrays).
%% 5/3/2002 - He+ inclusion now optional. 
%%

function [eden,oxden,hyden,heden]=plasisc11(lonlsh,latmv,dismatv,timmat,day,month,year,nop,iricf,gms,hes,sms)

% lonlsh = magnetic longitude of slice (rads) (scalar)
% latmv = magnetic latitude vector of grid (rads)
% dismatmv = radial distance vector of grid (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)
% 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)
%
% eden = electron density (e-/m^3)
% oxden = O+ density (O+/m^3)
% hyden = H+ density (H+/m^3)
% heden = He+ density (He+/m^3)

re=6371;                                          %Earth radi
r0=re+350;                                        %base radial distance
dno=daynum(day,month,year);                       %day number
timm=timut2ml2(timmat,lonlsh,dno);                %magnetic local time of slice (MLT)
dismv=dismatv;                                   
dismv(dismv<=(re+350))=[];                        %radial vector above 350km altitude
latmv=latmv(:);

[lshcp,lshmx]=plasires2(dno,year);

nopx=length(dismv);                               %number of radial points
nopy=length(latmv);                               %number of latitudinal points
svm=ones(nopy,nopx);

%Result matrices               
oxden=0*svm;
hyden=oxden;
heden=oxden;

[dismat,latmatm]=meshgrid(dismv,latmv);           %radial distance and magnetic latitude matrices of grid
[xmat,ymat]=pol2cart(latmatm,dismat);             %cartesian co-ordinates of grid
xmat(xmat<1e-3)=0;                                %correction for round off error
[lonmat,latmat]=mag2geo3(lonlsh*svm,latmatm);     %geographic co-ordinates of grid

%% Calculate the L-shell of each point in the grid
latmv2=latmv;                                     %copy of latmv to calculate L-shell          
rd89=deg2rad(89);
latmv2(latmv2>rd89)=rd89;                         %reset values above 89 degrees to 89 degrees
latmv2(latmv2<-rd89)=-rd89;                       %reset values below -89 degrees to -89 degrees
lshm=lshpar(-svm,latmv2*ones(1,nopx),dismat);     %L-shell of grid points

%% Determine in what region each point on the grid lies. These are the plasmasphere
%% inside L=min(7,lshcp) (northern and southern hemispheres) and outside (northern
%% and southern hemispheres).
poi1=lshm<min(7,lshcp);                           %points inside plasmasphere
poi11=poi1 & latmatm>=0;                          %northern plasmasphere
poi12=poi1 & latmatm<0;                           %southern plasmasphere
poi2=~poi1;                                       %'polar regions'
poi21=poi2 & latmatm>=0;                          %northern polar regions
poi22=poi2 & latmatm<0;                           %southern polar regions

%% Construct an L-shell vector at which to determine the H+ and He+ equatorial densities
%% including points outside the plasmapause. Therefore need to replace the -1 terms
%% in the look up data tables with the density calculated from tube content divided
%% by tube volume.
lshps1=[375:25:600 650:50:5000];  
lshps=[(re+lshps1)/re 1.8:0.05:2 2.1:0.1:6 6.2:0.2:lshmx];

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

%% Calculate the ion densities along each of the L-shells inside L=min(7,lshcp)
lshps2=lshps;                           %copy of lshps
hyeqp2=hyeqp;                           %copy of hyeqp
heeqp2=heeqp;                           %copy of heeqp
oxeqp2=oxeqp;                           %copy of oxeqp
lshps2(lshps>lshcp)=[];                 %remove terms beyond L=lshcp
hyeqp2(lshps>lshcp)=[];
oxeqp2(lshps>lshcp)=[];
if hes==1
 heeqp2(lshps>lshcp)=[];
end

[hydenf,hedenf,oxdenf,hydenfc,hedenfc,oxdenfc,latmatmf,dismatf]=plasisc12(hyeqp2,heeqp2,oxeqp2,lshps2,lonlsh,timmat,year,month,day,nop,iricf,gms,hes,sms);

%% Use the MATLAB griddata command to calculate the ion 
%% densities at the points on the grid inside L=min(7,lshcp).
[xmatn,ymatn]=pol2cart(latmatmf,dismatf);                             %northern hemisphere field lines cartesian co-ords
[xmats,ymats]=pol2cart(-latmatmf,dismatf);                            %southern hemisphere field lines cartesian co-ords
xmatn(xmatn<1e-3)=0;                                                  %correction for round off error
xmats(xmats<1e-3)=0;                                                  %correction for round off error
oxdpn=griddata(xmatn,ymatn,log10(oxdenf),xmat(poi11),ymat(poi11));    %northern O+ densities
oxdps=griddata(xmats,ymats,log10(oxdenfc),xmat(poi12),ymat(poi12));   %southern O+ densities
hydpn=griddata(xmatn,ymatn,log10(hydenf),xmat(poi11),ymat(poi11));    %northern H+ densities
hydps=griddata(xmats,ymats,log10(hydenfc),xmat(poi12),ymat(poi12));   %southern H+ densities
oxden(poi11)=10.^oxdpn;                                               %record northern O+ densities 
oxden(poi12)=10.^oxdps;                                               %record southern O+ densities 
hyden(poi11)=10.^hydpn;                                               %record northern H+ densities 
hyden(poi12)=10.^hydps;                                               %record southern H+ densities 
if hes==1
 hedpn=griddata(xmatn,ymatn,log10(hedenf),xmat(poi11),ymat(poi11));   %northern He+ densities
 hedps=griddata(xmats,ymats,log10(hedenfc),xmat(poi12),ymat(poi12));  %southern He+ densities
 heden(poi11)=10.^hedpn;                                              %record northern He+ densities 
 heden(poi12)=10.^hedps;                                              %record southern He+ densities 
end

%% To determine the ion densities in the polar regions need to determine
%% the field lines along which to determine profiles. As L-shell values 'blow
%% out' in the polar regions, do this using magnetic latitude at 350km starting
%% from inside lshcp up to required pole ie -+90 degrees in 1 degree steps.
latmpol=round(rad2deg(lshpar((lshcp-0.5),-1,r0)));                    %magnetic latitude to start polar region
latmpvd=[latmpol:80 80.5:0.5:90];                                     %polar base magnetic latitude vector (degrees)
latmpv=deg2rad(latmpvd);                                              %polar base magnetic latitude vector
[oxdpoln,hydpoln,hedpoln,latmpoln,dispoln]=plasisc13(hyeqp,heeqp,lshps,latmpv,lshcp,lonlsh,timmat,year,month,day,nop,iricf,gms,hes,sms);
[oxdpols,hydpols,hedpols,latmpols,dispols]=plasisc13(hyeqp,heeqp,lshps,-latmpv,lshcp,lonlsh,timmat,year,month,day,nop,iricf,gms,hes,sms);
[xmatnp,ymatnp]=pol2cart(latmpoln,dispoln);                           %northern hemisphere field lines cartesian co-ords
[xmatsp,ymatsp]=pol2cart(latmpols,dispols);                           %southern hemisphere field lines cartesian co-ords
xmatnp(xmatnp<1e-3)=0;                                                %correction for round off error
xmatsp(xmatsp<1e-3)=0;                                                %correction for round off error
oxdpn=griddata(xmatnp,ymatnp,log10(oxdpoln),xmat(poi21),ymat(poi21)); %northern polar O+ densities
oxdps=griddata(xmatsp,ymatsp,log10(oxdpols),xmat(poi22),ymat(poi22)); %southern polar O+ densities
hydpn=griddata(xmatnp,ymatnp,log10(hydpoln),xmat(poi21),ymat(poi21)); %northern polar H+ densities
hydps=griddata(xmatsp,ymatsp,log10(hydpols),xmat(poi22),ymat(poi22)); %southern polar H+ densities
oxden(poi21)=10.^oxdpn;                                               %record northern polar O+ densities 
oxden(poi22)=10.^oxdps;                                               %record southern polar O+ densities 
hyden(poi21)=10.^hydpn;                                               %record northern polar H+ densities 
hyden(poi22)=10.^hydps;                                               %record southern polar H+ densities 
if hes==1
 hedpn=griddata(xmatnp,ymatnp,log10(hedpoln),xmat(poi21),ymat(poi21));     %northern polar He+ densities
 hedps=griddata(xmatsp,ymatsp,log10(hedpols),xmat(poi22),ymat(poi22));     %southern polar He+ densities
 heden(poi21)=10.^hedpn;                                                   %record northern polar He+ densities 
 heden(poi22)=10.^hedps;                                                   %record southern polar He+ densities 
end

%Check for any NaN's and use interpolation to fill them in
pnan=isnan(oxden);
if any(any(pnan))
 xmn=xmat;
 ymn=ymat;
 oxn=oxden;
 hyn=hyden;
 xmn(pnan)=[];
 ymn(pnan)=[];
 oxn(pnan)=[];
 hyn(pnan)=[];
 oxl=griddata(xmn,ymn,oxn,xmat(pnan),ymat(pnan));
 hyl=griddata(xmn,ymn,hyn,xmat(pnan),ymat(pnan));
 oxden(pnan)=oxl;
 hyden(pnan)=hyl;
 if hes==1
  hen=heden;
  hen(pnan)=[];
  hel=griddata(xmn,ymn,hen,xmat(pnan),ymat(pnan));
  heden(pnan)=hel;
 end
end 

%Check for any remaining NaN's in last column. Replace with 
%nearest inner (radially distance) value.
pnan2=isnan(oxden(:,nopx));
if any(pnan2)
 pnan3=find(pnan2);
 oxden(pnan3,nopx)=oxden(pnan3,nopx-1);
 hyden(pnan3,nopx)=hyden(pnan3,nopx-1);
 if hes==1
  heden(pnan3,nopx)=heden(pnan3,nopx-1);
 end
end

pnan=isnan(oxden);
if any(any(pnan))
 disp('Phil - some NaN still in plasisc11.m')
 keyboard
end 

%Calculate electron densities
if hes==1
 eden=oxden+hyden+heden;
else
 eden=oxden+hyden;
end

if hes==0
 heden=nan*ones(size(hyden));
end

