%% teccal1.m calculates the TEC between the required altitudes
%% bands and mean ionospheric height for the given data sets.
%% For example, if 800 is entered for the altitude band, then
%% the TEC values returned are for 0 to 800km and 800km to the
%% maximum radial distance. If zero is entered as the only band 
%% altitude, then the total TEC is returned. teccal1.m assumes
%% that each row in the input matrices is for a seperate ray 
%% path. It does however assume that the radial distance of 
%% each column on the rays paths is the same, given by the 
%% radial distance vector. The output has the required TECs in
%% the corresponding rows of the eden matrix.

function [tec,mih,tectot]=teccal1(eden,dismv,altv)

% eden = electron density matrix (e-/m^3)
% dismv = radial distance vector (km)
% altv = altitude bands vector (km)
%
% tec = TEC (TECU's)
% mih = mean ionospheric height (km)
% tectot = Total TEC for each ray path
%

re=6371;                      %Earth radi (km)
sz=size(eden);                %size of input
altv=altv(:);                 %make sure altv a column vector
dismv=dismv(:);               %make sure dismv a column vector

if length(altv)==1 & altv==0
 zoc=0;                       %zero check
else
 zoc=1;
end 

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

%%Calculte required TECs
if zoc==0
 tec=tecmcs(:,sz(2));                                  %record final cumsum value it only total TEC required
else
 tec=zeros(sz(1),length(altv)+1);                      %TEC results matrix
 bsz=100;                                              %block size
 nol=ceil(sz(1)/bsz);                                  %number of loops required
 for i=1:nol
  ip1=1+bsz*(i-1);                                     %lower index of i'th block
  if i~=nol
   ip2=bsz*i;                                          %upper index of i'th block
  else 
   ip2=sz(1);
  end 
  tecv=interp1(dismv,tecmcs(ip1:ip2,:)',(altv+re));    %cumulative TEC at required altitudes
  tec(ip1:ip2,:)=diff([zeros(size(tecv,2),1) tecv' tecmcs(ip1:ip2,sz(2))],1,2);    %calculate required TEC's in alttude bands
 end
end 

%%Calculate the mean ionospheric height along each field line
tecmx=tecmcs(:,sz(2));                            %total TEC along ray path
poi=sum(tecmcs<(0.5*tecmx*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
disb=dismv(poi);                                  %dismv value below 1/2 maximum 
disa=dismv(poi+1);                                %dismv value above 1/2 maximum 
teccsb=tecmcs(poim);                              %cummulative TEC below 1/2 maximum                           
teccsa=tecmcs(poim+sz(1));                        %cummulative TEC above 1/2 maximum                           
gd=(disa-disb)./(teccsa-teccsb);                  %gradient term
mih=disb+gd.*(0.5*tecmx-teccsb)-re;               %mean ionospheric height
