%% plasigb11.m is designed to be called by plasigb1.m. It is designed to
%% calculate the electron and ion densities above 350km, from which to 
%% determine the global TEC maps. 
%%
%% 17/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]=plasigb11(lonmat,latmat,disv,timmat,day,month,year,nop,iricf,gms,hes,sms)

% lonmat = geographic longitude matrix of grid (rads)
% latmat = geographic latitude matrix of grid (rads)
% disv = radial distance vector of grid (km)
% timmat = univerisal time (fractional hours) (scalar)
% day = day (scalar)
% month = month (scalar)
% year = year (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)
% sms = switch to determine which version of MSIS to use (0 = MSIS90, 1 = NRLMSIS2000)(scalar)
%
% eden = e- 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
lonmv=deg2rad(-180:15:165);                       %magnetic longitude vector at which to calculate 'magnetic slices' (rads)     
latmv=deg2rad([-90:5:-35 -30:2.5:30 35:5:90]);    %magnetic latitude vector at which to use griddate on 'magnetic slices' (rads)
nopx=length(lonmv);
nopy=length(latmv);
nopz=length(disv);
[lonmg,latmg]=geo2mag3(lonmat,latmat);            %magnetic co-ords of geographic grid
dno=daynum(day,month,year);                       %day number
[lshcp,lshmx]=plasires2(dno,year);
kp=solarval(timmat,dno,year);                     %Kp index

%% Set up results matrices for the ion densities. Create one O+ and H+ 2D array 
%% for each required longitude. (Trying to do the indexing in 3D was giving me a
%% headache!)
resm=ones(nopy,nopz);
for i=1:nopx
 eval(['oxden',int2str(i),'=resm;'])
 eval(['hyden',int2str(i),'=resm;'])
 eval(['heden',int2str(i),'=resm;'])
end 

%% Calculate the L-shell of each point in the region bound by latmv by disv.
latmatm=latmv'*ones(1,nopz);                                %magnetic latitude marix
nopz=length(disv);
lshm=lshpar(-ones(nopy,nopz),latmatm,ones(nopy,1)*disv);    %L-shell of grid points
lshl=(re+375)/re;                                           %minimum allowed L-shell
lshm(lshm<lshl)=lshl;                                       %reset L-shells below lshl to lshl

%% Determine in what region each point in the region bound by latmv by disv. These
%% are the plasmasphere inside L=lshcp (northern and southern hemispheres) and outside
%% L=lshcp (northern and southern hemispheres).
poi1=lshm<=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. 
timm=timut2ml2(timmat*ones(1,nopx),lonmv,dno*ones(1,nopx));           %MLT
lshps1=[375:25:500 600:50:1000 1100:100:2500 2750:250:5000];  
lshps=[(re+lshps1)/re 1.8:0.05:2 2.1:0.1:6 6.2:0.2:lshmx];            %L-shells at which to calculate field lines     
latmb=lshpar(lshps,-ones(size(lshps)),r0*ones(size(lshps)));          %magnetic latitudes of field lines

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

%% Determine the cartesian co-ordinates of the region bound by latmv by disv.
[xmat,ymat]=pol2cart(latmatm,ones(nopy,1)*disv);

%% Calculate the ion densities along each of the L-shells inside L=lshcp
lshps2=lshps;                           %copy of lshps
hyeqp2=hyeqp;                           %copy of hyeqp
heeqp2=heeqp;                           %copy of heeqp
oxeqp2=oxeqp;                           %copy of heeqp
lshps2(lshps>lshcp)=[];                 %remove terms beyond L=lshcp
hyeqp2(:,lshps>lshcp)=[];
if hes==1
 heeqp2(:,lshps>lshcp)=[];
else
 heeqp2=nan*ones(size(hyeqp2));
end
oxeqp2(:,lshps>lshcp)=[];
lpl=length(lshps2);
lonm2=lonmv'*ones(1,lpl);
lshps2=ones(nopx,1)*lshps2;

%% Break the fitting into two parts, so as not to overload the memory. About 1500 
%% field lines is all it can handle without having to start writing to the hard drive.
lplh=ceil(lpl/2);
ind1=1:(lplh);                       %1st index terms
ind2=(lplh+1):lpl;                   %2nd index terms
[hydenpsn1,hedenpsn1,oxdenpsn1,hydenpss1,hedenpss1,oxdenpss1,latmatmps1,dismatps1]=plasisc12(hyeqp2(:,ind1),heeqp2(:,ind1),oxeqp2(:,ind1),lshps2(:,ind1),lonm2(:,ind1),timmat,year,month,day,nop,iricf,gms,hes,sms);
[hydenpsn2,hedenpsn2,oxdenpsn2,hydenpss2,hedenpss2,oxdenpss2,latmatmps2,dismatps2]=plasisc12(hyeqp2(:,ind2),heeqp2(:,ind2),oxeqp2(:,ind2),lshps2(:,ind2),lonm2(:,ind2),timmat,year,month,day,nop,iricf,gms,hes,sms);
hydenpsn=[hydenpsn1;hydenpsn2];
oxdenpsn=[oxdenpsn1;oxdenpsn2];
hydenpss=[hydenpss1;hydenpss2];
oxdenpss=[oxdenpss1;oxdenpss2];
if hes==1
 hedenpsn=[hedenpsn1;hedenpsn2];
 hedenpss=[hedenpss1;hedenpss2];
