%% modmcilpp2.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.
%%
%% Unlike modmcilpp.m, modmcilpp2.m determines the last closed
%% potential by conducting search for the location of the minimum
%% velocity (given by sqrt(vth^2+vrd^2)) using fminsearch, rather
%% than asuming it was located at MLT=18 the location where the
%% radial velocity is zero. Also, it only allows user to the 
%% Liemohn scaling factor (sfs = 1).
%%
%% Author: Phillip A. Webb
%% Last updated: September 16 2003 
%%

function [lshb,epb]=modmcilpp2(mlt,kp)

% mlt = MLT values at which to determine plasmapause
% kp = Kp value of interest (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

%% To determine the potential of the last closed loop, search for the location
%% of the minimum resultant velocity given by sqrt(vth^2+vrd^2).
c0=[18 6];               %initial starting location (MLT and L-shell)
options=optimset('Display','off');
[c,velmin]=fminsearch('e5dvfit',c0,options,kp);        %search for minimum
mlt0=c(1);                                             %MLT of minimum
lsh0=c(2);                                             %L-shell of minimum
epb=modmcle(mlt0,lsh0,kp,1);                           %potential of minimum velocity point

%% To determine the plasmapause at the required MLT values, loop through each
%% and determine its L-shell using interpolation.
lshp=1.5:0.1:9;                                   %required L-shell values
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,1);                    %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
   
end
  
lshb=reshape(lshb,sz);                            %resize lshb to original input size
