%% plsd1.m is designed to be called by the plasdyn#.m models. It calculates the
%% oxygen density along the field lines from the base. This will be done from r0
%% to along the field line to the equator. Starting at r0, estimate the O+ density
%% at each point of interest using the previous value below, starting from r0
%% (where chemical equilibrium has been used to determine O+).  
%% NOTE for points below rb, the inputted value for the oxygen is used. The 
%% diffusion calculation is only used for points above rb.
%%
%% plsd1.m is based on the earlier plsd8.m, but now includes He+.
%%

function [oxden,oxdenc,shoxm,shoxmc]=plsd1(oxden,oxdenc,hyden,hydenc,heden,hedenc,dflm,dismat,latmatm,g0,tmpi,tmpe,tmpt,tmpic,tmpec,tmptc,rb,rbc)

% oxden = oxygen ion density profiles along field lines (initial column has value at base of field line) (O+/m^3)
% oxdenc = conjugate ion oxygen density profiles along field lines (initial column has value at base of field line) (O+/m^3)
% hyden = hydrogen ion density profile (H+/m^3)
% hydenc = conjugate hydrogen ion density profile (H+/m^3)
% heden = helium ion density profile (He+/m^3)
% hedenc = conjugate helium ion density profile (He+/m^3)
% dflm = distance along field line of points on field line from base (km)
% dismat = radial distance of points on field line from base (km)
% latmatm = magnetic latitude of points on field line from base (rads)
% g0 = gravity of points on field line (m/s^2)
% tmpi = ion temperature of points on field line (K)
% tmpe = electron temperature of points on field line (K)
% tmpt = total temperature of points on field line (K)
% tmpic = conjugate ion temperature of points on field line (K)
% tmpec = conjugate electron temperature of points on field line (K)
% tmptc = conjugate total temperature of points on field line (K)
% rb = radial distance below which inputted O+ values to be used (km)
% rbc = radial distance below which inputted conjugate O+ values to be used (km)
%
% oxden = calculated oxygen density values (O+/m^3)
% oxdenc = calculated conjugate oxygen density values (O+/m^3)
% shoxm = calculated oxygen scale height (km)
% shoxmc = calculated conjugate oxygen scale height (km)

%Required constants
re=6371;                                     %Earth's radius
k=1.38e-23;                                  %Boltzman constant (J/K) 
amu=1.66e-27;                                %atomic mass unit (kg)

[p1,dvl]=size(oxden);                        %required size parameters
shoxm=zeros(p1,dvl);                         %oxygen scale heights
shoxmc=shoxm;                                %conjugate oxygen scale heights