end
latmatmps=[latmatmps1;latmatmps2];
dismatps=[dismatps1;dismatps2];

clear hydenpsn1 hedenpsn1 oxdenpsn1 hydenpss1 hedenpss1 oxdenpss1 latmatmps1 dismatps1
clear hydenpsn2 hedenpsn2 oxdenpsn2 hydenpss2 hedenpss2 oxdenpss2 latmatmps2 dismatps2

for i=1:nopx                                         %loop through each magnetic longitude
 ix=i:nopx:(i+nopx*(size(lonm2,2)-1));               %index of i'th rows
 if i==1
  [xmatps,ymatps]=pol2cart(latmatmps(ix,:),dismatps(ix,:));                     %field lines cartesian co-ords
 end 
 xmatps(xmatps<1e-3)=0;                                                         %correction for round off error

 oxdp=griddata(xmatps,ymatps,log10(oxdenpsn(ix,:)),xmat(poi11),ymat(poi11));    %northern O+ densities
 hydp=griddata(xmatps,ymatps,log10(hydenpsn(ix,:)),xmat(poi11),ymat(poi11));    %northern H+ densities
 if hes==1
  hedp=griddata(xmatps,ymatps,log10(hedenpsn(ix,:)),xmat(poi11),ymat(poi11));   %northern He+ densities
 end
 eval(['oxden',int2str(i),'(poi11)=10.^oxdp;']);                                %record northern O+ densities for i'th longitude
 eval(['hyden',int2str(i),'(poi11)=10.^hydp;']);                                %record northern H+ densities for i'th longitude
 if hes==1
  eval(['heden',int2str(i),'(poi11)=10.^hedp;']);                               %record northern He+ densities for i'th longitude
 end

 oxdp=griddata(xmatps,-ymatps,log10(oxdenpss(ix,:)),xmat(poi12),ymat(poi12));   %southern O+ densities
 hydp=griddata(xmatps,-ymatps,log10(hydenpss(ix,:)),xmat(poi12),ymat(poi12));   %southern H+ densities
 if hes==1
  hedp=griddata(xmatps,-ymatps,log10(hedenpss(ix,:)),xmat(poi12),ymat(poi12));  %southern He+ densities
 end
 eval(['oxden',int2str(i),'(poi12)=10.^oxdp;']);                                %record southern O+ densities for i'th longitude
 eval(['hyden',int2str(i),'(poi12)=10.^hydp;']);                                %record southern H+ densities for i'th longitude
 if hes==1
  eval(['heden',int2str(i),'(poi12)=10.^hedp;']);                               %record southern He+ densities for i'th longitude
 end

 %Check for any NaN's in last column. Replace with nearest inner (radially distance) value.
 eval(['pnan2=isnan(oxden',int2str(i),'(:,nopz));'])
 if any(pnan2)
  pnan3=find(pnan2);
  eval(['oxden',int2str(i),'(pnan3,nopz)=oxden',int2str(i),'(pnan3,nopz-1);'])
  eval(['hyden',int2str(i),'(pnan3,nopz)=hyden',int2str(i),'(pnan3,nopz-1);'])
  if hes==1
   eval(['heden',int2str(i),'(pnan3,nopz)=heden',int2str(i),'(pnan3,nopz-1);'])
  end
 end
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.1),-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
lshpol=lshpar(-ones(size(latmpv)),latmpv,r0*ones(size(latmpv)));      %L-shell of polar field lines

