%% plasipt13.m calculates the O+ up the specified field lines. 
%% Note that for the input oxeq, if 0 is entered then no equatorial
%% correction is done (for use in polar regions), otherwise enter a 
%% column vector with the equatorial O+ densities to correct O+ 
%% profiles.
%% The pcs is useful for the polar regions, since the H+ and He+ profiles 
%% are fixed and therefore the O+ density profile can be determine exactly
%% on the first iteration. In that case, the returned scale height
%% is 0 for all points on the polar field lines.
%%

function [oxden,shoxm]=plasipt13(oxden,oxeq,shoxm,hyden,heden,lsh,dismat,dflm,latmatm,g0,tmpi,tmpe,tmpt,rb,pcs)

% oxden = oxygen ion density profiles along field lines (initial column has value at base of field line) (O+/m^3)
% oxeq = oxygen ion density at apex of field line (O+/m^3) or 0 if no correction required
% shoxm = previously calculated oxygen scale height (km)
% hyden = hydrogen ion density profile (H+/m^3)
% heden = helium ion density profile (He+/m^3)
% lsh = L-shell of given field line (column vector)
% 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)
% rb = radial distance below which inputted O+ values to be used (km) (column vector)
% pcs = polar calculation switch (0 = use previous calculated values for polar regions, other = calculate O+) (scalar)
%
% oxden = calculated oxygen density values (O+/m^3)
% shoxm = calculated 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

poip=lsh>13 & pcs==0;                        %polar O+ profiles that don't need to be calculated again

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  
 poi1=~poip & dismat(:,ip)>rb & oxden(:,i)>1;     %terms were diffusion to be used
 poi2=~poip & dismat(:,ip)<=rb;                   %terms below rb where only scale height required 
 poi3=~poip & ~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 
 
end

%% Now need to combine the O+ densities from each field line with its apex 
%% value, 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. Note that no change is made if equatorial O+ density
%% is 0 O+/m^3;
if oxeq(1)~=0
 latmb=latmatm(:,1);                              %base magnetic latitude
 wto=abs(latmatm./(0.5*latmb*ones(1,dvl)));       %weight for local term
 wto(wto>1)=1;                                    %set maximum value to 1                                     
 oxden=wto.*oxden+(1-wto).*(oxeq*ones(1,dvl));    %corrected oxygen density at points of interest 
end 

