%% tranht.m calculates the O+/H+ transition height and density from the inputted
%% oxygen and hydrogen ion density profiles. It assumes that each row of input
%% corresponds to separate density profile of interest. Note that outputted 
%% transition height corresponds to the  inputted distance matrix ie if altitude
%% is inputted, then the returned transition height is also altitude etc. The 
%% same is true for the density ie has the same units as input. If transition 
%% height is outside given outputted values zero is returned for both outputs.
%%
%% NOTE - the transition height returned is the lowest one along each profile
%%        ie it ignores any transitions that may occur at higher altitude/radial 
%%        distance.

function [trh,trd]=tranht(dist,oxden,hyden)

% dist = distance matrix
% oxden = oxygen ion density
% hyden = hydrogen ion density
%
% trh = transition height
% trd = transition density

nop=size(dist);                    %size of input
poi=oxden>=hyden;                  %points where O+ density is greater than H+

%The following is a check through poi matrix to sweep through and remove poi=1 
%terms that occur after the original O+/H+ transition point (ie removes higher
%trnasition events if they occur).
spv=ones(nop(1),1);                %'sweep' vector
for i=1:nop(2)                     %loop through each column
 spi=poi(:,i);                     %poi values in i'th column
 poic=spv~=spi;                    %check for terms that are not equal
 if any(poic)
  poi(poic,i)=0;                   %reset poi that are not equal to zero
  spv(poic)=0;                     %reset corresponding spv also to zero
 end 
end 

pois=sum(poi,2);                  %sum of poi
poib=(pois-1)*nop(1)+(1:nop(1))'; %index of required terms below transition height
poia=poib+nop(1);                 %index of required terms above transition height 

%Check for terms where transition height is above the maximum distance value
%for the field line. In these cases the output will be returned as zero.
tc=pois~=nop(2);                   %terms with transition height NOT above max distance
trh=zeros(nop(1),1);               %transition height values
trd=trh;                           %transition density values

%required parameters above and below the transition height
y1=dist(poib(tc));
y2=dist(poia(tc));
x1=hyden(poib(tc));
x2=hyden(poia(tc));
x3=oxden(poib(tc));
x4=oxden(poia(tc));

%required gradient terms
m1=(y2-y1)./(x2-x1);
m2=(y2-y1)./(x4-x3);

%required output
trh(tc)=y1+((m1.*m2)./(m1-m2)).*(x1-x3);
trd(tc)=x1+(m2./(m1-m2)).*(x1-x3);
