%% orbpath3.m calculates the new location of the points of interest
%% in the magnetic equatorial plane after the specific time 
%% intereval. It uses the approach set out by Wolf in Chapter 10,
%% section 10.5.7 of "Introduction to Space Physics" to determine 
%% the drift potential, which then allows the radial and 
%% longitudinal velocites to be calculated. The dawn to dusk 
%% electric field is determined using the  paper by Gallagher et al.
%% "On the azimuthal variation of core plasma in the equatorial 
%%  magnetosphere".
%% orbpath3.m is based on orbpath.m, but does not include the rotation
%% term for the magnetospheric flow. Thus the bulge is always at 18 MLT.
%% 
%% 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) dts must be divideable by 50 (exactly)

function [lshn,mltn,vth,vrd]=orbpath3(lsh,mlt,kp,dts)

% 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)
%
% lshn = new L-shell of points of interest
% mltn = new MLT
% vth = final angular velocity (km/s)
% vrd = final radial velocity (km/s)

loni=mlt*(pi/12);                            %convert mlt from time to longitude
re=6371;                                     %Earth radi
bo=30.4e-6;                                  %equatorial magnetic field (T)
wo=(2*pi)/(24*3600);                         %angular velocity of Earth (rad/s)
eo=bo*wo*(1e3*re)/(2*(5.6-0.46*kp))^2;       %dawn to dusk electric field (V/m)
dt=50;                                       %time increment of each iteration (s) 
nop=dts/dt;                                  %number of iterations required

%Loop through each required iteration 
for i=1:nop
 bf=bo./lsh.^3;                              %magnetic field at i'th point (T)
 vth=wo*(1e3*re*lsh)+(eo./bf).*sin(loni);    %longitudinal velocity (m/s)
 vrd=-(eo./bf).*cos(loni);                   %radial velocity (m/s)
 lsh=lsh+dt*(vrd/1e3)/re;                    %new L-shell value
 lsh(lsh<1.05)=1.05;                         %correction for points less than L=1.05
 loni=loni+dt*(vth/1e3)./(lsh*re);           %new longitude
 loni=longch(loni);                          %correct longitude values if required
end 

lshn=lsh;                                    %new L-shell
mltn=(12/pi)*loni;                           %new MLT
mltn=timch(mltn);                            %check MLT times are between 0 and 24

bf=bo./lsh.^3;                                    %magnetic field at i'th point (T)
vth=(wo*(1e3*re*lsh)+(eo./bf).*sin(loni))/1e3;    %final longitudinal velocity (km/s)
vrd=-(eo./bf).*cos(loni)/1e3;                     %final radial velocity (km/s)
