%% tbcck.m (tube content check) is designed to be called by the plasdyn#.m models. 
%% It calculates the total ions of each speices number along the tube from a base height
%% of 350km and compares it to the known tube content of that species. It then corrects 
%% the speices 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.
%% tbcck.m is based on the earlier tbcck2.m, but now includes He+ (again!)  
%%

function [hydenf,hydenfc,hyeqf,hedenf,hedenfc,heeqf,csa]=tbcck(hyden,hydenc,heden,hedenc,latmatm,dismat,dflm,hytot,hyeq,hetot,heeq,hes)

% hyden = hydrogen ion density at points of interest (m^-3)
% hydenc = conjugate hydrogen ion density at points of interest (m^-3)
% heden = helium ion density at points of interest (m^-3)
% hedenc = conjugate helium 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)
% hetot = total He+ of field line (He+)
% heeq = equatorial He+ density of field line (He+/m^3)
% hes = He+ simulated (0 = no, 1 = yes) (scalar)
%
% hydenf = corrected H+ density profile
% hydenfc = corrected conjugate H+ density profile
% hyeqf = corrected equatorial H+ density of field line (H+/m^3)
% hedenf = corrected He+ density profile
% hedenfc = corrected conjugate He+ density profile
% heeqf = corrected equatorial He+ density of field line (He+/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

if hes==1
 hedens=vol.*(heden(:,2:sz(2))+heden(:,1:(sz(2)-1)))/2;     %He+ content of each section (He+)
 hedensc=vol.*(hedenc(:,2:sz(2))+hedenc(:,1:(sz(2)-1)))/2;  %conjugate He+ content of each section (He+)
 hetotn=sum(hedens,2)+sum(hedensc,2);                       %total He+ content of each field line
end

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)

if hes==1
 cf=hetot./hetotn;                                      %correction factor
 hedenf=heden;                                          %corrected He+ density matrix
 hedenf=(cf*ones(1,sz(2))).*hedenf;                     %corrected He+ density (m^-3)
 hedenfc=hedenc;                                        %corrected conjuage He+ density matrix
 hedenfc=(cf*ones(1,sz(2))).*hedenfc;                   %corrected conjuage He+ density (m^-3)
 heeqf=cf.*heeq;                                        %corrected equatorial He+ density of field line (He+/m^3)
else
 hedenf=heden;
 hedenfc=hedenc;
 heeqf=heeq;
end