%% orbpath.m calculates the new location of the points of interest
%% in the magnetic equatorial plane after the specific time 
%% intereval. 
%% orbpath.m is based on the earlier orbpath5.m.
%% orbpath.m uses the electric field potential given by 
%% McIlwain Adv. Space Res. Vol 6 1986. It uses either a magnetic 
%% dipole field, or the M2 magnetic field also in the McIlwain paper. 
%% The field is choosen using the bsf switch. The V5 constant has been 
%% changed from 92 to 89.75 so low L-shell point of interest rotates 
%% through 2pi in 24hrs to within 0.02%
%% 
%% The radial and angular velocities are then given by 
%% v = {B X del(potential)}/B^2 as before, except the equations
%% are more complex due to the derivatives of the shielding function.
%% 
%% NOTE: 1) for all points of interest orbpath uses a single Kp and 
%%          dt value.
%%       2) all points below L=1.05 are returned as L=1.05.
%%       3) all points above L=15 are returned as L=15
%%       4) dts must be divideable by 50 (exactly)
%% 
%% NOTE - add in a switch tmd, which if set to 1, results in orbpath.m 
%%        calculating the flow backwards in time. Useful for determining
%%        where a point has come from, rather than where it is going 
%%        (normal case).
%%

function [lshn,mltn,vth,vrd]=orbpath(lsh,mlt,kp,dts,bfs,tmd,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)
% bfs = magnetic field (1 = M2 field, other = dipole)
% tmd = time direction (1 = run time backwards, other = normal forward flow)
% sfs = scaling factor switch (0 = original McIllwain, 1 = Liemohn) (scalar)
%
% lshn = new L-shell of points of interest
% mltn = new MLT
% vth = final angular velocity (km/s)
% vrd = final radial velocity (km/s)

re=6371;                 %Earth radi
dt=50;                   %time increment of each iteration (s)
nop=dts/dt;              %number of iterations required

if tmd==1
 dt=-dt;             
end
 
%Loop through each required iteration 
for i=1:nop 
 %Determine the angular and radial velocities (km/s)
 if sfs==0
  [vth,vrd,vthr]=modmcilvel(lsh,mlt,kp,0,bfs,0);
 else
  [vth,vrd,vthr]=modmcilvel2(lsh,mlt,kp,bfs);
 end

 %Determine the new L-shell and magnetic longitude of point of interest 
 lsh=lsh+dt*vrd/re;                          %new L-shell value
 lsh(lsh<1.05)=1.05;                         %correction for points less than L=1.05
 lsh(lsh>15)=15;                             %correction for points greater than L=15
 mlt=mlt+dt*(12/pi)*vth./(lsh*re);           %new MLT
 mlt=timch(mlt);                             %correct MLT if required
end 

lshn=lsh;                                    %new L-shell
mltn=mlt;                                    %new MLT
if sfs==0
 [vth,vrd]=modmcilvel(lsh,mlt,kp,0,bfs,0);   %final angular and radial velocities (km/s)
else
 [vth,vrd]=modmcilvel2(lsh,mlt,kp,bfs);
end