%% parktc.m is a conversion program that converts between the tube
%% volume for a given field line used by GPID and that used by 
%% Park. The later is defined as a tube having a cross section of
%% 1 cm^2 at 1000km altitude along the given field line to the 
%% magnetic equator. 

function hytf=parktc(hyti,lsh,sp)

% hyti = tube content of interest
% lsh = L-shell of tube
% sp = conversion switch (0 = convert from Park to GPID, otherwise GPID to Park) (scalar)
% hytf = required tube content

re=6371;                           %Earth radi (km)
r0=re+350;                         %base hieght used by GPID (km)

%Calculate the cross sectional area at 1000km along the field line using the GPID approach.
%GPID uses a standard flux based on 1 m^2 at r0 along the L=3 field line.
sv=ones(size(lsh));
disf=re+1000*sv;                                  %radial distance of 1000km point
latmf=lshpar(lsh,-sv,disf);                       %magnetic latitude of 1000km point
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)
bf=bo*(re./disf).^3.*sqrt(1+3*sin(latmf).^2);     %magnetic field at flux point (T)
arat=fls./bf;                                     %area ratio of 1000km point

%Convert tube content depending on the sp setting
if sp==0
 hytf=2e4*hyti.*arat; 
else
 hytf=0.5e-4*hyti./arat;
end