for i=1:(dvl-1) 
 ip=i+1;                                     %i plus one
 
 % Determine the location of the top of the i'th section. This is done using
 % distance along the field line from the base of field line.
 ddis=dflm(:,ip)-dflm(:,i);                  %step size increment for i'th section 
 
 %% First consider the Northern Hemisphere Case
 poi1=dismat(:,ip)>rb & oxden(:,i)>1;        %terms were diffusion to be used
 poi2=dismat(:,ip)<=rb;                      %terms below rb where only scale height required 
 poi3=~poi1 & ~poi2;                         %points where O+ already reduced to minimum value of 1 O+/m^3
 
 %diffusion calculation
 if any(poi1)                     
  mb=(16*oxden(poi1,i)+4*heden(poi1,i)+hyden(poi1,i))./(oxden(poi1,i)+heden(poi1,i)+hyden(poi1,i));      %mean ion density
  shoxm(poi1,i)=1e-3*(k*tmpi(poi1,i)./(16*amu*g0(poi1,i)))./(1-mb.*tmpe(poi1,i)./(16*tmpt(poi1,i)));     %oxygen scale height (km)  
  oxdeni=oxden(poi1,i).*(tmpt(poi1,i)./tmpt(poi1,ip)).*exp(-ddis(poi1)./shoxm(poi1,i));                  %oxygen density at top of section 
  oxdeni(oxdeni<1)=1;                                                                                    %reset values below 1 O+/^3 to 1
  oxden(poi1,ip)=oxdeni;                                                                                 %record i'th O+ densities
 end 
   
 %scale height only calculation
 if any(poi2)                                                                                    
  mb=(16*oxden(poi2,i)+4*heden(poi2,i)+hyden(poi2,i))./(oxden(poi2,i)+heden(poi2,i)+hyden(poi2,i));      %mean ion density
  shoxm(poi2,i)=1e-3*(k*tmpi(poi2,i)./(16*amu*g0(poi2,i)))./(1-mb.*tmpe(poi2,i)./(16*tmpt(poi2,i)));     %oxygen scale height (km)  
 end
 
 %no change in values
 if any(poi3)
  oxden(poi3,ip)=1;                          %1 O+/m^3 case
  shoxm(poi3,i)=shoxm(poi3,(i-1));           %record previous scale height for current value
 end 
 
 %% Second consider the Southern Hemisphere case
 poi1c=dismat(:,ip)>rbc & oxdenc(:,i)>1;     %conjugate terms were diffusion to be used
 poi2c=dismat(:,ip)<=rbc;                    %conjugate terms below rb where only scale height required 
 poi3c=~poi1c & ~poi2c;                      %conjugate points where O+ already reduced to minimum value of 1 O+/m^3

 %diffusion calculation
 if any(poi1c)                     
  mb=(16*oxdenc(poi1c,i)+4*hedenc(poi1c,i)+hydenc(poi1c,i))./(oxdenc(poi1c,i)+hedenc(poi1c,i)+hydenc(poi1c,i));    %mean ion density
  shoxmc(poi1c,i)=1e-3*(k*tmpic(poi1c,i)./(16*amu*g0(poi1c,i)))./(1-mb.*tmpec(poi1c,i)./(16*tmptc(poi1c,i)));      %oxygen scale height (km)  
  oxdeni=oxdenc(poi1c,i).*(tmptc(poi1c,i)./tmptc(poi1c,ip)).*exp(-ddis(poi1c)./shoxmc(poi1c,i));                   %oxygen density at top of section 
  oxdeni(oxdeni<1)=1;                                                                                              %reset values below 1 O+/^3 to 1
  oxdenc(poi1c,ip)=oxdeni;                                                                                         %record i'th O+ densities
 end 
   
 %scale height only calculation
 if any(poi2c)                                                                                    
  mb=(16*oxdenc(poi2c,i)+4*hedenc(poi2c,i)+hydenc(poi2c,i))./(oxdenc(poi2c,i)+hedenc(poi2c,i)+hydenc(poi2c,i));    %mean ion density
  shoxmc(poi2c,i)=1e-3*(k*tmpic(poi2c,i)./(16*amu*g0(poi2c,i)))./(1-mb.*tmpec(poi2c,i)./(16*tmptc(poi2c,i)));      %oxygen scale height (km)  
 end
 
 %no change in values
 if any(poi3c)
  oxdenc(poi3c,ip)=1;                        %1 O+/m^3 case
  shoxmc(poi3c,i)=shoxmc(poi3c,(i-1));       %record previous scale height for current value
 end 

end 

%% Now need to combine the O+ densities from each point of interest and its
%% conjugate point, so that the values agree at the field apex (ie magnetic 
%% equator). Use a weighting based on magnetic latitude, but choosen so that
%% weighting does not come into play until latitude reduced to half that of
%% the base latitude.
latmb=latmatm(:,1);                                         %base magnetic latitude
wt1=0.5*(abs(latmatm./(0.5*latmb*ones(1,dvl)))+1);          %weight for local term
wt1(wt1>1)=1;                                               %set maximum value to 1                                      
wt2=1-wt1;                                                  %weight for conjugate term 
oxdeno=oxden;                                               %copy of original oxygen density
oxdenco=oxdenc;                                             %copy of original conjugate oxygen density
oxden=wt1.*oxdeno+wt2.*oxdenco;                             %corrected oxygen density at points of interest
oxdenc=wt1.*oxdenco+wt2.*oxdeno;                            %corrected conjugate oxygen density at points of interest
