%% plsd41.m is designed to be called by plsd4.m. It determines the distance
%% along the field line from r0 to the point at which the required ratio of 
%% H+ chemical loss/diffusion loss occurs.
%% 

function z=plsd41(dismat,dflm,tmpi,tmpn,oxden,shox,oxnt,shhy,r0,bet)

% dismat = distance from center of Earth of points of interest (km).
% dflm = distance along field lines of points of interest from r0 (km).
% tmpi = ion temperature at points of interest (K)
% tmpn = neutral temperature at points of interest (K)
% oxden = oxygen ion density at points of interest (O+/m^3)
% shox = oxygen ion scale heights at points of interest (km)
% oxnt = neutral oxygen density at points of interest (O/m^3)
% shhy = hydrogen ion scale heights at points of interest (km)
% r0 = base radial distance of field lines (km)
% bet = required ratio of chemical loss/diffusion loss
%
% z = distance along field line from r0 to H+ diffusion height boundary (km)

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

%Take shhy absolute value
shhy=abs(shhy);

%% Estimate z using at iterative process, which starts at r0 and estimates z,
%% then uses the new z to determine the next estimate of z and so on, It has
%% be found that 5 iterations is adequate to determine an estimation of z.
nop=size(dflm,1);                                               %number of field lines
dfi=size(dflm,2);                                               %number of points of interest per field line
sv1=ones(1,dfi);
r0c=r0-re;                                                      %altitude value of r0

for i=1:5
 if i==1                                                        
  z=zeros(nop,1);                                               %initial z value        
  shoxi=shox(:,1);                                              %required ion oxygen scale heights at r0
  g0=9.8*(re./r0).^2;                                           %required gravity values (in m/s^2)
  shoxni=1e-3*k*tmpn(:,1)./(m1*g0);                             %required neutral oxygen scale heights at r0
  oxdeni=oxden(:,1);                                            %oxygen density at r0
  oxnti=oxnt(:,1);                                              %neutral oxygen density at r0
  tmpii=tmpi(:,1);                                              %ion temperature at r0
  shhyi=shhy(:,1);                                              %required ion hydrogen scale heights at r0
 else 
  %Estimate values the ion and electron temperature and the oxygen ion density
  %using linear interpolation of of each of these values. 
  poi=dflm<=(z*sv1);                                            %finds points in dsm matrix below or equal to z value
  pois1=sum(poi,2);                                             %number of points in each row below or equal to z value 
  pois2=pois1+1;                                                %pois1 plus 1
  pois2(pois2==(dfi+1))=dfi;                                    %reset required for points at 950km
  po1=(pois1-1)*nop+(1:nop)';                                   %index of pois1 points
  po2=(pois2-1)*nop+(1:nop)';                                   %index of pois2 points  
  dfv=dflm(po2)-dflm(po1);                                      %difference vector between points of interest
  dfzf=z-dflm(po1);                                             %difference vector between z value and lower point of interest
    
  %Estimate ion temperature
  m=(tmpi(po2)-tmpi(po1))./dfv;
  tmpii=m.*dfzf+tmpi(po1);
    
  %Estimate neutral temperature
  m=(tmpn(po2)-tmpn(po1))./dfv;
  tmpni=m.*dfzf+tmpn(po1);
    
  %Estimate the oxygen ion density
  m=(log(oxden(po2))-log(oxden(po1)))./dfv;
  oxdeni=exp(m.*dfzf).*oxden(po1);  
  
  %Estimate oxygen ion scale heights 
  m=(shox(po2)-shox(po1))./dfv;
  shoxi=m.*dfzf+shox(po1);
  
  %Estimate hydrogen ion scale heights 
  m=(shhy(po2)-shhy(po1))./dfv;
  shhyi=m.*dfzf+shhy(po1);
    
  %Estimate the required neutral oxygen density
  m=(log(oxnt(po2))-log(oxnt(po1)))./dfv;
  oxnti=exp(m.*dfzf).*oxnt(po1);
  
  %Calculate the required neutral oxygen scale height
  m=(dismat(po2)-dismat(po1))./dfv;
  disi=m.*dfzf+dismat(po1);                           
  g0=9.8*(re./disi).^2;                                         %required gravity values (in m/s^2)
  shoxni=1e-3*k*tmpni./(m1*g0);                                 %neutral oxygen scale heights (km)
 end  
 
 if i~=5
  %Determine new estimate of z
  cet=shoxi.*shoxni./(shoxi+shoxni);
  z=z-cet.*(log(7.5e20.*tmpii.^2./(shhyi.^2.*oxdeni.*oxnti))-log(bet));   
  
  %Correction values for z
  z(z<0)=0;
  poc=z>=dflm(:,dfi);
  if any(poc)
   z(poc)=dflm(poc,dfi)-1;
  end  
 end 
end 
