%% modmcilvel2.m is a modified version of modmcilvel.m that uses the
%% Liemohn JGR 2001 scaling factor in the E5D model. It also does not 
%% use the plasmapasue formation velocity.
%% modmcilvel2.m calculates the angular and radial velocities (km/s),
%% the angular velocity (rads/s) in the magnetic equatorial plane 
%% using the McIlwain electric field model. 
%% It 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
%% with the bfs switch value. The V5 constant 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) all points above L=15 are returned as 15. Their 
%%          corresponding magnetic local times should be ignored. 
%%       2) all points below L=1.05 are returned as L=1.05 
%%

function [vth,vrd,vthr]=modmcilvel2(lsh,mlt,kp,bfs)

% lsh = L-shell of points of interest (km)
% mlt = MLT, magnetic local time (fractional hours)
% kp = Kp (scalar)
% bfs = magnetic field (1 = M2 field, other = dipole)
%
% vth = angular velocity (km/s)
% vrd = radial velocity (km/s)
% vthr = angular velocity (rad/s)

loni=mlt*(pi/12);                            %convert mlt from time to longitude

bo=30.4e-6;                                  %equatorial magnetic field (T)
re=6371;                                     %Earth radi

%Auroral ring constant values
s1=9.8;
s2=-1.4;
s3=-0.9;
s4=-0.3;

kr=kp./(1+0.1*kp);            %Kr parameter

%Determine the magnetic field (T)
if bfs==1
 bf=modmclb(mlt,lsh);
else 
 bf=bo./lsh.^3;                              
end 

%Geoelectric field potential constant values
v1=0.8;
v2=0.2;
v3=3;
v4=0.3;
v5=89.75;
v6=0.8;
gm=8;

rar=s1+s2*cos(loni)+(s3+s4*cos(loni)).*kr;             %auroral ring radi (Earth radi)
hi=(v6*rar./lsh).^gm;                                  %internal value from shielding function
h=1./(1+hi);                                           %shielding term
g1=lsh.*(v1*sin(loni)+v2*cos(loni))+v3;                %1st potential term
g2=modmclie(kp);                                       %Liemohn potential term

dg1dr=v1*sin(loni)+v2*cos(loni);                       %1st potential term r derivative
dg3dr=v5./lsh.^2;                                      %3rd potential term r derivative
dhdr=(gm./lsh).*hi.*h.^2;                              %shielding term r derivative
dgdr=(dg1dr.*g2.*h)+(g1.*g2.*dhdr)+dg3dr;              %radial derivative of the electric field potential function (keV/re)

dg1dt=lsh.*(v1*cos(loni)-v2*sin(loni));                %1st potential term theta derivative
dhdt=gm*sin(loni).*h.^2.*hi.*(s2+kr*s4)./rar;          %shielding term theta derivative
dgdt=(dg1dt.*g2.*h)+(g1.*g2.*dhdt);                    %theta derivative of the electric field potential function (keV/re)

%Determine the angular and radial velocities
vth=(1./(1e3*re*bf)).*dgdr;                            %angular velocity (km/s)
vrd=-(1./(1e3*re*lsh.*bf)).*dgdt;                      %radial velocity (km/s)
vthr=vth./(re*lsh);                                    %angular velocity (rad/s)
