%
%  Simple time-dependent model including variable outflow temperature
%   K. Emanuel, October-November 2010
%   Last modified April 2019
%---------------------------------------------------------------------
%
%  Constants and control parameters
%  (Number in parentheses are default values)
%
endtime=12;     % End time (days) of integration (12)
outerwind='y';  % Do you want outer windfield added? ('y' or 'n')
                % Note:  When set to 'y' a red dot will appear in some
                % graphs inidicating transition to outer profile
%
Cd=1.0e-3;      % Drag coefficient (1.0e-3)
Ck=1.0e-3;      % Enthalpy exchange coefficient (1.0e-3)
vcap=600;       % Capping wind speed (m/s) in surface enthalpy flux (600, means no cap)
wsmin=0;        % Gustiness (m/s) in surface flux formulation (0)
shear=0;        % Environmental wind shear (kts; parameterized effect; default=0)
Tt=200;         % Minimum outflow temperature (tropopause temperature) (200.0)
Tb=300;         % Surface temperature (300)
the0=370;       % Sea surface theta_e (370)
themax=345.15;  % Peak theta_e in starting vortex (345.15; Should be only slightly greater than theenv)
theenv=345;     % Ambient boundary layer theta_e (345)
thedecay=100;   % Radial decay rate of theta_e in starting vortex (100)
themid=320;     % Ambient mid-level entropy
rcool=0.5;      % Constant radiative cooling rate (deg/day) (0.5)
wcool=2;        % Clear sky radiative subsidence rate (mm/s) (2)
tauout=1.5;     % Time scale (days) to relax outflow temperature (1.5)
weight=1;       % PBL entropy advection weighting; 0.5 = centered, 1.0 = upstream (1)
eyediff=0.7;    % Sets temperature of eye: 1 = near sold-body rotation; 0 = no adjustment (0.7) 
%
h=5000.0;       % Depth (m) over which surfaces fluxes converge (5000)
r0=600.0e3;     % Outer radius (m) of starting vortex (600e3)
rt=15e3;        % Radius (m) needed in critical outflow formulation (15e3; Note:  Modified later in code)
Ric=1;          % Critical Richardson Number (1)
Diss=0;         % Coefficient of dissipative heating (1=full heating; 0=no heating)(0)
isoexpand=0;    % Coefficient governing pressure dependence of surface saturation entropy (0-1)(0)
%
dt=100.0;       % Time step (s) (100)
afilt=0.1;      % Asselin filter value (0.1)
nr=400;         % Number of radial nodes (M-coordinate model) (400)
%
cp=1005.0;      % Heat capacity at constant pressure (1005)
f=5.0e-5;       % Coriolis parameter (5e-5)
%
%-------------------------------------------------------------------------
%
fwait = waitbar(0,'Program running'); % Displays progress bar
%
rt=rt*min(Cd/Ck,1);  %  Empirical modification of rt to help with smooth time evolution
hi=1./h;
cfac=Ck/Cd;
if cfac == 2         % Just to avoid (removable) singularity
    cfac=1.99;
