%% teccal2.m is a modified version teccal1.m which calculates
%% the TEC above and below the specified altitude (one altitude
%% value per ray path). Designed to be used when the altitude
%% value is varying, if constant for each raypath then teccal1
%% should be used. Note that as with teccal1 that teccal2 assumes
%% each row corresponds to a seperate ray path and that the radial
%% distance of each point along each ray path is the same. The 
%% output has two columns, the first value is the TEC below the
%% required altitude, the second is the TEC above. It also returns
%% the total TEC in a seperate vector.

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

% eden = electron density matrix (e-/m^3)
% dismv = radial distance vector (km)
% altv = altitude vector (km)
%
% tec = TEC (TECU's)
% tectot = total TEC

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 row vector

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

altmp=altv*ones(1,length(dismv));                      %altitude vector of point of interest
altm=ones(length(altv),1)*(dismv-re);                  %altitude comparison matrix
poi=altm<=altmp;                                       %points below point of interest on each ray path
px=sum(poi,2);                                         %x index of poi
px2=px+1;                                              %x index of poi plus one
px2(px2==sz(2))=sz(2);                                 %reset values above maximum to maximum
li=sz(1)*(px-1)+(1:sz(1))';                            %lower index terms
ui=sz(1)*(px2-1)+(1:sz(1))';                           %upper index terms
tec1=tecmcs(li);                                       %lower TEC values
tec2=tecmcs(ui);                                       %upper TEC values
m=(tec2-tec1)./dsdf(px)';                              %TEC gradient term
tecp=tec1+m.*(altv-(dismv(px)'-re));                   %TEC at point of interest on each ray path
tec=[tecp tectot-tecp];                                %required TEC values
