%% iriap9.m estimates the electron density at the point of interest and at
%% the conjugate points. The inputted values correspond to the required 
%% parameters at 350km (see iriap8.m for one way to generate values), from
%% which the required parameters are calculated using a parabolic layer 
%% for points below HmF2 and the Bent Model (see IRI2000) for points above 
%% HmF2. For points below 350km altitude, the 350km value is returned. For
%% points 100km above the HmF2 zero is returned.
%%
%% This approach assumes that the field lines are vertical (which is a close
%% approximation at low altitudes) and that latitudinal variation is obtained
%% by mixing the point of interest and its conjugate point, with the weighting
%% being based on magnetic latitude (ie most change is therfore around the
%% apex of the field line)
%%
%% See iriap8.m for description of format types. 

function [eden,edenc,hmf2p,hmf2pc]=iriap9(timi,dismat,latmatm,eden350,hmf2,nmf2,beta,eta,delta,zeta,eden350c,hmf2c,nmf2c,betac,etac,deltac,zetac)

% timi = time index to be considered (scalar)
% dismat = radial distance of points of interests along field lines
% latmatm = magnetic latitude of points of interest along field line
% eden350 = electron density at 350km (m^-3)
% hmf2 = HmF2 (km)
% nmf2 = NmF2 (m^-3)
% beta, eta, delta, zeta = Bent model parameters
% #c = conjugate terms
%
% eden = electron density at point of interest (m^-3)
% edenc = conjugate electron density (m^-3)
% hmf2p = HmF2 of points of interest
% hmf2pc = conjuage HmF2 of points of interest

re=6371;                                %Earth radi (km)

%Results matrices
eden=zeros(size(latmatm));              
edenc=eden;

sv=ones(1,size(latmatm,2));

dismat(dismat<(re+350))=re+350;         %reset values below 350km altitude to 350km
altmat=dismat-re;                       %altitude matrix (km)

%Select out values for required time increment
edn=eden350(:,timi);
ednc=eden350c(:,timi);
h2=hmf2(:,timi);
h2c=hmf2c(:,timi);
n2=nmf2(:,timi);
n2c=nmf2c(:,timi);
bt=beta(:,timi);
btc=betac(:,timi);
et=eta(:,timi);
etc=etac(:,timi);
del=delta(:,timi);
delc=deltac(:,timi);
zt=zeta(:,timi);
ztc=zetac(:,timi);

%Record rerquired HmF2
hmf2p=h2;
hmf2pc=h2c;

%Break points of interest into regions depending on HmF2 of field line
poi1=altmat<(h2*sv);                              %points below HmF2
poi1c=altmat<(h2c*sv);                            %conjugate points below HmF2
poi2=altmat==(h2*sv);                             %points equal to HmF2
poi2c=altmat==(h2c*sv);                           %conjugate points equal to HmF2
poi3=altmat>(h2*sv) & altmat<=((h2+100)*sv);      %points between HmF2 and HmF2+100km altitude     
poi3c=altmat>(h2c*sv) & altmat<=((h2c+100)*sv);   %conjugate points between HmF2 and HmF2+100km altitude     
poi4=h2<=250;                                     %check for HmF2 below 250km
poi4c=h2c<=250;                                   %check for HmF2 below 250km

%Find row and column index of poi's
[p1y,p1x]=find(poi1);              
[p1yc,p1xc]=find(poi1c);              
[p2y,p2x]=find(poi2);              
[p2yc,p2xc]=find(poi2c); 
[p3y,p3x]=find(poi3);              
[p3yc,p3xc]=find(poi3c); 

%Points below HmF2
if any(p1y)
 sqtm=sqrt(1-edn(p1y)./n2(p1y));                            %square root term
 sqtm(sqtm<1e-3)=1e-3;                                      %set lower limit
 ym=(h2(p1y)-350)./sqtm;                                    %half width
 alt=altmat(poi1);                                          %altitude terms
 eden(poi1)=n2(p1y).*(1-((alt(:)-h2(p1y))./ym).^2);         %electron density
end 

%Conjugate points below HmF2
if any(p1yc)
 sqtm=sqrt(1-ednc(p1yc)./n2c(p1yc));                        %square root term
 sqtm(sqtm<1e-3)=1e-3;                                      %set lower limit
 ymc=(h2c(p1yc)-350)./sqtm;                                 %half width
 alt=altmat(poi1c);                                         %altitude terms
 edenc(poi1c)=n2c(p1yc).*(1-((alt(:)-h2c(p1yc))./ymc).^2);  %electron density
end 

%Points equal to HmF2
if any(p2y)
 eden(poi2)=n2(p2y);          
end

%Conjugate points equal to HmF2
if any(p2yc)
 edenc(poi2c)=n2c(p2yc);      
end

%Points above HmF2 and below HmF2+100km altitude
if any(p3y)
 alt=altmat(poi3);            
 eden(poi3)=modbent(alt(:),h2(p3y),n2(p3y),bt(p3y),et(p3y),del(p3y),zt(p3y)); 
end 

%Conjugate points above HmF2 and below HmF2+100km altitude
if any(p3yc)
 alt=altmat(poi3c);            
 edenc(poi3c)=modbent(alt(:),h2c(p3yc),n2c(p3yc),btc(p3yc),etc(p3yc),delc(p3yc),ztc(p3yc)); 
end 

%Case when HmF2 is below 250km
if any(poi4)
 eden(poi4,1)=edn(poi4);
end 

%Case when conjugate HmF2 is below 250km
if any(poi4c)
 edenc(poi4c,1)=ednc(poi4c);
end 

%Add in a small correction for small L-shells whose apex altitude is
%less than 1500km, so that electron densities agree at magnetic equator.
altmx=max(altmat,[],2);                                %apex altitude of field line
poia=altmx<1500;                                       %apex terms less than 1500km
if any(poia)
 wt=0.5*(1+latmatm(poia,:)./(latmatm(poia,1)*sv));     %weighting term
 wtc=(1500-altmx(poia))./800;                          %altitude correction term
 wtc(wtc>1)=1;                                         %upper wtc value
 wt1=wt.*(wtc*sv);                                     %1st weighting term
 wt2=1-wt1;                                            %2nd weighting term
 edeno=eden(poia,:);                                   %initial electron densities
 edenoc=edenc(poia,:);                                 %initial conjugate electron densities
 edeno(edenoc==0)=0;                                   %reset edeno values to zero if edenoc is zero
 edenoc(edeno==0)=0;                                   %reset edenoc values to zero if edeno is zero
 eden(poia,:)=wt1.*edeno+wt2.*edenoc;                  %corrected electron densities
 edenc(poia,:)=wt2.*edeno+wt1.*edenoc;                 %corrected conjugate electron densities
 hmf2pmin=min(hmf2p(poia),hmf2pc(poia));               %minimum HmF2 along field line
 hmf2p(poia)=hmf2pmin;                                 %reset HmF2 to minimu value
 hmf2pc(poia)=hmf2pmin;                                %reset conjugate HmF2 to minimu value
end 
