%% tbcck2.m (tube content check) is designed to be called by the 
%% plasdyn#.m models. It calculates the total H+ number along
%% the tube from a base height of 350km and compares it to the
%% known tube content. It then corrects the H+ density along
%% the tube (% correction) so that the two agree. tbcck.m 
%% assumes that each row corresponds to separate field line.
%% It also outputs the cross sectional area of each point
%% along the field line.
%% tbcck2.m is based on tbcck.m, but does not calculate the 
%% correction for He+.

function [hydenf,hydenfc,hyeqf,csa]=tbcck2(hyden,hydenc,latmatm,dismat,dflm,hytot,hyeq)

% hyden = hydrogen ion density at points of interest (m^-3)
% hydenc = conjugate hydrogen ion density at points of interest (m^-3)
% latmatm = magnetic latitude of profile points (rads)
% dismat = radial displacement of profile points (km)
% dflm = distance along field line from base (km)
% hytot = total H+ of field line (H+)
% hyeq = equatorial H+ density of field line (H+/m^3)
%
% hydenf = corrected H+ density profile
% hydenfc = corrected conjugate H+ density profile
% hyeqf = corrected equatorial H+ density of field line (H+/m^3)
% csa = cross sectional area of each point on the field line (m^2)

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(hyden);                                             %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(2))+csa(:,1:(sz(2)-1)))/2;             %volume of each section (m^3)

hydens=vol.*(hyden(:,2:sz(2))+hyden(:,1:(sz(2)-1)))/2;          %H+ content of each section (H+)
hydensc=vol.*(hydenc(:,2:sz(2))+hydenc(:,1:(sz(2)-1)))/2;       %conjugate H+ content of each section (H+)
hytotn=sum(hydens,2)+sum(hydensc,2);                            %total H+ content of each field line

cf=hytot./hytotn;                                       %correction factor
hydenf=hyden;                                           %corrected H+ density matrix
hydenf=(cf*ones(1,sz(2))).*hydenf;                      %corrected H+ density (m^-3)
hydenfc=hydenc;                                         %corrected conjuage H+ density matrix
hydenfc=(cf*ones(1,sz(2))).*hydenfc;                    %corrected conjuage H+ density (m^-3)

hyeqf=cf.*hyeq;                                         %corrected equatorial H+ density of field line (H+/m^3)
