%% fltpr (flux tube parameters) determines the following 
%% parameters of the specified flux tube 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.
%% 1) volume (m^3)
%% 2) cross sectional area at re+350km (m^2)
%% 3) cross sectional area at equatorial plane (m^2)
%%

function [vol,a350,aeq]=fltpr(lsh)

% lsh = L-shell of point of interest
%
% vol = volume of flux tube (m^3)
% a350 = cross sectional area at re+350km (m^2)
% aeq = cross sectional area at equatorial plane (m^2)

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)

%Calculate the volume of the flux tube
vt=r0./(re*lsh);
vol=fls*((32*(1e3*re)*lsh.^4)/(35*bo)).*sqrt(1-vt).*(1+(1/2)*vt+(3/8)*vt.^2+(5/16)*vt.^3);

%Calculate the cross sectional area at r0
thb=acos(sqrt((1./lsh)*(r0/re)));            %latitude of base (rads)
bfb=bo*(re/r0)^3*sqrt(1+3*sin(thb).^2);      %magnetic field of base (T)
a350=fls./bfb;                               %cross section area of base (m^2)

%Calculate the cross sectional area at equator
bfe=bo*(1./lsh).^3;                          %magnetic field at equator (T)
aeq=fls./bfe;                                %cross section area at equator (m^2)