end    
M0=0.5*f*r0^2;       % Outermost value of angular momentum
dM=M0/nr;
dMi=1./dM;
dti=1./dt;
M=0:dM:M0;           % Specify M grid
nr=max(size(M));
endtime=endtime*24.0*3600.0;
tauouti=1./(24.*3600.*tauout);
rcool=rcool.*cp./(Tb.*24.*3600);
dampnorm=5000.*hi;
shear=shear*0.02/500;
nt=1+floor(endtime./dt);
%
% Base value of surface saturation mixing ratio (added February 2018)
%
Tcs=Tb-273.15;
es0=6.112*exp(17.67*Tcs./(243.5+Tcs));
qs0=0.622*es0./(1000-es0);
%
%  Set initial vortex saturation entropy
%
s0=cp.*log(the0);
se=cp.*(log(theenv));
smid=cp.*log(themid);
s1=cp.*log(themax)-cp.*log(themax./theenv).*(1-exp(-thedecay.*M.^2./M0.^2));
s2=s1;
%
% Set other initial conditions
%
t=0;
dssdM=zeros(1,nr);
dsdM=zeros(1,nr);
%
dsdM1=zeros(1,nr);
dsdM2=zeros(1,nr);
sf=zeros(1,nr);
%
vt=zeros(1,nt);
pt=zeros(1,nt);
rmt=zeros(1,nt);
jmax=10;
r=zeros(1,nr)+50e3;
V=zeros(1,nr)+1;
adv=zeros(1,nr);
u1=zeros(1,nr);
u2=zeros(1,nr);
dMdr1=zeros(1,nr);
dMdr2=dMdr1;
dr3=zeros(1,nr);
jmin=nr;
Vmax=0;
clear w psi
psi=zeros(1,nr);
n=0;
T=zeros(1,nr)+Tt;
T1=T;
T2=T;
x=zeros(1,nr);
clear vtemp
rifac=sqrt(Cd*Ric/Ck);
int=1;
%
% Begin time marching
%
while t <= endtime
    n=n+1;
    t=t+dt;
    %
    if t > int
        int=int+3*24*3600;
        waitbar(int/endtime,fwait,'Program running');
    end    
    %
    smax=s1(jmax);  % Value of entropy at radius of maximum winds
    %
    dsdM1(1:nr-1)=(s1(2:nr)-s1(1:nr-1)).*dMi; % M gradient of PBL entropy
    dsdM2(2:nr)=dsdM1(1:nr-1); % M gradient of PBL entropy
    dsdM1(nr)=dsdM1(nr-1);
    dsdM2(nr)=dsdM2(nr-1);
    dssdM(2:nr-1)=0.5*(s1(3:nr)-s1(1:nr-2)).*dMi; % M gradient of Sat. entropy
    dssdM(1)=0;
    dssdM(nr)=dssdM(nr-1);
    [dum2, jss]=min((Tb-T1).*dssdM);    % Max value of -ds*/dM * (Tb-T)
    dssmin=dssdM(jss);
    dssdM(1:jss-1)=eyediff*dssmin+(1-eyediff).*(dssdM(1:jss-1));  % Set -ds*/dM to max value inside location of max value
    dssdM=min(dssdM,0);     % Make sure -ds*/dM is positive or zero
    %
    % Find outflow temperature by integrating outflow layer critical
    % equation for dT/dM
    %
    if n == 1
        T=Tt+(Tb-Tt).*(smax-s1)./(smax-se); %Used only in first time step
        T1=T;
        T2=T;
    else
        T(1:jmax)=Tt; % Outflow T at and inside radius of maximum wind
        x(jmax)=log(Tb-Tt); % Log of (T_b-T) at radius of maximum winds
        %
        % March dx/dM outward from jmax
        %
        for j=jmax+1:nr
            vrmean=0.5.*(V(j)./r(j)+V(j-1)./max(r(j-1),1));
            x(j)=x(j-1)-(Ric./rt.^2.)*(M(j)-M(j-1))./(vrmean+0.5.*f);
        end  
        T(jmax+1:nr)=Tb-exp(x(jmax+1:nr)); % Calculate T from x
    end
    %
    T=min(T,Tb);    % Do not allow outflow T to exceed Tb
    T=max(T,Tt);    % Do not allow outflow T to be less than Tt
    Tmax=T(jmax);   % Outflow temperature at radius of max winds
    %
    r=sqrt(M./(0.5.*f-(Tb-T2).*dssdM)); % Find radius
    for i=nr-1:-1:1                    % Condition radii to be monotonic
        r(i)=min(r(i),r(i+1));
    end    
    %
    V=M./r-0.5.*f.*r;                  % Find Azimuthal velocity
    V(nr)=2.*V(nr-1)-V(nr-2);
    %Vtr(n,:)=V;
    V=max(V,0);
    [Vmax, jmax]=max(V);               % Find max value of V and location of max value
    [dum, jmin]=min(V(2:nr));  
    %
    vt(n)=Vmax;                 % Time series of Vmax for plotting
    rmt(n)=0.001.*r(jmax);      % Time series of radius of maximum winds
    rtnew=r(jmax)*rifac;        % Revised estimate of rt
    rt=0.7*rt+0.3*rtnew;
    %
    % Calculate forcing of boundary layer entropy and outflow temperature
    %
    Tf=(T-T1).*tauouti;
    %
    % The following makes a smooth transition to vcap, to avoid abrupt changes
    %
    ny=5;
    vy=(V.*vcap)./(V.^ny+vcap.^ny).^(1./ny);
    vy(1)=vy(2);
    %
    ws=sqrt(vy.^2+wsmin.^2);
    ws1=sqrt(V.^2+wsmin.^2);
    %
    %  Calculate surface pressure, after first calculating saturation
    %  entropy
    %
    ss=s1;
    ss(1:jmax)=ss(jmax)+(Vmax.^2./(Tb-Tmax)).*(1.0-M(1:jmax)./...
       M(jmax));
    ss=max(ss,s1);
    %
    sx=zeros(1,nr);
    for i=nr-1:-1:1
        sx(i)=sx(i+1)-(Tb-0.5*(T(i)+T(i+1)))*(ss(i+1)-ss(i));
    end    
    alogp=-(0.5.*V.*V+0.5*f*(M-M0)+sx)./(287.*Tb); % Inserted 2/18
    pp0=exp(alogp);
    pt(n)=pp0(1);
    %
    if isoexpand ~= 0
        %
        % Account for variable To along inflow trajectory (added 2/18)
        %
        s0mod=s0+isoexpand*qs0.*2.5e6.*(1./pp0-1)./Tb-287*isoexpand.*alogp; % Changed 0.02 to qs0 2/2018
    else
        s0mod=s0;
    end        
    %
    % calculate forcing (sf) of boundary layer entropy equation
    %
    term1=Cd.*r.*V.*ws1;
    %
    u2(1:nr-1)=0.5*(term1(2:nr)+term1(1:nr-1));
    u1(2:nr)=u2(1:nr-1);
    adv=weight*u2.*dsdM1+(1-weight)*u1.*dsdM2;
    %
    sf=(adv+Ck.*ws.*(s0mod-s1)+Diss.*Cd.*ws1.*V.^2./Tb).*hi;
    sf=sf-shear.*ws.*(s1-smid).*hi;  % Parameterized shear effect
    %
    %  Add thermal damping 
    %
    sf=sf-rcool;
    %
    % Time stepping (leap frog with Asselin filter)
    %
    s3=s1+2.*dt.*sf;
    s3=max(s3,se);
    s1=s2+afilt.*(s1+s3-2.*s2);
    s1=max(s1,se);
    s2=s3;
    T3=T1+2.*dt.*Tf;
    T3=max(T3,Tt);
    T3=min(T3,Tb);
    T1=T2+afilt.*(T1+T3-2.*T2);
    T2=T3;
    %
    %  Some boundary conditions
    %
    s2(nr)=2.*s2(nr-1)-s2(nr-2);
    s1(nr)=2.*s1(nr-1)-s1(nr-2);
    s1(1)=2.*s1(2)-s1(3);
    s2(1)=2.*s2(2)-s2(3);
    %
