%% plasirg11.m is designed to be called by plasirg1.m. It is designed to
%% calculate the electron, oxygen and hydrogen ion densities above 350km,
%% from which to determine the regional 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]=plasirg11(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)
% hes = He+ simulated (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:10:180);                       %magnetic longitude vector at which to calculate 'magnetic slices' (rads)     
lonmvb=deg2rad(0:10:360);                         %2nd version of lonmv
latmv=deg2rad([-90:5:-35 -30:2.5:30 35:5:90]);    %magnetic latitude vector at which to use griddata on 'magnetic slices' (rads)
dno=daynum(day,month,year);                       %day number
[lshcp,lshmx]=plasires2(dno,year);
kp=solarval(timmat,dno,year);                     %Kp index

%% Determine the magnetic longitudes of the boundary region required for which
%% to calculate the required ion densities. Note that because of 'wrap around' 
%% occuring with the longitude values, need to be a little tricky to select out
%% the correct lonmv2 vector.
[lonmg,latmg]=geo2mag3(lonmat,latmat);            %magnetic co-ords of geographic grid
lonmg(lonmg==pi)=-pi;                             %reset +180 degree values to -180 degree
lonmg1=lonmg;                                     %first copy of lonmg
lonmg2=lonmg;                                     %second copy of lonmg
lonmg2=longch2(lonmg2);                           %reset lonmg2 to between 0 and 2pi
lonmgmin1=min(min(lonmg1));                       %1st minimum magnetic longitude (rads)
lonmgmax1=max(max(lonmg1));                       %1st maximum magnetic longitude (rads)
lonmvli1=sum(lonmv<=lonmgmin1);                   %1st index of lower lonmv value required
lonmvui1=length(lonmv)-sum(lonmv>lonmgmax1)+1;    %1st index of upper lonmv value required
lonmv21=lonmv(lonmvli1:lonmvui1);                 %1st magnetic longitudes at which to calculate densities
lonmgmin2=min(min(lonmg2));                       %2nd minimum magnetic longitude (rads)
lonmgmax2=max(max(lonmg2));                       %2nd maximum magnetic longitude (rads)
lonmvli2=sum(lonmvb<=lonmgmin2);                  %2nd index of lower lonmv value required
lonmvui2=length(lonmvb)-sum(lonmvb>lonmgmax2)+1;  %2nd index of upper lonmv value required
lonmv22=lonmvb(lonmvli2:lonmvui2);                %2nd magnetic longitudes at which to calculate densities
lonmv22=longch(lonmv22);                          %correct 2nd longitude values to between -pi and pi
lons1=length(lonmv21);                            %length of magnetic longitude vector using -pi to pi values
lons2=length(lonmv22);                            %length of magnetic longitude vector using -pi to pi values
%Assume that the correct longitude vector to be used is the smallest of the two calculated
if lons1<=lons2
 lonmv2=lonmv21;
else
 lonmv2=lonmv22;
end
lonmv2(lonmv2==pi)=-pi;
nopx=length(lonmv2);

%% Determine the magnetic latitudes of the boundary region required for which
%% to calculate the required ion densities. 
latmgmin=min(min(latmg));                         %minimum magnetic latitude (degrees)
latmgmax=max(max(latmg));                         %maximum magnetic latitude (degrees)
latmvli=sum(latmv<=latmgmin);                     %index of lower latmv value required
latmvui=length(latmv)-sum(latmv>latmgmax)+1;      %index of upper latmv value required
latmv2=latmv(latmvli:latmvui);                    %magnetic latitudes at which to calculate densities
nopy=length(latmv2);

%% Set up results matrices for the ion densities. Create one O+, H+ and He+ 2D array for
%% each required longitude. (Trying to do the indexing in 3D was giving me a headache!)
nopz=length(disv);
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 latmv2 by disv.
latmatm=latmv2'*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
lshmn=lshm(latmatm>=0);                                     %northern hemisphere L-shells
if ~isempty(lshmn)
 lshminn=min(min(lshmn));                                   %minimum northern L-shell
 lshmaxn=max(max(lshmn));                                   %maximum northern L-shell
else
 lshminn=0;                            
 lshmaxn=12;
end 
lshms=lshm(latmatm<0);                                      %southern hemisphere L-shells
if ~isempty(lshms)
 lshmins=min(min(lshms));                                   %minimum southern L-shell
 lshmaxs=max(max(lshms));                                   %maximum southern L-shell
else
 lshmins=0;                            
 lshmaxs=12;
end 

%% Calculate the magnetic latitudes of the minimum and maximum L-shells
if lshminn~=0
 latminn=lshpar(lshminn,-1,r0);
else
 latminn=0;
end 
latmaxn=lshpar(lshmaxn,-1,r0);
if lshmins~=0
 latmins=-lshpar(lshmins,-1,r0);
else
 latmins=0;
end 
latmaxs=-lshpar(lshmaxs,-1,r0);
latrt=lshpar(lshcp,-1,r0);

