%% modmcilvel.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. 
%% If ac==1 it calculates the required modification to the radial 
%% velcity between 23 MLT and 6 MLT using the "plasmapause formation"
%% radial velocity terms for points outside the specific L-shell for 
%% the given Kp.
%% If the switch ac is set to ac==2 it calculates the angular 
%% acceleration (rads/s^2). Set ac==3 to calculate both of these 
%% paramters.
%% 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) for all points of interest orbpath uses a single Kp and 
%%          dt value.
%%       2) all points above L=15 are returned as 15. Their 
%%          corresponding magnetic local times should be ignored. 
%%       3) all points below L=1.05 are returned as L=1.05 
%%
%% NOTE: For information on the parameters used to determine the 
%%       plasmapause formation velocity, see pwpp1.m and pwpp2.m.
%%

function [vth,vrd,vthr,vrdpp,ath]=modmcilvel(lsh,mlt,kp,ac,bfs,lshr)

% lsh = L-shell of points of interest (km)
% mlt = MLT, magnetic local time (fractional hours)
% kp = Kp (scalar)
% ac = extra calculations switch (1 = plasmapause formation velocity, 2 = angular acceleration, 3 = both, other = neither) (scalar)
% bfs = magnetic field (1 = M2 field, other = dipole)
% lshr = round minimum L-shell plasmapasue formation location to nearest 0.1 (1 = yes, other = no)
%
% vth = angular velocity (km/s)
% vrd = radial velocity (km/s)
% vthr = angular velocity (rad/s)
% vrdpp = radial plasmapause formation velocity (km/s)
% ath = angular acceleration (rad/s^2)

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;
alp=1+kr*v4;

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=1+v4*kr;                                            %2nd 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)

%Calculate the angular acceleration if required
if ac==2 | ac==3
 d2g1drdt=v1*cos(loni)-v2*sin(loni);                                                 %1st potential term r and t derivative
 d2hdrdt=h.^2.*hi.*(gm^2./(lsh.*rar)).*sin(loni).*(s2+kr*s4).*(2*h.*hi-1);           %shielding term r and t derivative
 d2gdrdt=(d2g1drdt.*g2.*h)+(dg1dr.*g2.*dhdt)+(dg1dt.*g2.*dhdr)+(g1.*g2.*d2hdrdt);    %potential term r and t derivative
 ath2=(1./(1e6*re^2*lsh.*bf)).*d2gdrdt;                %angular acceleration (km/s^2)
 ath=ath2./(re*lsh);                                   %angular acceleration (rad/s)
else
 ath=zeros(size(vthr));
end

%Calculate the radial plasmapause formation velocity
vrdpp=zeros(size(vthr));
if ac==1 | ac==3
 cf1=[-2.81e-4 0.0061 0.0199];                         %1st L max angular velocity coefficent
 cf2=[4.06e-4 -0.0077 0.0728];                         %2nd L max angular velocity coefficent
 cf3=[0.0304 -0.7001 5.7852];                          %3rd L max angular velocity coefficent
 cf4=[1.17e-4 -0.0020 3.38e-4];                        %1st L/hr coefficent coefficents
 cf5=[-0.0012 0.0228 0.0345];                          %2nd L/hr coefficent coefficents
 cf6=[0.0011 -0.0192 0.0733];                          %3rd L/hr coefficent coefficents
 
 lshmin=10*ones(size(mlt));                            %minimum L-shell to apply plasmapause formation velocity
 poi1=mlt>=23 | mlt<6;                                 %points of interest between 23 MLT and 6 MLT
 if any(any(poi1))
  clm1=polyval(cf1,kp);                                %1st L coefficent
  clm2=polyval(cf2,kp);                                %2nd L coefficent
  clm3=polyval(cf3,kp);                                %3rd L coefficent
  lshmin(poi1)=polyval([clm1 clm2 clm3],mlt(poi1));    %minimum L-shell
 end 
 
 if lshr==1
  lshmin=round(10*lshmin)/10;                          %round lshmin to nearest 0.1 if required
 end

 poi2=poi1 & lsh>=lshmin;                              %points with non-zero plasmapause formation velocity
 if any(any(poi2))
  cpp1=polyval(cf4,kp);                                          %1st L/hr coefficent
  cpp2=polyval(cf5,kp);                                          %2nd L/hr coefficent
  cpp3=polyval(cf6,kp);                                          %3rd L/hr coefficent
  vrdpp(poi2)=(re/3600)*polyval([cpp1 cpp2 cpp3],mlt(poi2));     %radial plasmapause formation velocity (km/s)  
 end
end 