end
%
% Calculate streamfunction from PBL angular momentum balance
%   and then calculate vertical velocity
%
psi(2:nr-1)=-0.5*Cd.*(r(2:nr-1).*V(2:nr-1)).^2.*(r(3:nr)-r(1:nr-2)).*dMi;
w(2:nr-1)=-2.*(psi(3:nr)-psi(1:nr-2))./(r(3:nr).^2-r(1:nr-2).^2);
w(1)=w(2);
w(nr)=0;
w(3:nr-2)=smooth134(w);
%
% Calculate saturation entropy for plotting
%
ss=s1;
ss(1:jmax)=ss(jmax)+(Vmax.^2./(Tb-Tmax)).*(1.0-M(1:jmax)./...
    M(jmax));
ss=max(ss,s1);
%
nt=max(size(vt));
time=dt:dt:nt*dt;
time=time./(24.*3600);
%
% Calcualte potential intensity from theory
s0p=s0;
if isoexpand ~= 0
    [~,imax]=max(V);
    s0p=s0mod(imax);
end    
vtheory=sqrt((Diss.*Tb+(1-Diss).*Tt).*Ck.*(Tb-Tt).*(s0p-se)./(Tt.*(Cd))); % Raw potential intensity (diagnostic only)
[vtheory] = vpotcalc( vtheory,Ck,Cd,vcap );  % Calculates correction to potential intensity
vtheory=vtheory/sqrt(2);
%
r=0.001*r;  % Reset radius to km
%
%  Add outer windfield if desired
%
if strcmp(outerwind,'y')
    rinner=r;
    Vinner=V;
    psiinner=psi;
    winner=w;
    [~,ivmax]=min(V(3:end));
    %
    [r,V,ro,int] = windmeld(r(1:ivmax),V(1:ivmax),wcool,Cd,f);
    qnew=max(size(r));
    w(int+1:qnew)=-0.001*wcool;
    psi(int+1:qnew)=-0.5*0.001*wcool*1e6*(ro^2-r(int+1:qnew).^2);
    psi(qnew+1:end)=[];
    w(qnew+1:end)=[];
