%% ftcont.m determines the total ion content of the inputted half flux tube (ie
%% content from base altitude to magnetic equator). Each input row is assumed to 
%% correspond to a seperate flux tube.

function [iontot,iondens]=ftcont(ionden,dismat,dflm,latmatm)

% ionden = ion densities at points on field line (m^-3)
% dismat = radial distance of points on field line (km)
% dflm = distance along feild line from base of points on field line (km)
% latmatm = magnetic latitude of points on field line (rads)
%
% iontot = hemispheric tube content
% iondens = ion cotent of each section

re=6371;                                               %Earth's radius
r0=re+350;                                             %Base radial distance

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)
fls=bo*(re/r0)^3*sqrt(1+3*sin(ths)^2);                 %standard magnetic flux (T.m^2)

sz=size(ionden,2);                                     %column size of input
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,1,2);                                %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)

iondens=vol.*(ionden(:,2:sz)+ionden(:,1:(sz-1)))/2;    %ion content of each section
iontot=sum(iondens,2);                                 %total ion content of each field line
