%% fltpr2 (flux tube parameters) determines the total tube
%% content from the base at 350km altitude to the equator 
%% for the inputted density profile. 
%% Calculated using the standard magnetic flux, which is 
%% the value at L=3 at r=re+350 km with a cross sectional 
%% value of 1m^2.
%%
 
function [itot,vols,vol]=fltpr2(lsh,dismat,latmatm,den)

% lsh = L-shell of point of interest (scalar)
% dismat = radial distance of points on field line (vector) (km)
% latmatm = magnetic latitude of points on field line (vector) (rads)
% den = density at points on field line (m^-3) 
%
% itot = total tube content
% volt = total tube volume (m^-3)
% vols = volum of each section (m^-3)

re=6371;                                     %Earth radi (m)
r0=re+350;                                   %reference height (km)
bo=30.4e-6;                                  %equatorial magnetic field (T)
lshs=3;                                      %standard L-shell

ths=acos(sqrt((1/lshs)*(r0/re)));            %latitude of standard point (rads)
bfs=bo*(re/r0)^3*sqrt(1+3*sin(ths)^2);       %magnetic field of standard point (T)
fls=bfs;                                     %standard magnetic flux (T.m^2)

disc=dismat<r0;                                        %terms to be removed
den(disc)=[];                                          %reset densities below r0 to 0
dismat(disc)=[];
latmatm(disc)=[];
sz=length(den);                                        %size of input
sv=ones(1,sz);

dflm=distfl2(latmatm,lsh*sv,r0*sv);                    %distance along field line of points

bf=bo*(re./dismat).^3.*sqrt(1+3*sin(latmatm).^2);      %magnetic field of points of interest (T)
csa=fls./bf;                                           %cross sectional area of point of interest (m^2)
ddf=1e3*diff(dflm);                                    %distance between points of interest along field line (m)
vol=ddf.*(csa(2:sz)+csa(1:(sz-1)))/2;                  %volume of each section (m^3)

dens=vol.*(den(2:sz)+den(1:(sz-1)))/2;                 %content of each section
itot=sum(dens);                                        %total content of each field line
vols=sum(vol);