end    
close(fwait)
%
% Plot various quantitites
%
indx=1;
while isempty(indx) == 0
    items={'Maximum wind speed vs. time','Radius of maximum wind vs. time', ...
        'Radial profile of V at final time', ...
        'Radial profile of u at final time','Radial profile of w at final time', ...
        'Radial profiles of entropy at final time'};
    [indx,tf] = listdlg('PromptString','Select a plot:',...
                               'SelectionMode','single',...
                               'CancelString','Stop', ...
                               'ListSize',[220 120], ...
                               'ListString',items);
    if indx == 1
        vtemp(1:nt)=vtheory;
        h2=plot(time,vt,time,vtemp);
        set(gca,'fontweight','bold')
        set(h2,'linewidth',2);
        xlabel('Time (days)')
        ylabel('Azimuthal wind (m/s)')
        title('Azimuthal Wind Time Series')
        %
        vt2=vtheory.*tanh(0.5*max(Cd,Ck)*vtheory*hi*time*24*3600);
        hold on
        j=plot(time,vt2,'--r','linewidth',2);
        hold off
        legend('Azimuthal Wind','Potential Intensity','Analytic approximation','location','best')
    elseif indx == 2
        h1=plot(time,rmt);
        set(gca,'fontweight','bold')
        set(h1,'linewidth',2);
        xlabel('Time (days)')
        ylabel('Radius of maximum wind (km)')
        title('Radius of Maximum Wind Time Series')
    elseif indx == 3
        h1=plot(r,V);
        if strcmp(outerwind,'y')
            hold on
            plot(r(int),V(int),'.r','markersize',18)
            hold off
        end    
        set(gca,'fontweight','bold') 
        set(h1,'linewidth',2);
        xlabel('Radius (km)')
        ylabel('Azimuthal wind (m/s)')
        title('Radial Profile of Azimuthal Wind at Final Time')
    elseif indx == 4
        u=0.001*hi.*psi./r;
        h3=plot(r,u);
        if strcmp(outerwind,'y')
            hold on
            plot(r(int),u(int),'.r','markersize',18)
            hold off
        end    
        set(gca,'fontweight','bold') 
        set(h3,'linewidth',2);
        xlabel('Radius (km)')
        ylabel('Radial wind (m/s)')
        title('Radial Profile of Radial Wind at Final Time')
    elseif indx == 5
        h4=plot(r,w);
        if strcmp(outerwind,'y')
            hold on
            plot(r(int),w(int),'.r','markersize',18)
            hold off
        end            
        set(gca,'fontweight','bold') 
        set(h4,'linewidth',2);
        xlabel('Radius (km)')
        ylabel('Vertical wind (m/s)')
        title('Radial Profile of Vertical Wind at Final Time')
    elseif indx == 6
        %
        h5=plot(rinner,exp(s1./cp),'b',rinner,exp(ss./cp),'r');
        set(gca,'fontweight','bold') 
        set(h5,'linewidth',2);
        xlabel('Radius (km)')
        ylabel('\theta_e (K)')
        legend('PBL \theta_e','\theta_e^*')
        title('Radial Profile of PBL \theta_e and \theta_e^* at Final Time')
    end
end   
%
% Make and display table of various values at final time
%
Vtheory=round(vtheory);
rmv=round(rmt(end));Vmax=round(max(V));pcmin=round(1000*pt(end));
[Umin,iumin]=min(0.001*hi.*psi./r);
Umin=0.1*round(10*Umin);rmu=round(r(iumin));
w(isinf(w))=NaN;
[Wmax,iwmax]=nanmax(w);
Wmax=0.01*round(100*Wmax);rmw=round(r(iwmax));
T=table(Vmax,Vtheory,rmv,Umin,rmu,Wmax,rmw,pcmin) %#ok<NOPTS>
%
%--------------------------------------------------------------------------
function [ vg ] = vpotcalc( vpot,Ck,Cd,vcap )
% 
% This script calculates actual potential intensity (vg) from raw potential
% intensity (vpot), the exchange coefficients, and any capping wind speed
%
vg=0.5*(vcap+vpot);
r=1;
n=0;
while r > 0.0001   % Need to find the solution iteratively
    n=n+1;
    ckn=Ck*min(vcap/vg,1);
    ckcd=ckn/Cd;
    if ckcd == 2
        ckcd=1.99;
    end    
    vgnew=sqrt(2)*vpot*(0.5*ckcd)^(0.5*ckcd/(2-ckcd));
    r=abs(vgnew-vg);
    vg=vgnew;