hyeqp2=interp1(lshps',hyeqp',lshpol');                                %required hyeqp values in polar regions
hyeqp2=hyeqp2';                                                       %return hyeqp2 to required size
hynan=isnan(hyeqp2);                                                  %check for NaN
lshpolm=ones(size(hynan,1),1)*lshpol;
[hybk,hebk]=bgden(lshpolm(hynan),kp);                                 %calculate NaN for background values
hyeqp2(hynan)=hybk;                                                   %reset NaN to background value
if hes==1
 heeqp2=interp1(lshps',heeqp',lshpol');                               %required heeqp values in polar regions
 heeqp2=heeqp2';                                                      %return heeqp2 to required size
 heeqp2(hynan)=hebk;                                                  %reset NaN to background value
end

latm=ones(nopx,1)*latmpv;                                             %required magnetic latitude values
lonm=lonmv(:)*ones(1,length(latmpv));                                 %required magnetic longitude values 
[oxdenpon,hydenpon,hedenpon,latmatmpo,dismatpo]=plasirg13(hyeqp2(:),heeqp2(:),lonm(:),latm(:),timmat,year,month,day,lshcp,nop,iricf,gms,hes,sms);
[oxdenpos,hydenpos,hedenpos]=plasirg13(hyeqp2(:),heeqp2(:),lonm(:),-latm(:),timmat,year,month,day,lshcp,nop,iricf,gms,hes,sms);
for i=1:nopx                                                          %loop through each magnetic longitude
 ix=i:nopx:(i+nopx*(size(lonm,2)-1));                                 %index of i'th rows
 if i==1
  [xmatpon,ymatpon]=pol2cart(latmatmpo(ix,:),dismatpo(ix,:));                        %northern hemisphere field lines cartesian co-ords
 end 
 xmatpon(xmatpon<1e-3)=0;                                                            %correction for round off error
 oxdpn=griddata(xmatpon,ymatpon,log10(oxdenpon(ix,:)),xmat(poi21),ymat(poi21));      %northern O+ densities
 hydpn=griddata(xmatpon,ymatpon,log10(hydenpon(ix,:)),xmat(poi21),ymat(poi21));      %northern H+ densities
 oxdps=griddata(xmatpon,-ymatpon,log10(oxdenpos(ix,:)),xmat(poi22),ymat(poi22));     %southern O+ densities
 hydps=griddata(xmatpon,-ymatpon,log10(hydenpos(ix,:)),xmat(poi22),ymat(poi22));     %southern H+ densities
 eval(['oxden',int2str(i),'(poi21)=10.^oxdpn;']);                                    %record northern O+ densities for i'th longitude
 eval(['hyden',int2str(i),'(poi21)=10.^hydpn;']);                                    %record northern H+ densities for i'th longitude
 eval(['oxden',int2str(i),'(poi22)=10.^oxdps;']);                                    %record southern O+ densities for i'th longitude
 eval(['hyden',int2str(i),'(poi22)=10.^hydps;']);                                    %record southern H+ densities for i'th longitude
 if hes==1
  hedpn=griddata(xmatpon,ymatpon,log10(hedenpon(ix,:)),xmat(poi21),ymat(poi21));     %northern He+ densities
  hedps=griddata(xmatpon,-ymatpon,log10(hedenpos(ix,:)),xmat(poi22),ymat(poi22));    %southern He+ densities
  eval(['heden',int2str(i),'(poi21)=10.^hedpn;']);                                   %record northern He+ densities for i'th longitude
  eval(['heden',int2str(i),'(poi22)=10.^hedps;']);                                   %record southern He+ densities for i'th longitude
 end
end

%% For NaNs being produced at the magnetic equator at very low altitudes (between
%% 350km and 375km). Replace with average of adjacent values.
for i=1:nopx
 eval(['poi=find(isnan(oxden',int2str(i),'));'])
 if any(poi)
  eval(['oxden',int2str(i),'(poi)=0.5*(oxden',int2str(i),'(poi-1)+oxden',int2str(i),'(poi+1));'])
  eval(['hyden',int2str(i),'(poi)=0.5*(hyden',int2str(i),'(poi-1)+hyden',int2str(i),'(poi+1));'])
  if hes==1
   eval(['heden',int2str(i),'(poi)=0.5*(heden',int2str(i),'(poi-1)+heden',int2str(i),'(poi+1));'])
  end
 end 
end 
 
%% Combine the results for each longitude slice into a single 3D matrix. 
oxdent=zeros(nopy,(nopx+1),nopz);
hydent=oxdent;
hedent=oxdent;
for i=1:nopx
 eval(['oxdent(:,i,:)=oxden',int2str(i),';'])
 eval(['hydent(:,i,:)=hyden',int2str(i),';'])
 if hes==1
  eval(['hedent(:,i,:)=heden',int2str(i),';'])
 end
end 

%% Copy -pi result to pi for lookup arrays.
oxdent(:,(nopx+1),:)=oxdent(:,1,:);
hydent(:,(nopx+1),:)=hydent(:,1,:);
if hes==1
 hedent(:,(nopx+1),:)=hedent(:,1,:);
end 
lonmv=[lonmv pi];

%% Use linear interpolation of each disv slice to determine the ion densities
%% at the required points of interest.
oxden=pwint2m(lonmv,latmv,oxdent,lonmg,latmg);
hyden=pwint2m(lonmv,latmv,hydent,lonmg,latmg);
if hes==1
 heden=pwint2m(lonmv,latmv,hedent,lonmg,latmg);
end

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

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