%% Determine in what region each point in the region bound by latmv2 by disv. These
%% are the plasmasphere inside L=min(7,lshcp) (northern and southern hemispheres) 
%% and outside L=min(7,lshcp) (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. Also calculate the O+ equatorial density, but as this is only
%% needed below lshcp don't need to replace -1 terms.
timm=timut2ml2(timmat*ones(1,nopx),lonmv2,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 latmv2 by disv.
[xmat,ymat]=pol2cart(latmatm,ones(nopy,1)*disv);

%% If any points within the region bound by latmv2 by disv are inside lshcp, then
%% calculate the required equatorial ion densities and then determine the densities
%% along the required field lines. Use these results to interpolate to the results
%% matrices.
% 'Plasmasphere Results'
if any(any(poi1))
 if any(any(poi11))                                    %consider northern hemisphere points
  li=sum(lshps<=lshminn);                              %lower index 
  ui=length(lshps)-sum(lshps>min(lshmaxn,lshcp))+1;    %upper index
  lshps2=ones(nopx,1)*lshps(li:ui);                    %required lshps values
  hyeqp2=hyeqp(:,li:ui);                               %required hyeqp values
  if hes==1
   heeqp2=heeqp(:,li:ui);                              %required heeqp values
  else
   heeqp2=zeros(size(hyeqp2));
  end
  oxeqp2=oxeqp(:,li:ui);                               %required oxeqp values
  lonm=lonmv2(:)*ones(1,size(hyeqp2,2));               %required magnetic longitude values 
  latm=ones(nopx,1)*latmb(li:ui);                      %required magnetic latitude values
  [oxdenpsn,hydenpsn,hedenpsn,dismatpsn,latmatmpsn]=plasirg12(hyeqp2(:),heeqp2(:),oxeqp2(:),lshps2(:),lonm(:),latm(:),timmat,day,month,year,nop,iricf,gms,hes,sms);
  for i=1:nopx                                         %loop through each magnetic longitude
   ix=i:nopx:(i+nopx*(ui-li));                         %index of i'th rows
   if i==1
    [xmatpsn,ymatpsn]=pol2cart(latmatmpsn(ix,:),dismatpsn(ix,:));                    %northern hemisphere field lines cartesian co-ords
   end 
   xmatpsn(xmatpsn<1e-3)=0;                                                          %correction for round off error
   oxdp=griddata(xmatpsn,ymatpsn,log10(oxdenpsn(ix,:)),xmat(poi11),ymat(poi11));     %northern O+ densities
   hydp=griddata(xmatpsn,ymatpsn,log10(hydenpsn(ix,:)),xmat(poi11),ymat(poi11));     %northern H+ densities
   if hes==1
    hedp=griddata(xmatpsn,ymatpsn,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
  end
 end
  
 if any(any(poi12))                                    %consider southern hemisphere points
  li=sum(lshps<=lshmins);                              %lower index 
  ui=length(lshps)-sum(lshps>min(lshmaxs,lshcp))+1;    %upper index
  lshps2=ones(nopx,1)*lshps(li:ui);                    %required lshps values
  hyeqp2=hyeqp(:,li:ui);                               %required hyeqp values
  if hes==1
   heeqp2=heeqp(:,li:ui);                              %required heeqp values
  end
  oxeqp2=oxeqp(:,li:ui);                               %required oxeqp values
  lonm=lonmv2(:)*ones(1,size(hyeqp2,2));               %required magnetic longitude values 
  latm=ones(nopx,1)*latmb(li:ui);                      %required magnetic latitude values
  [oxdenpss,hydenpss,hedenpss,dismatpss,latmatmpss]=plasirg12(hyeqp2(:),heeqp2(:),oxeqp2(:),lshps2(:),lonm(:),-latm(:),timmat,day,month,year,nop,iricf,gms,hes,sms);
  for i=1:nopx                                         %loop through each magnetic longitude
   ix=i:nopx:(i+nopx*(ui-li));                         %index of i'th rows
   if i==1
    [xmatpss,ymatpss]=pol2cart(-latmatmpss(ix,:),dismatpss(ix,:));                   %southern hemisphere field lines cartesian co-ords
   end 
   xmatpss(xmatpss<1e-3)=0;                                                          %correction for round off error
   oxdp=griddata(xmatpss,ymatpss,log10(oxdenpss(ix,:)),xmat(poi12),ymat(poi12));     %southern O+ densities
   hydp=griddata(xmatpss,ymatpss,log10(hydenpss(ix,:)),xmat(poi12),ymat(poi12));     %southern H+ densities
   if hes==1
    hedp=griddata(xmatpss,ymatpss,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
  end
 end
end 

%'Polar' results
latmpv=deg2rad([0:80 80.5:0.5:90]);                    %polar base magnetic latitude vector (degrees)
if any(any(poi2))
 if any(any(poi21))                                    %consider northern hemisphere points
  li=sum(latmpv<=max(latminn,latrt));                  %lower index 
  ui=length(latmpv)-sum(latmpv>latmaxn)+1;             %upper index
  ui(ui==length(latmpv))=length(latmpv);
  latmpol=latmpv(li:ui);                                              %required polar magnetic latitudes
  lshpol=lshpar(-ones(1,(ui-li+1)),latmpol,r0*ones(1,(ui-li+1)));     %L-shell of polar points

  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)*latmpol;                           %required magnetic latitude values
  lonm=lonmv2(:)*ones(1,length(latmpol));              %required magnetic longitude values 
  [oxdenpon,hydenpon,hedenpon,latmatmpon,dismatpon]=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*(ui-li));                         %index of i'th rows
   if i==1
    [xmatpon,ymatpon]=pol2cart(latmatmpon(ix,:),dismatpon(ix,:));                    %northern hemisphere field lines cartesian co-ords
   end 
   xmatpon(xmatpon<1e-3)=0;                                                          %correction for round off error
   oxdp=griddata(xmatpon,ymatpon,log10(oxdenpon(ix,:)),xmat(poi21),ymat(poi21));     %northern O+ densities
   hydp=griddata(xmatpon,ymatpon,log10(hydenpon(ix,:)),xmat(poi21),ymat(poi21));     %northern H+ densities
   if hes==1
    hedp=griddata(xmatpon,ymatpon,log10(hedenpon(ix,:)),xmat(poi21),ymat(poi21));    %northern He+ densities
   end
   eval(['oxden',int2str(i),'(poi21)=10.^oxdp;']);     %record northern O+ densities for i'th longitude
   eval(['hyden',int2str(i),'(poi21)=10.^hydp;']);     %record northern H+ densities for i'th longitude
   if hes==1
    eval(['heden',int2str(i),'(poi21)=10.^hedp;']);     %record northern He+ densities for i'th longitude
   end
  end
 end
 
 if any(any(poi22))                                    %consider southern hemisphere points
  li=sum(-latmpv>=min(latmins,-latrt));                %lower index 
  ui=length(latmpv)-sum(-latmpv<latmaxs)+1;            %upper index
  ui(ui>length(latmpv))=length(latmpv);
  latmpol=-latmpv(li:ui);                                             %required polar magnetic latitudes
  lshpol=lshpar(-ones(1,(ui-li+1)),latmpol,r0*ones(1,(ui-li+1)));     %L-shell of polar points

  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)*latmpol;                           %required magnetic latitude values
  lonm=lonmv2(:)*ones(1,length(latmpol));              %required magnetic longitude values 
  [oxdenpos,hydenpos,hedenpos,latmatmpos,dismatpos]=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*(ui-li));                         %index of i'th rows
   if i==1
    [xmatpos,ymatpos]=pol2cart(latmatmpos(ix,:),dismatpos(ix,:));                    %southern hemisphere field lines cartesian co-ords
   end 
   xmatpos(xmatpos<1e-3)=0;                                                          %correction for round off error
   oxdp=griddata(xmatpos,ymatpos,log10(oxdenpos(ix,:)),xmat(poi22),ymat(poi22));     %southern O+ densities
   hydp=griddata(xmatpos,ymatpos,log10(hydenpos(ix,:)),xmat(poi22),ymat(poi22));     %southern H+ densities
   if hes==1
    hedp=griddata(xmatpos,ymatpos,log10(hedenpos(ix,:)),xmat(poi22),ymat(poi22));    %southern He+ densities
   end 
   eval(['oxden',int2str(i),'(poi22)=10.^oxdp;']);     %record southern O+ densities for i'th longitude
   eval(['hyden',int2str(i),'(poi22)=10.^hydp;']);     %record southern H+ densities for i'th longitude
   if hes==1
    eval(['heden',int2str(i),'(poi22)=10.^hedp;']);     %record southern He+ densities for i'th longitude
   end
  end
 end 
end

%% Combine the results for each longitude slice into a single 3D matrix
oxdent=zeros(nopy,nopx,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 

%% Use linear interpolation of each disv slice to determine the ion densities
%% at the required points of interest.

%Calculate whether to use longch2 or not. Use if maximum longitude difference is reduced.
lonf1=max(abs(diff(lonmv2)));
lonf2=max(abs(diff(longch2(lonmv2))));
if lonf1<=lonf2
 oxden=pwint2m(lonmv2,latmv2,oxdent,lonmg,latmg);
 hyden=pwint2m(lonmv2,latmv2,hydent,lonmg,latmg);
 if hes==1
  heden=pwint2m(lonmv2,latmv2,hedent,lonmg,latmg);
 end 
else 
 oxden=pwint2m(longch2(lonmv2),latmv2,oxdent,longch2(lonmg),latmg);
 hyden=pwint2m(longch2(lonmv2),latmv2,hydent,longch2(lonmg),latmg);
 if hes==1
  heden=pwint2m(longch2(lonmv2),latmv2,hedent,longch2(lonmg),latmg);
 end
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