end  
%
end
%
function [xs] = smooth134(x)
%
%  This applies a smoother over 5 points, weightin 1-3-4-3-1. It returns an array 4 points shorter.
%
m=max(size(x));
xs=zeros(1,m-4);
for i=3:m-2
   xs(i-2)=(x(i-2)+3.*x(i-1)+4.*x(i)+3.*x(i+1)+x(i+2))./12;
end
%
end
%--------------------------------------------------------------------------
function [rnew,Vnew,ro,int] = windmeld(r,V,wcool,Cd,f)
%
% This function takes an inner wind profile given by V as a function
% of r and melds it with an outer windfield according to the Emanuel
% (2004; Lilly volume; see also Chavas et al, J. Atmos. Sci., 2015).
%
% Input:
%    r (radius in km)
%    V (azimuthal wind in m/s)
%    wcool (radiative subsidence velocity in mm/s)
%    Cd (drag coefficient)
%    f (Coriolis parameter, in s^{-1})
%
% Output:
%    rnew (melded radius in km)
%    Vnew (melded azimuthal velocity in m/s)
%    ro   (Outer radius in km)
%    int  (radial index of transition from inner to outer windfield)
%
% Note that this function requires two subroutines: vouternew (written by
% K. Emanuel) and intersections (written bi Douglas Scwartz and available
% at the MATLAB file exchange)
%
%  K. Emanuel, April, 2019
%--------------------------------------------------------------------------
racc=2;  % Accuracy of determination of outer radius (km)
%--------------------------------------------------------------------------
Vmax=max(V);
[~, iend]=min(V(3:end));
ro=6*r(iend); % Initial guess for outer radius (km)
ri=[]; 
n=0;
while isempty(ri)
    n=n+1; % Keep track of number of iterations
    ro=ro-racc;
    %
    [ Vout,rout ] = vouternew(Vmax,f,ro,wcool,Cd );
    [ri,~,iout,jout] = intersections(r,V,rout,Vout,1); % Copyright Douglas Schwartz (see MATLAB File Exchange)
end    
%
itemp=max(size(iout));
jtemp=max(size(jout));
iinner=floor(iout(itemp)); % Inner intersection index of second intersection
iouter=ceil(jout(jtemp));  % Outer intersection index of second intersection
%
% Now meld two curves
%
rnew=[r(1:iinner),rout(iouter:end)];
Vnew=[V(1:iinner),Vout(iouter:end)];
int=iinner;
%
end
%--------------------------------------------------------------------------
function [ V,r ] = vouternew(vm,f,ro,wc,CD )
% Numerical integration of the outer wind profile from simple ODE
% Inputs: vm (maximum wind in m/s); f (Coriolis in s^-1), ro (outer radius
% in km), wc (radiative subsidence rate in mm/s), and drag coefficient.
%
% Outputs:  V (in m/s) versus radius r (im km)
%
% Created January, 2017
%----------------------------------------------------------------------
ass=0.1;    % Asselin filter coefficient
dr=0.002;   % Nondimensional radius increment used for outer wind
%-----------------------------------------------------------------------
ro=ro*1000;
wc=wc*0.001;
v2=max(vm,f*ro);
chi=2*CD*v2/wc;
rdim=vm/f;
rond=ro/rdim;
rnd=0:dr:rond;
q=max(size(rnd));
%
% Integrate outer wind ODE
%
m=zeros(1,q);
v=zeros(1,q);
m(q-1)=0.25.*(rond^2-rnd(q-1).^2);
v(q-1)=m(q-1)/(rond-dr);
for i=q-2:-1:1
    m(i)=m(i+2)-2*dr*(chi*m(i+1)^2/(rond^2-rnd(i+1)^2)-rnd(i+1));
    m(i+1)=m(i+1)+ass.*(m(i)+m(i+2)-2*m(i+1));
    v(i)=m(i)/rnd(i);
    if v(i) > 1
        imax=i;
        break
    end
