%% hytgrid2.m calculates the new total H+ and He+ (see 23/11/2001 note) ion 
%% content of the flux tubes on the given grid after the required time 
%% increment.
%% hytgrid2.m is based on the earlier hytgrid6.m.
%%
%% During this time increment, the original tubes are moved along 
%% their potential surfaces and their content is modified to reflect
%% tube volume compression/rarefraction. It then interpolates
%% the  of the new locations back to the original grid.
%%
%% NOTE - before hytgrid2.m can be call, hytgrid1.m needs to have
%%        been called to generate the required lookup tables for the
%%        required L-shell and MLT vectors.
%%
%% 25/1/2002 - remove code that looks for stationary points and replaces them
%%             with the background content before interpolation and then 
%%             re-enters them afterwards (orginally included to stop stationary
%%             points acting as a continious source of plasma).
%%

function hytotn=hytgrid2(hytot,lshv,mltv,kp,lsht,mltt,ratt,kpt,cprf,bkd,iontp)

% hytot = original total ion content matrix (size length(lshv) by length(mltv))
% lshv = L-shell value vector 
% mltv = magnetic local time vector (fractional hours)
% kp = Kp value (scalar)
% lsht = L-shell lookup table array
% mltt = MLT lookup table array
% ratt = volume ratio lookup table array
% kpt = Kp vector for table lookup arrays
% cprf = include compression/rarefaction of plasma ( 1 = yes, other = no)
% bkd = background density ( 1 = bgden.m, other = constant value of one (use only for display of idealized plasmasphere))
% iontp = ion type (0 = H+, 1 = He+) (scalar)
%
% hytotn = new total ion content matrix 

lshv=lshv(:);                      %convert L-shell vector to column vector
mltv=mltv(:);                      %convert MLT vector to column vector
nopy=length(lshv);                 %number of L-shell values
nopx=length(mltv);                 %number of MLT values
svy=ones(nopy,1);
svx=ones(1,nopx);
lshm=lshv*svx;                     %L-shell matrix
mltm=svy*mltv';                    %MLT matrix

xvm=lshm.*cos(mltm*pi/12);         %x-coordinates
yvm=lshm.*sin(mltm*pi/12);         %y-coordinates

lshmx=max(lshv);                   %maximum L-shell

kpt2=kpt+0.5*(kpt(2)-kpt(1));      %lookup vector for Kp values
kpi=1+sum(kpt2<kp);                %index of required Kpt value (uses nearest one to inputted Kp)
lshn=squeeze(lsht(:,:,kpi));       %new L-shell
mltn=squeeze(mltt(:,:,kpi));       %new MLT
ratn=squeeze(ratt(:,:,kpi));       %ratio of final area to original area

lshn(lshn>lshmx)=nan;              %reset points with L-shells greater than the maximum L-shell to NaN
mltn(lshn>lshmx)=nan;              %reset points with L-shells greater than the maximum L-shell to NaN

xvmn=lshn.*cos(mltn*pi/12);        %new x-coordinates
yvmn=lshn.*sin(mltn*pi/12);        %new y-coordinates
  
%Add in a ring of plasma at inner point-0.1 L-shell at each MLT (or L=1 if inner 
%point-0.1 is inside that) and a ring of plasma at L=lshmx+0.1 to allow GRIDDATA to 
%correctly determine the tube content. Take the L=lshmx+0.1 value as being equal to 
%background space density, to allow GRIDDATA to correctly determine the tube 
%content.  
lshin=lshn(1,:)-0.1;                              %inner L-shell
lshin(lshin<1)=1;                                 %reset values less than L=1 to L=1
xvm2=lshin.*cos(mltn(1,:)*pi/12);                 %inner ring x-coordinates
yvm2=lshin.*sin(mltn(1,:)*pi/12);                 %inner ring y-coordinates
hytot2=hytot(1,:);                                %inner ring tube content

%Background tube content at lshmx+0.1
if bkd==1
 [hyeq,heeq,hytbk,hetot]=bgden((lshmx+0.1),kp);
 if iontp==1
  hytbk=hetot;
 end
else
 hytbk=1;
end

xvm3=(lshmx+0.1)*cos(mltm(1,:)*pi/12);            %outer ring x-coordinates
yvm3=(lshmx+0.1)*sin(mltm(1,:)*pi/12);            %outer ring y-coordinates
hytot3=hytbk*svx;                                 %outer ring tube content
 
%Add in a further semi-ring just outside the maximum new L-shells
lshr=lshn(nopy,:);                                %maximum new L-shell for each MLT
mltr=mltn(nopy,:);                                %corresponding MLT
xvm4=(lshr+0.05).*cos(mltr*pi/12);                %semi-ring x-coordinate
yvm4=(lshr+0.05).*sin(mltr*pi/12);                %semi-ring y-coordinate

%Background tube content at semi-ring
if bkd==1
 [hyeq,heeq,hytbk,hetot]=bgden(lshr,kp);
 if iontp==1
  hytbk=hetot;
 end
 hytot4=hytbk;                                    %outer semi-ring tube content
else
 hytot4=svx;
end
  
%Convert all vector to column vectors 
xvm5=[xvmn(:);xvm2(:);xvm3(:);xvm4(:)];
yvm5=[yvmn(:);yvm2(:);yvm3(:);yvm4(:)];

%Determine if any NaN terms exist. Remove any that do.
pnn=isnan(xvm5);
xvm5(pnn)=[];
yvm5(pnn)=[];
 
%Included effects of tube volume compression/rarefraction if required.
if cprf==1
 hytot=hytot./ratn;                 
end

%Use GRIDDATA to interpolate data back to original grid
hytot5=[hytot(:);hytot2(:);hytot3(:);hytot4(:)];
hytot5(pnn)=[];
hytotn=griddata(xvm5,yvm5,hytot5,xvm,yvm,'nearest');  



