%% plasdyn11.m is designed to be called by plasdyn1.m. Based on iriap9.m,
%% it calculates the required IRI2001 parameters along the given field lines.
%%

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

% dismat = radial distance of points of interests along field lines
% latmatm = magnetaic 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 modelta parametaers
% #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;         %reseta values below 350km altitude to 350km
altmat=dismat-re;                       %altitude matrix (km)

%Break points of interest into regions depending on HmF2 of field line
poi1=altmat<(hmf2*sv);                                 %points below HmF2
poi1c=altmat<(hmf2c*sv);                               %conjugate points below HmF2
poi2=altmat==(hmf2*sv);                                %points equal to HmF2
poi2c=altmat==(hmf2c*sv);                              %conjugate points equal to HmF2
poi3=altmat>(hmf2*sv) & altmat<=((hmf2+100)*sv);       %points betaween HmF2 and HmF2+100km altitude     
poi3c=altmat>(hmf2c*sv) & altmat<=((hmf2c+100)*sv);    %conjugate points betaween HmF2 and HmF2+100km altitude     
poi4=hmf2<=250;                                        %check for HmF2 below 250km
poi4c=hmf2c<=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-eden350(p1y)./nmf2(p1y));                           %square root term
 sqtm(sqtm<1e-3)=1e-3;                                           %seta lower limit
 ym=(hmf2(p1y)-350)./sqtm;                                       %half width
 alt=altmat(poi1);                                               %altitude terms
 if size(alt,1)~=size(ym,1);
  alt=alt(:);
 end 
 eden(poi1)=nmf2(p1y).*(1-((alt-hmf2(p1y))./ym).^2);             %electron density
end 

%Conjugate points below HmF2
if any(p1yc)
 sqtm=sqrt(1-eden350c(p1yc)./nmf2c(p1yc));                       %square root term
 sqtm(sqtm<1e-3)=1e-3;                                           %seta lower limit
 ymc=(hmf2c(p1yc)-350)./sqtm;                                    %half width
 alt=altmat(poi1c);                                              %altitude terms
 if size(alt,1)~=size(ymc,1);
  alt=alt(:);
 end
 edenc(poi1c)=nmf2c(p1yc).*(1-((alt-hmf2c(p1yc))./ymc).^2);   %electron density
end 

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

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

%Points above HmF2 and below HmF2+100km altitude
if any(p3y)
 alt=altmat(poi3);            
 eden(poi3)=modbent(alt(:),hmf2(p3y),nmf2(p3y),beta(p3y),eta(p3y),delta(p3y),zeta(p3y)); 
end 

%Conjugate points above HmF2 and below HmF2+100km altitude
if any(p3yc)
 alt=altmat(poi3c);            
 edenc(poi3c)=modbent(alt(:),hmf2c(p3yc),nmf2c(p3yc),betac(p3yc),etac(p3yc),deltac(p3yc),zetac(p3yc)); 
end 

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

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

%Add in a small correction for small L-shells whose apex altitude is
%less than 1500km, so that electron densities agree at magnetaic 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;                                   %reseta edeno values to zero if edenoc is zero
 edenoc(edeno==0)=0;                                   %reseta 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;                                 %reseta HmF2 to minimu value
 hmf2pc(poia)=hmf2pmin;                                %reseta conjugate HmF2 to minimu value
end 