end
V=vm*v(imax:end);
r=0.001*vm*rnd(imax:end)/f;
%
end
%--------------------------------------------------------------------------
function [x0,y0,iout,jout] = intersections(x1,y1,x2,y2,robust)
%INTERSECTIONS Intersections of curves.
%   Computes the (x,y) locations where two curves intersect.  The curves
%   can be broken with NaNs or have vertical segments.
%
% Example:
%   [X0,Y0] = intersections(X1,Y1,X2,Y2,ROBUST);
%
% where X1 and Y1 are equal-length vectors of at least two points and
% represent curve 1.  Similarly, X2 and Y2 represent curve 2.
% X0 and Y0 are column vectors containing the points at which the two
% curves intersect.
%
% ROBUST (optional) set to 1 or true means to use a slight variation of the
% algorithm that might return duplicates of some intersection points, and
% then remove those duplicates.  The default is true, but since the
% algorithm is slightly slower you can set it to false if you know that
% your curves don't intersect at any segment boundaries.  Also, the robust
% version properly handles parallel and overlapping segments.
%
% The algorithm can return two additional vectors that indicate which
% segment pairs contain intersections and where they are:
%
%   [X0,Y0,I,J] = intersections(X1,Y1,X2,Y2,ROBUST);
%
% For each element of the vector I, I(k) = (segment number of (X1,Y1)) +
% (how far along this segment the intersection is).  For example, if I(k) =
% 45.25 then the intersection lies a quarter of the way between the line
% segment connecting (X1(45),Y1(45)) and (X1(46),Y1(46)).  Similarly for
% the vector J and the segments in (X2,Y2).
%
% You can also get intersections of a curve with itself.  Simply pass in
% only one curve, i.e.,
%
%   [X0,Y0] = intersections(X1,Y1,ROBUST);
%
% where, as before, ROBUST is optional.

% Version: 2.0, 25 May 2017
% Author:  Douglas M. Schwarz
% Email:   dmschwarz=ieee*org, dmschwarz=urgrad*rochester*edu
% Real_email = regexprep(Email,{'=','*'},{'@','.'})


% Theory of operation:
%
% Given two line segments, L1 and L2,
%
%   L1 endpoints:  (x1(1),y1(1)) and (x1(2),y1(2))
%   L2 endpoints:  (x2(1),y2(1)) and (x2(2),y2(2))
%
% we can write four equations with four unknowns and then solve them.  The
% four unknowns are t1, t2, x0 and y0, where (x0,y0) is the intersection of
% L1 and L2, t1 is the distance from the starting point of L1 to the
% intersection relative to the length of L1 and t2 is the distance from the
% starting point of L2 to the intersection relative to the length of L2.
%
% So, the four equations are
%
%    (x1(2) - x1(1))*t1 = x0 - x1(1)
%    (x2(2) - x2(1))*t2 = x0 - x2(1)
%    (y1(2) - y1(1))*t1 = y0 - y1(1)
%    (y2(2) - y2(1))*t2 = y0 - y2(1)
%
% Rearranging and writing in matrix form,
%
%  [x1(2)-x1(1)       0       -1   0;      [t1;      [-x1(1);
%        0       x2(2)-x2(1)  -1   0;   *   t2;   =   -x2(1);
%   y1(2)-y1(1)       0        0  -1;       x0;       -y1(1);
%        0       y2(2)-y2(1)   0  -1]       y0]       -y2(1)]
%
% Let's call that A*T = B.  We can solve for T with T = A\B.
%
% Once we have our solution we just have to look at t1 and t2 to determine
% whether L1 and L2 intersect.  If 0 <= t1 < 1 and 0 <= t2 < 1 then the two
% line segments cross and we can include (x0,y0) in the output.
%
% In principle, we have to perform this computation on every pair of line
% segments in the input data.  This can be quite a large number of pairs so
% we will reduce it by doing a simple preliminary check to eliminate line
% segment pairs that could not possibly cross.  The check is to look at the
% smallest enclosing rectangles (with sides parallel to the axes) for each
% line segment pair and see if they overlap.  If they do then we have to
% compute t1 and t2 (via the A\B computation) to see if the line segments
% cross, but if they don't then the line segments cannot cross.  In a
% typical application, this technique will eliminate most of the potential
% line segment pairs.


% Input checks.
if verLessThan('matlab','7.13')
	error(nargchk(2,5,nargin)) %#ok<NCHKN>
else
	narginchk(2,5)
end

% Adjustments based on number of arguments.
switch nargin
	case 2
		robust = true;
		x2 = x1;
		y2 = y1;
		self_intersect = true;
	case 3
		robust = x2;
		x2 = x1;
		y2 = y1;
		self_intersect = true;
	case 4
		robust = true;
		self_intersect = false;
	case 5
		self_intersect = false;
end

% x1 and y1 must be vectors with same number of points (at least 2).
if sum(size(x1) > 1) ~= 1 || sum(size(y1) > 1) ~= 1 || ...
		length(x1) ~= length(y1)
	error('X1 and Y1 must be equal-length vectors of at least 2 points.')
end
% x2 and y2 must be vectors with same number of points (at least 2).
if sum(size(x2) > 1) ~= 1 || sum(size(y2) > 1) ~= 1 || ...
		length(x2) ~= length(y2)
	error('X2 and Y2 must be equal-length vectors of at least 2 points.')
