%% teccal3.m is a modified version teccal1.m and teccal2.m which 
%% calculates the total TEC along each ray path (assumed to be 
%% specificed by each row in the input), the 'mean ionospheric height'
%% (which doesn't actually have to be in the ionosphere!) and the
%% longitude and latitude of this point. The mean ionospheric height
%% is not calculated if switch mihc==0.
%%

function [tec,mih,lonmih,latmih]=teccal3(eden,dismat,lonmat,latmat,mihc)

% eden = electron density matrix (e-/m^3)
% dismat = radial distance matrix (km)
% lonmat = geographic longitude (rads)
% latmat = geographic latitude (rads)
% mihc = 0=>don't calculate mean ionospheric height, other=> do calculate (scalar)
%
% tec = Total Electron Content (TECU)
% mih = 'mean ionospheric height' vector (km)
% lonmih = 'mean ionospheric height' longitude vector (rads)
% latmih = 'mean ionospheric height' latitude vector (rads)

re=6371;                      %Earth radi (km)
sz=size(eden);                %size of input

%% Since teccal3 does not assume ray paths are radial in direction, need to 
%% calculate the actual distance between points along the ray path.
[xc,yc,zc]=sph2cart(lonmat,latmat,dismat);                  %cartesian co-ords
dsdf=sqrt((xc(:,2:sz(2))-xc(:,1:(sz(2)-1))).^2+(yc(:,2:sz(2))-yc(:,1:(sz(2)-1))).^2+(zc(:,2:sz(2))-zc(:,1:(sz(2)-1))).^2);                                              %distance between points along ray paths

eddf=0.5*(eden(:,1:(sz(2)-1))+eden(:,2:sz(2)));   %mean electron density matrix
tecm=[zeros(sz(1),1) eddf.*dsdf/1e13];            %electron density difference matrix (TECU)
tecmcs=cumsum(tecm,2);                            %cumulative sum of TEC matrix
tec=tecmcs(:,sz(2));                              %total TEC for each ray path

%%Calculate the mean ionospheric height along each field line if required
if mihc==0
 mih=[];
 lonmih=[];
 latmih=[];
else
 poi=sum(tecmcs<(0.5*tec*ones(1,sz(2))),2);       %index of points below 1/2 maximum TEC
 poim=(poi-1)*sz(1)+(1:sz(1))';                   %matrix index of points below 1/2 maximum TEC
 xcb=xc(poim);                                    %x co-ords below 1/2 maximum
 xca=xc(poim+sz(1));                              %x co-ords above 1/2 maximum
 ycb=yc(poim);                                    %y co-ords below 1/2 maximum
 yca=yc(poim+sz(1));                              %y co-ords above 1/2 maximum
 zcb=zc(poim);                                    %z co-ords below 1/2 maximum
 zca=zc(poim+sz(1));                              %z co-ords above 1/2 maximum
 teccsb=tecmcs(poim);                             %cummulative TEC below 1/2 maximum                           
 teccsa=tecmcs(poim+sz(1));                       %cummulative TEC above 1/2 maximum                           
 xgd=(xca-xcb)./(teccsa-teccsb);                  %x co-ords gradient term
 ygd=(yca-ycb)./(teccsa-teccsb);                  %y co-ords gradient term
 zgd=(zca-zcb)./(teccsa-teccsb);                  %z co-ords gradient term
 xmih=xcb+xgd.*(0.5*tec-teccsb);                  %x co-ords of 1/2 TEC maximum
 ymih=ycb+ygd.*(0.5*tec-teccsb);                  %y co-ords of 1/2 TEC maximum
 zmih=zcb+zgd.*(0.5*tec-teccsb);                  %z co-ords of 1/2 TEC maximum
 [lonmih,latmih,mihr]=cart2sph(xmih,ymih,zmih);   %return to polar co-ordinates
 mih=mihr-re;                                     %mean ionospheric height (altitude)
end