%% plasisc1.m determines the electron densities for the specified magnetic
%% longitudinal 'slice' at the given UT, day, month and year using data files
%% produced by plasimod(2).m. 
%% The user specifies the number of angular seperations, the minimum and maximum
%% altitude to be used and the number of points to break the radial profile
%% into (using a log scale so smaller divisions closer to Earth).
%%
%% 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 - the y to x ratio required to plot the results with a 'round' Earth
%%        is 7.5/18.5.
%%
%% NOTE - minimum allowed altitude is 100km, to a maximum given by the rounded
%%        down lshcp value from the plasmod(2) results. 
%%
%% 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,latmatm,dismat,hes]=plasisc1(lonlsh,latmnop,altmin,altmax,radnop,timmat,day,month,year,nop)

% lonlsh = magnetic longitude of slice (rads) (scalar)
% latmnop = number of radial paths to use ie how many steps to break magnetic latitude into (scalar)
% altmin = minimum altitude of each radial path (km) (scalar)
% altmax = maximum altitude of each radial path (km) (scalar)
% radnop = number of points to break radial direction to use (scalar)
% timmat = univerisal time (fractional hours) (scalar)
% day = day (scalar)
% month = month (scalar)
% year = year (scalar)
% nop = number of points per field line (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)
% latmatm = magnetic latitude of grid
% dismat = radial distance of grid
% hes = He+ simulated (0 = no, 1 = yes) (scalar)

re=6371;                                %Earth radi

%Check altitude values are within allowed values
dno=daynum(day,month,year);                       %day number
lshcp=plasires2(dno,year);
altmx=1e3*floor((lshcp-1)*re/1e3);
altmin(altmin<100)=100;
altmax(altmax>altmx)=altmx;

%Construct grid of points at which to calculate densities
latmatmv=linspace(-pi/2,pi/2,latmnop);            %magnetic latitude vector
dv=logspace(-1.4,0,radnop);                        
dv=dv-dv(1);       
dv(radnop)=1;
dismatv=re+altmin+(altmax-altmin)*dv;             %radial distance vector
latmatm=latmatmv'*ones(1,radnop);                 %magnetic latitude matrix
dismat=ones(latmnop,1)*dismatv;                   %radial distance matrix

%Result matrices
eden=zeros(latmnop,radnop);               
oxden=eden;
hyden=eden;
heden=eden;

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

%% Calculate the electron densities below 350km using the IRI2001
poi=dismat<=(re+350);                                       %points below 350km
if any(any(poi))
 svi=ones(sum(sum(poi)),1);
 [lonmat,latmat]=mag2geo3(lonlsh*svi,latmatm(poi));         %geographic co-ordinates     
 eden(poi)=modiri5(lonmat,latmat,timmat*svi,day*svi,month*svi,year*svi,dismat(poi),500,iricf,gms);
end

%% Calculate the ion and electron densites above 350km
if any(any(~poi))
 [eden(~poi),oxden(~poi),hyden(~poi),heden(~poi)]=plasisc11(lonlsh,latmatmv,dismatv,timmat,day,month,year,nop,iricf,gms,hes,sms);
end
