%% plasigb1.m determines the global total electron content (TEC). If the user inputs
%% one or more TEC altitude bands, the returned TEC is 3D, which each altitude bands
%% TEC per layer. It also returns the electron and ion densities in a 3D grid with
%% dimensions latitude x longitude x altitude. Also calculated is the mean ionospheric
%% height, O+/H+ transition height and density and the TEC above and below the O+/H+
%% transition height.
%%
%% NOTE - while the electron densities are returned are returned for all
%%        points, the oxygen and hydrogen densities are only returned above
%%        350km (zero is returned below 350km). 
%%
%% NOTE - minimum allowed altitude is 100km, to a maximum given by the rounded
%%        down lshcp value from the plasmod(2) results. 
%%       
%% NOTE - for global TEC outputs, use plasigb1.m.
%%
%% NOTE - plasigb1 uses a standard size longitude and latitude steps, the values inputted
%%        by the user are used to interpolate the densities from the standard grid. 
%%        Consequently computation time basically independent of grid size given by user.
%%
%%
%% 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 [tec,tectot,mih,hytrh,hytrd,hetrh,hetrd,tecoh,techyp,techep,lonmat,latmat,eden,oxden,hyden,heden,disv,hes]=plasigb1(tecalt,lonss,latss,altmin,altmax,radnop,timmat,day,month,year,lshcp,nop)

% tecalt = TEC altitude vector (0 = total TEC) (km)
% lonss = longitude grid step size (degrees)
% latss = latitude grid step size (degrees)
% altmin = minimum altitude of each radial path (km) (scalar)
% altmax = maximum altitude of each radial path (km) (scalar)
% radnop = number of points to use in radial direction  (scalar)
% timmat = univerisal time (fractional hours) (scalar)
% day = day (scalar)
% month = month (scalar)
% year = year (scalar)
% lshcp = maximum L-shell for calculation of profiles (scalar)
% nop = number of points per field line (scalar)
%
% tec = Total Electron Content results matrix (TECU)
% tectot = Total TEC (same as tec if tecalt==0)
% mih = mean ionospheric height (km)
% hytrh = O+/H+ transition height (km)
% hytrd = O+/H+ transition density (m^-3)
% hetrh = O+/He+ transition height (km)
% hetrd = O+/He+ transition density (m^-3)
% tecoh = TEC below and above the O+/H+ transition height (TECU)
% techyp = percentage of TEC due to H+
% techep = percentage of TEC due to He+
% lonmat = geographic longitude matrix of grid (rads)
% latmat = geographic latitude matrix of grid (rads)
% eden = electron density (e-/m^3)
% oxden = O+ density (O+/m^3)
% hyden = H+ density (H+/m^3)
% heden = H+ density (He+/m^3)
% disv = radial distance vector of grid (km)
% hes = He+ simulated (0 = no, 1 = yes) (scalar)

re=6371;                                %Earth radi

%Check altitude values are within allowed values
altmx=1e3*floor((lshcp-1)*re/1e3);
altmin(altmin<100)=100;
altmax(altmax>altmx)=altmx;

%Construct vectors describing grid of points at which to calculate densities
latv=deg2rad(-90:latss:90);                       %geographic latitude vector (rads)
lonv=deg2rad(-180:lonss:(180-lonss));             %geographic longitude vector (rads)
[lonmat,latmat]=meshgrid(lonv,latv);              %longitude and latitude matrices
dv=logspace(-1.4,0,radnop);                        
dv=dv-dv(1);       
dv(radnop)=1;
disv=re+altmin+(altmax-altmin)*dv;                %radial distance vector

%Result matrices
latnop=length(latv);
lonnop=length(lonv);
eden=zeros(latnop,lonnop,radnop);               
oxden=eden;
hyden=eden;
heden=eden;

%Determine parameters used in plasimod.m
dno=daynum(day,month,year);
[gms,iricf,hes,sms]=plasires4(dno,year);

%% Calculate the electron densities below 350km using the IRI2001
poi=disv<=(re+350);
if any(poi)
 dis95v=disv(poi)';
 lon95=ones(length(dis95v),1)*lonmat(:)';
 lat95=ones(length(dis95v),1)*latmat(:)';
 dis95=dis95v*ones(1,latnop*lonnop);
 sv95=ones(size(dis95));
 eden95=modiri5(lon95,lat95,timmat*sv95,day*sv95,month*sv95,year*sv95,dis95,500,iricf,gms);
 for i=1:length(dis95v)
  eden(:,:,i)=reshape(eden95(i,:),latnop,lonnop);
 end 
