%% ftdiff.m (flux tube difference) calculates the difference in
%% the radial and angular width between to flux tubes on the same
%% potential surface. The difference is a comparison between the
%% expected size of the flux tube in the equatorial plane at the
%% initial time and at the final time and the size that results
%% from rotating the initial times area's corners along the their
%% potentials to the final time. The area of this rotated points
%% are then compared to the expected area at the new location.
%%
%% ftdiff.m is based on the earlier ftdiff4.m.
%% 

function [lshf,mltf,rat]=ftdiff(lsh,mlt,kp,dts,sfs)

% lsh = L-shell of points of interest
% mlt = MLT, magnetic local time (fractional hours)
% kp = Kp (scalar)
% dts = time step in seconds to advance points of interest (must be divideable by 50 seconds)
% sfs = scaling factor switch (0 = original McIllwain, 1 = Liemohn) (scalar)
%
% lshf = final L-shell
% mltf = final MLT
% rat = ratio of final area to original area

re=6371;                                     %Earth radi (km)
re2=1e3*re;                                  %Earth radi (m)
r0=re+350;                                   %Base radial distance

sz=size(lsh);                                %size of input
lsh=lsh(:);                                  %convert lsh to a column vector
mlt=mlt(:);                                  %convert mlt to a column vector

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*(1./lsh).^3;                           %magnetic field of points of interest (T)
csa=fls./bf;                                 %cross sectional area of point of interest (m^2)

%Determine the location of the vertices of the cross sectional area in the equatorial
%plane, taking it to be square in shape.
dms=sqrt(csa);                               %dimensions of square (m)
dth=dms./(re2*lsh);                          %angular width of square (rads)
lsht=[lsh lsh+dms/re2 lsh];                  %total L-shell matrix
mltt=[mlt mlt mlt+(12/pi)*dth];              %total MLT matrix

%Determine the new co-ordinates of the points of interest
[lshn,mltn]=orbpath(lsht,mltt,kp,dts,0,0,sfs);
 
%Calculate the area in the equatorial plane of the of the new points of interest
%after rotation. Compare this the expected area for the standard magnetic unit flux.
%Note that to get accurate results, need to convert to cartesian co-ordinates to
%calculate the area of the rotated points using the cross product rule.
bf2=bo*(1./lshn(:,1)).^3;                    %magnetic field after rotation (T)
csa2=fls./bf2;                               %cross sectional area after rotation (m^2)
[xn,yn]=pol2cart(mltn*pi/12,lshn*re2);       %cartesian co-ordinates of new points
dax=xn(:,2)-xn(:,1);                         %first area vector x-coordinates
day=yn(:,2)-yn(:,1);                         %first area vector y-coordinates
dbx=xn(:,3)-xn(:,1);                         %second area vector x-coordinates
dby=yn(:,3)-yn(:,1);                         %second area vector y-coordinates
csan=abs((dax.*dby)-(day.*dbx));             %area of rotated points
rat=csan./csa2;                              %ratio of new area to area of vertices to expected area      

%Reshape to original input size
lshf=reshape(lshn(:,1),sz);
mltf=reshape(mltn(:,1),sz);
rat=reshape(rat,sz);