end


% Force all inputs to be column vectors.
x1 = x1(:);
y1 = y1(:);
x2 = x2(:);
y2 = y2(:);

% Compute number of line segments in each curve and some differences we'll
% need later.
n1 = length(x1) - 1;
n2 = length(x2) - 1;
xy1 = [x1 y1];
xy2 = [x2 y2];
dxy1 = diff(xy1);
dxy2 = diff(xy2);


% Determine the combinations of i and j where the rectangle enclosing the
% i'th line segment of curve 1 overlaps with the rectangle enclosing the
% j'th line segment of curve 2.

% Original method that works in old MATLAB versions, but is slower than
% using binary singleton expansion (explicit or implicit).
% [i,j] = find( ...
% 	repmat(mvmin(x1),1,n2) <= repmat(mvmax(x2).',n1,1) & ...
% 	repmat(mvmax(x1),1,n2) >= repmat(mvmin(x2).',n1,1) & ...
% 	repmat(mvmin(y1),1,n2) <= repmat(mvmax(y2).',n1,1) & ...
% 	repmat(mvmax(y1),1,n2) >= repmat(mvmin(y2).',n1,1));

% Select an algorithm based on MATLAB version and number of line
% segments in each curve.  We want to avoid forming large matrices for
% large numbers of line segments.  If the matrices are not too large,
% choose the best method available for the MATLAB version.
if n1 > 1000 || n2 > 1000 || verLessThan('matlab','7.4')
	% Determine which curve has the most line segments.
	if n1 >= n2
		% Curve 1 has more segments, loop over segments of curve 2.
		ijc = cell(1,n2);
		min_x1 = mvmin(x1);
		max_x1 = mvmax(x1);
		min_y1 = mvmin(y1);
		max_y1 = mvmax(y1);
		for k = 1:n2
			k1 = k + 1;
			ijc{k} = find( ...
				min_x1 <= max(x2(k),x2(k1)) & max_x1 >= min(x2(k),x2(k1)) & ...
				min_y1 <= max(y2(k),y2(k1)) & max_y1 >= min(y2(k),y2(k1)));
			ijc{k}(:,2) = k;
		end
		ij = vertcat(ijc{:});
		i = ij(:,1);
		j = ij(:,2);
	else
		% Curve 2 has more segments, loop over segments of curve 1.
		ijc = cell(1,n1);
		min_x2 = mvmin(x2);
		max_x2 = mvmax(x2);
		min_y2 = mvmin(y2);
		max_y2 = mvmax(y2);
		for k = 1:n1
			k1 = k + 1;
			ijc{k}(:,2) = find( ...
				min_x2 <= max(x1(k),x1(k1)) & max_x2 >= min(x1(k),x1(k1)) & ...
				min_y2 <= max(y1(k),y1(k1)) & max_y2 >= min(y1(k),y1(k1)));
			ijc{k}(:,1) = k;
		end
		ij = vertcat(ijc{:});
		i = ij(:,1);
		j = ij(:,2);
	end
	
elseif verLessThan('matlab','9.1')
	% Use bsxfun.
	[i,j] = find( ...
		bsxfun(@le,mvmin(x1),mvmax(x2).') & ...
		bsxfun(@ge,mvmax(x1),mvmin(x2).') & ...
		bsxfun(@le,mvmin(y1),mvmax(y2).') & ...
		bsxfun(@ge,mvmax(y1),mvmin(y2).'));
	
else
	% Use implicit expansion.
	[i,j] = find( ...
		mvmin(x1) <= mvmax(x2).' & mvmax(x1) >= mvmin(x2).' & ...
		mvmin(y1) <= mvmax(y2).' & mvmax(y1) >= mvmin(y2).');
	
end


% Find segments pairs which have at least one vertex = NaN and remove them.
% This line is a fast way of finding such segment pairs.  We take
% advantage of the fact that NaNs propagate through calculations, in
% particular subtraction (in the calculation of dxy1 and dxy2, which we
% need anyway) and addition.
% At the same time we can remove redundant combinations of i and j in the
% case of finding intersections of a line with itself.
if self_intersect
	remove = isnan(sum(dxy1(i,:) + dxy2(j,:),2)) | j <= i + 1;
else
	remove = isnan(sum(dxy1(i,:) + dxy2(j,:),2));
end
i(remove) = [];
j(remove) = [];

% Initialize matrices.  We'll put the T's and B's in matrices and use them
% one column at a time.  AA is a 3-D extension of A where we'll use one
% plane at a time.
n = length(i);
T = zeros(4,n);
AA = zeros(4,4,n);
AA([1 2],3,:) = -1;
AA([3 4],4,:) = -1;
AA([1 3],1,:) = dxy1(i,:).';
AA([2 4],2,:) = dxy2(j,:).';
B = -[x1(i) x2(j) y1(i) y2(j)].';

% Loop through possibilities.  Trap singularity warning and then use
% lastwarn to see if that plane of AA is near singular.  Process any such
% segment pairs to determine if they are colinear (overlap) or merely
% parallel.  That test consists of checking to see if one of the endpoints
% of the curve 2 segment lies on the curve 1 segment.  This is done by
% checking the cross product
%
%   (x1(2),y1(2)) - (x1(1),y1(1)) x (x2(2),y2(2)) - (x1(1),y1(1)).
%
% If this is close to zero then the segments overlap.

% If the robust option is false then we assume no two segment pairs are
% parallel and just go ahead and do the computation.  If A is ever singular
% a warning will appear.  This is faster and obviously you should use it
% only when you know you will never have overlapping or parallel segment
% pairs.

if robust
	overlap = false(n,1);
	warning_state = warning('off','MATLAB:singularMatrix');
	% Use try-catch to guarantee original warning state is restored.
	try
		lastwarn('')
		for k = 1:n
			T(:,k) = AA(:,:,k)\B(:,k);
			[unused,last_warn] = lastwarn; %#ok<ASGLU>
			lastwarn('')
			if strcmp(last_warn,'MATLAB:singularMatrix')
				% Force in_range(k) to be false.
				T(1,k) = NaN;
				% Determine if these segments overlap or are just parallel.
				overlap(k) = rcond([dxy1(i(k),:);xy2(j(k),:) - xy1(i(k),:)]) < eps;
			end
		end
		warning(warning_state)
	catch err
		warning(warning_state)
		rethrow(err)
	end
	% Find where t1 and t2 are between 0 and 1 and return the corresponding
	% x0 and y0 values.
	in_range = (T(1,:) >= 0 & T(2,:) >= 0 & T(1,:) <= 1 & T(2,:) <= 1).';
	% For overlapping segment pairs the algorithm will return an
	% intersection point that is at the center of the overlapping region.
	if any(overlap)
		ia = i(overlap);
		ja = j(overlap);
		% set x0 and y0 to middle of overlapping region.
		T(3,overlap) = (max(min(x1(ia),x1(ia+1)),min(x2(ja),x2(ja+1))) + ...
			min(max(x1(ia),x1(ia+1)),max(x2(ja),x2(ja+1)))).'/2;
		T(4,overlap) = (max(min(y1(ia),y1(ia+1)),min(y2(ja),y2(ja+1))) + ...
			min(max(y1(ia),y1(ia+1)),max(y2(ja),y2(ja+1)))).'/2;
		selected = in_range | overlap;
	else
		selected = in_range;
	end
	xy0 = T(3:4,selected).';
	
	% Remove duplicate intersection points.
	[xy0,index] = unique(xy0,'rows');
	x0 = xy0(:,1);
	y0 = xy0(:,2);
	
	% Compute how far along each line segment the intersections are.
	if nargout > 2
		sel_index = find(selected);
		sel = sel_index(index);
		iout = i(sel) + T(1,sel).';
		jout = j(sel) + T(2,sel).';
	end
else % non-robust option
	for k = 1:n
		[L,U] = lu(AA(:,:,k));
		T(:,k) = U\(L\B(:,k));
	end
	
	% Find where t1 and t2 are between 0 and 1 and return the corresponding
	% x0 and y0 values.
	in_range = (T(1,:) >= 0 & T(2,:) >= 0 & T(1,:) < 1 & T(2,:) < 1).';
	x0 = T(3,in_range).';
	y0 = T(4,in_range).';
	
	% Compute how far along each line segment the intersections are.
	if nargout > 2
		iout = i(in_range) + T(1,in_range).';
		jout = j(in_range) + T(2,in_range).';
	end
end
end
% Plot the results (useful for debugging).
% plot(x1,y1,x2,y2,x0,y0,'ok');

function y = mvmin(x)
% Faster implementation of movmin(x,k) when k = 1.
y = min(x(1:end-1),x(2:end));
end
function y = mvmax(x)
% Faster implementation of movmax(x,k) when k = 1.
y = max(x(1:end-1),x(2:end));
end