end

%% Calculate the ion and electron densites above 350km
poi2=disv>(re+350);
if any(poi2)
 [eden35,oxden35,hyden35,heden35]=plasigb11(lonmat,latmat,disv(poi2),timmat,day,month,year,nop,iricf,gms,hes,sms);
 oxden(:,:,poi2)=oxden35;
 hyden(:,:,poi2)=hyden35;
 heden(:,:,poi2)=heden35;
 eden(:,:,poi2)=eden35;
end 

%% Calculate the TEC and Mean Ionospheric height for the region.
eden2=zeros(latnop*lonnop,radnop);                %storage matrix to convert from 3D to 2D
for i=1:lonnop                                    %loop through each longitude slice and place in eden2
 li=1+(i-1)*latnop;
 ui=i*latnop;
 eden2(li:ui,:)=squeeze(eden(:,i,:));
end

[tecm,mih,tectot]=teccal1(eden2,disv,tecalt);     %calculate TEC and mean ionospheric height
mih=reshape(mih,latnop,lonnop);                   %reshape mih to latitude - longitude grid
tectot=reshape(tectot,latnop,lonnop);             %reshape tectot to latitude - longitude grid
if tecalt==0
 tec=tectot;                                      %only Total TEC required
else 
 tecl=length(tecalt)+1;                           %number of TEC altitude bands
 tec=zeros(latnop,lonnop,tecl);                   %3D TEC results matrix
 for i=1:tecl
  tec(:,:,i)=reshape(tecm(:,i),latnop,lonnop);    %reshape i'th tec altitude band to latitude - longitude grid
 end
end 

%% Calculate the O+/H+ and O+/He+ transition heights and densities.
re=6371;                                          %Earth radi (km)
oxden2=zeros(latnop*lonnop,radnop);               %storage matrix to convert from 3D to 2D
hyden2=oxden2;
heden2=oxden2;
for i=1:lonnop                                    %loop through each longitude slice and place in oxden2
 li=1+(i-1)*latnop;
 ui=i*latnop;
 oxden2(li:ui,:)=squeeze(oxden(:,i,:));
 hyden2(li:ui,:)=squeeze(hyden(:,i,:));
 if hes==1
  heden2(li:ui,:)=squeeze(heden(:,i,:));
 end
end
[hytrh,hytrd]=tranht(ones(latnop*lonnop,1)*disv,oxden2,hyden2);      %calculate O+/H+ transition height and density
hytrh=reshape(hytrh,latnop,lonnop)-re;                               %reshape hytrh to latitude - longitude grid and convert to altitude
hytrd=reshape(hytrd,latnop,lonnop)-re;                               %reshape hytrd to latitude - longitude grid

if hes==1
 [hetrh,hetrd]=tranht(ones(latnop*lonnop,1)*disv,oxden2,hyden2);     %calculate O+/He+ transition height and density
 hetrh=reshape(hetrh,latnop,lonnop)-re;                              %reshape hytrh to latitude - longitude grid and convert to altitude
 hetrd=reshape(hetrd,latnop,lonnop)-re;                              %reshape hytrd to latitude - longitude grid
else
 hetrh=nan*ones(size(hytrh));
 hetrd=nan*ones(size(hytrh));
end

%% Calculate the TEC above and below the O+/H+ transition height
tecoh2=teccal2(eden2,disv,hytrh(:));                   %calculate the TEC above and below O+/H+ transition height
tecoh=zeros(latnop,lonnop,2);                          %results array
tecoh(:,:,1)=reshape(tecoh2(:,1),latnop,lonnop);       %reshape lower tecoh values to latitude - longitude grid
tecoh(:,:,2)=reshape(tecoh2(:,2),latnop,lonnop);       %reshape upper tecoh values to latitude - longitude grid

%% Calculate the percentage of the TEC due to H+ and He+
techy=teccal1(hyden2,disv,0);
techyp=100*reshape(techy,latnop,lonnop)./tectot;
if hes==1
 teche=teccal1(heden2,disv,0);
 techep=100*reshape(teche,latnop,lonnop)./tectot;
else
 teche=nan*ones(size(techy));
 techep=nan*ones(size(techyp));
end