%% modmcilpp.m is designed to work out the location of the
%% 'plasmapause' for the McIlwain 1986 electric potential 
%% model. The 'plasmapause' in this sense refers of the last
%% closed potential for the given Kp value.
%%
%% 28/12/2001 modify to allow user to select use of original McIllwain scaling
%%            factor, or new approach based on Liemohn et al.
%%

function [lshb,epb]=modmcilpp(mlt,kp,sfs)

% mlt = MLT values at which to determine plasmapause
% kp = Kp value of interest (scalar)
% sfs = scaling factor switch (0 = original McIllwain, 1 = Liemohn) (scalar)
% lshb = plasmapasue L-shell values for given MLT
% epb = electric potential of last closed field line (V)

re=6371;                 %Earth radi
sz=size(mlt);
mlt=mlt(:);

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

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

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

%% To determine the potential of the last closed loop, assume it oocurs
%% at MLT=18 when the radial potential gradient derivative is zero. Find this
%% by iterative steps starting at L=2. Note that MLT=18 corresponds to 
%% lon = (3/2)pi => sin(lon)= -1 and cos(lon)=0.
lonip=18*(pi/12);                                 %convert 18 MLT from time to longitude
lshp=1.5:0.1:9;
rar=s1+s3*kr;                                     %auroral ring radi (Earth radi)
hi=(v6*rar./lshp).^gm;                            %internal value from shielding function
h=1./(1+hi);                                      %shielding term
g1=lshp.*(-v1)+v3;                                %1st potential term
if sfs==0
 g2=1+v4*kr;                                      %2nd potential term
else
 g2=modmclie(kp); 
end

dg1dr=-v1;                                        %1st potential term r derivative
dg3dr=v5./lshp.^2;                                %3rd potential term r derivative
dhdr=(gm./lshp).*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)

lsh0=interp1(dgdr,lshp,0,'spline');                %L-shell of zero derivative 
epb=modmcle(18,lsh0,kp,sfs);                      %potential of zero derivative

%% To determine the plasmapause at the required MLT values, loop through each
%% and determine its L-shell using interpolation.
lshb=zeros(size(mlt));                            %results matrix of plasmapause location      
nopp=length(lshp);                                %number of points in lshp
for i=1:length(mlt)
 lshpi=lshp;                                      %copy of lshp
 ep=modmcle(mlt(i),lshp,kp,sfs);                  %electric potential at given MLT as a function of L-shell
 epc=sum(diff(ep)<0);                             %check for 'turning' points
 
 if epc~=0                                        %remove terms beyond 'turning' point so interp1 will work
  epo=ep;
  lshpio=lshpi;
  ep((nopp-epc+1):nopp)=[];
  lshpi((nopp-epc+1):nopp)=[];
 end 
 lshb(i)=interp1(ep,lshpi,epb,'spline');           %use interpolation to determine plasmapause at given MLT
 
 ii=1
 keyboard
 
 %If interp1 returns NaN, then assume occurs at maximum ep value
 if isnan(lshb(i))
  %lshvf=lshp(1):0.01:max(lshp);
  lshvf=(lsh0-2):0.01:(lsh0+2);
  epof=interp1(lshp,epo,lshvf,'spline');  
  lshb(i)=lshvf(max(find(max(epof)==epof)));
 end
  
end

%% If any NaN, replace with MLT=18 L-shell value
%lshc=isnan(lshb);
%if any(lshc)
% lshb(lshc)=lsh0;
%end 

lshb=reshape(lshb,sz);                            %resize lshb to original input size
