%% modmcle.m is the equatorial electric potential model from McIlwain 1986.
%%
%% 28/12/2001 modify to allow user to select use of original McIllwain scaling
%%            factor, or new approach based on Liemohn et al.
%%

function [ep,rar]=modmcle(mlt,lsh,kp,sfs)

% mlt = magnetic local time (fractional hours)
% lsh = L-shell
% kp = Kp value 
% sfs = scaling factor switch (0 = original McIllwain, 1 = Liemohn) (scalar)
% ep = equatorial electric potential (V)
% rar = Auroral ring radi (Earth radi)

%% First need to determine the location of the auroral ring
%constant values
s1=9.8;
s2=-1.4;
s3=-0.9;
s4=-0.3;

mlt2=(pi/12)*mlt;             %MLT in rads

%auroral ring radi
kr=kp./(1+0.1*kp);            %Kr parameter
rar=s1+s2*cos(mlt2)+(s3+s4*cos(mlt2)).*kr;

%% Determine the geoelectric field potential
%constant values
v1=0.8;
v2=0.2;
v3=3;
v4=0.3;
v5=92;
v6=0.8;

h=1./(1+(v6*rar./lsh).^8);              %shielding function

ep1=v1*sin(mlt2)+v2*cos(mlt2);
if sfs==0
 ep2=1+v4*kr;
else
 ep2=modmclie(kp);
end
ep=1e3*(((lsh.*ep1+v3).*ep2.*h)-v5./lsh);


