%
% This script creates maps of historical tropical cyclone tracks, ordered by
% maximum lifetime intensity, by maximum wind speed at a point of interest
% or the intersection points of tracks with a specified set of line
% segments, or randomly. 
%
% September, 2021
% Copyright 2021 WindRiskTech L.L.C.
%
%------------------------------------------------------------------------
%
clf('reset')
params    %  Load parameters
load temp eset yearstore shape vcrit vstore
%
bas='gb';
bestproc
[nsize,jsize]=size(vbest); 
%
% Convert wind speed units and define S-S category boundaries
%
vbest=vbest*wunitsfac;
dc(1)=1;dc(2)=64;dc(3)=83;dc(4)=96;dc(5)=113;dc(6)=137;
dc=dc*wunitsfac;
vcmax=10*round(0.1*wunitsfac*180);
vcmin=10*round(0.1*wunitsfac*20);
vmax=max(vbest,[],2);
%
if strcmp(bas,'GB')
       projection=gproject;
end
if strcmp(bas,'MT')
    for i=1:nsize
        for j=1:jsize
            if longbest (i,j) > 300
                longbest(i,j)=longbest(i,j)-360; %#ok<SAGROW>
            end
        end
    end
end    
if strcmp(bas,'SA')
    for i=1:nsize
        for j=1:jsize
            if longbest (i,j) < 100
                longbest(i,j)=longbest(i,j)+360; 
            end
        end
    end
end  
%
if exist('storm','var') == 0 || isempty(storm)
    prompt = {'Order by maximum at POI (1), maximum near POI (2), lifetime maximum (3), or random (4)?', ...
        'Starting track number ', ...
        'Ending track number (hit return to use whole set)'};
    dlg_title = 'Track Parameters';
    num_lines = 1;
    answer = inputdlg(prompt,dlg_title,num_lines);
    if isempty(answer) || isempty(answer{1})
        return
    end    
    nswitch=str2double(answer{1});
    if (strcmp(shape,'poly') || city_radius >= 2000) && nswitch == 1
        warndlg('Warning: Winds at POI not meaningful for this set; ordering by maximum near POI')
        nswitch=2;
    end        
    nt=str2double(answer{2});
    net=str2double(answer{3});
else    
   nswitch=3;
end   
if isnan(net)
    net=nsize;
end  
net=min(net,nsize);
%
[~,jmax]=min(vbest,[],2);
jmax(vstore(:,end)~=0)=jsize;
jmax=jmax-1;
jmax=max(jmax,3);
%
hold off
%
clear col
col=cell(1,2000);
col(1)=cellstr('r');
col(2)=cellstr('b');
col(3)=cellstr('g');
col(4)=cellstr('c');
col(5)=cellstr('m');
%
for i=6:2000
   col(i)=col(i-5);
end   
colormap(trackcolormap);
vcolor=colormap;
[cosize,~]=size(vcolor);
da=0.15*cosize;
dcolor=da+(cosize-da)*(dc-dc(1))/(dc(end)-dc(1));
dcolor=round(dcolor);
if colorscheme == 4
    vcolor=flipud(vcolor);
end    
%
ymin=90;
ymax=-90;
xmin=360;
xmax=0;
%
% Determine ordering of tracks
%
nt3=zeros(1,net);
rng('shuffle')   %  Set seed of random number generator based on time
for ntt=nt:net
   if nswitch == 3
   %
   % Ordered by most intense over lifetime
   %
   nt3(ntt)=nsortb(ntt);
   %
   elseif nswitch == 2
   %
   % Ordered by most intense near POI
   %
   nt3(ntt)=nsortpb(ntt);
   %
   elseif nswitch == 1
   %
   % Ordered by most intense at POI
   %
   nt3(ntt)=nsortcb(ntt);
   %
   elseif nswitch == 4
   %
   % Randomly ordered
   %
   nt3(ntt)=1+floor((max(nsortb)-1).*rand);
   %
   end
   %
   for i=1:jsize
      if latbest(nt3(ntt),i) ~= 0
         y2=latbest(nt3(ntt),i);
         %
         % Preprocess longbest to handle Greenwich crossings
         %
         if longbest(nt3(ntt),1) < 30 && longbest(nt3(ntt),i) > 200
             longbest(nt3(ntt),i)=longbest(nt3(ntt),i)-360; 
         end
         if longbest(nt3(ntt),1) > 200 && longbest(nt3(ntt),i) < 100
             longbest(nt3(ntt),i)=longbest(nt3(ntt),i)+360; 
         end
         %
         x2=longbest(nt3(ntt),i);
         xmin=min(xmin,x2);
         xmax=max(xmax,x2);
         ymin=min(ymin,y2);
         ymax=max(ymax,y2);
      end
   end   
end  
ymax=ymax+dellat;
ymin=ymin-dellat;
xmax=xmax+dellong;
xmin=xmin-dellong;
%
if strcmp(mapmode,'manu')
    xmin=longmin;
    xmax=longmax;
    ymin=latmin;
    ymax=latmax;
end
%
%  Specify map projection (see m_map documentation)
%
m_proj(projection,'long',[xmin xmax], ...
    'lat',[ymin ymax]);
%
%  Low resolution coastline (faster)
%
%m_coast('patch',landcolor,'edgecolor','none');
%
%  Higher resolution coastlines:
%
if strcmp(bas,'SA')
    m_coast('patch',landcolor,'edgecolor','none');
else    
    m_gshhs_l('patch',landcolor,'edgecolor',stcolor);
end    
%
if char(pstates) == char('y')
    m_states(xmin, xmax, ymin, ymax,'color',stcolor)
end   
%
hold on
mmax=0;   
%
for ntt=net:-1:nt
   if nswitch == 2
        if vmax(nt3(ntt)) < lpthresh
            break
        end
   end    
   nt3s=num2str(nt3(ntt));
   ntcolor=ntt-nt+1;
   clear x y s
   % 
   y=latbest(nt3(ntt),1:jmax(nt3(ntt)));
   x=longbest(nt3(ntt),1:jmax(nt3(ntt)));
   v=vbest(nt3(ntt),1:jmax(nt3(ntt)));
   %
   % Color and plot tracks
   %
   if colorscheme == 0
        s=char(col(ntcolor));
   elseif colorscheme == 1
        s=tcolor;
   elseif colorscheme == 2 || colorscheme == 3
        s=vcolor(1,:);
   elseif colorscheme == 4
        s=vcolor(cosize,:);
   end 
   sar=s;
   if colorscheme > 1
    xt=zeros(1,2);
    yt=zeros(1,2);
    for i=2:jmax(nt3(ntt))
      xt(1)=x(i-1);
      xt(2)=x(i);
      yt(1)=y(i-1);
      yt(2)=y(i);
      vta=v(i-1);
      if colorscheme == 2
         if vta < dc(2)
            vtcolor=vcolor(dcolor(1),:);
         elseif vta >= dc(2) && vta < dc(3)
            vtcolor=vcolor(dcolor(2),:);
         elseif vta >= dc(3) && vta < dc(4)
            vtcolor=vcolor(dcolor(3),:);
         elseif vta >= dc(4) && vta < dc(5)
            vtcolor=vcolor(dcolor(4),:);
         elseif vta >= dc(5) && vta < dc(6)
            vtcolor=vcolor(dcolor(5),:);
         else
            vtcolor=vcolor(dcolor(6),:);
         end
      else   
         if colorscheme == 3 
            vtr=round(cosize*(vta-vcmin)/(vcmax-vcmin));
            vtr=max(vtr,1);
         else
            vtr=round(0.5.*cosize*(vta-900)/64);
         end   
         vtr=min(vtr,cosize);
         vtr=max(vtr,1);
         vtcolor=vcolor(vtr,:);
      end   
      sar=vtcolor;
      m_plot(xt,yt,'color',vtcolor,'linewidth',twidth)
    end  
   else
     m_plot(x,y,'color',s,'linewidth',twidth) 
   end 
   %
   % Put arrowheads on tracks
   %
   mmax=jmax(nt3(ntt));mmax=max(mmax,2);
   ua=x(mmax)-x(mmax-1);
   va=y(mmax)-y(mmax-1);
   xq=x(mmax-1);
   yq=y(mmax-1);
   scale=afactor/sqrt(ua^2+va^2);
   ua=ua*scale;
   va=va*scale;
   g1=m_quiver(xq,yq,ua,va,0);
   set(g1,'maxheadsize',4.5);
   set(g1,'showarrowhead','off')
   set(g1,'showarrowhead','on')
   set(g1,'linewidth',twidth,'color',sar);
   xb=x(1);
   yb=y(1);
   if strcmp(startdot,'y')
        g2=m_plot(xb,yb,'.','color',s,'markersize',12);
   end     
   if strcmp(textlabel,'y')
       %  Undo comment on next line if date labels at start of tracks are
       %  desired
       %nt3s=['  ',char(monthtext(monthbest(nt3(ntt),1))),' ',char(num2str(yearbest(nt3(ntt))))];
       g2=m_text(xb,yb,nt3s);
       set(g2,'fontname','arial','fontweight','bold','fontsize',9,'color',s)
   end    
end
hold off
m_grid('fontname',axisfont,'fontsize',axisfontsize,'fontweight',axisfontweight, ...
    'linestyle',gridline,'backgroundcolor',oceancolor)
%
if exist('storm','var') == 0 || isempty(storm)
    %datea(1)=upper(datea(1));
    %tit1=datea;
    tit1=eset;
else
    tit1=[storm(1:2) storm(3:4) ' ' datea(1:4) ' ' datea(5:6) ' ' datea(7:8)...
    ' ' datea(9:10) ' GMT'];
end    
%title(tit1,'fontsize',gfontsize,'fontweight','bold','interpreter','none')
%
% Add colorbar
%
if strcmpi(cbar,'y') && colorscheme > 1
    intv=20;
    if strcmpi(wunits,'kmh')
        intv=40;
    end        
    vcmin=intv*round(vcmin/intv);
    vcmax=intv*round(vcmax/intv);
    rts=colorbar('horiz');
    if colorscheme == 3
        aa=vcmin:intv:vcmax;
        bb=(aa-vcmin+1)./(vcmax-vcmin+1);
        set(rts,'xtick',bb,'xticklabel',aa)
        set(get(rts,'xlabel'),'String',wunitslabel,'fontweight','bold')
    end    
    if colorscheme == 2
        ucolor=[vcolor(dcolor(1),:);vcolor(dcolor(2),:);vcolor(dcolor(3),:); ...
          vcolor(dcolor(4),:);vcolor(dcolor(5),:);vcolor(dcolor(6),:)];
        aa={'TS','Cat 1','Cat 2','Cat3','Cat4','Cat 5',' '};
        bb=0:1/6:1;
        set(rts,'xtick',bb,'xticklabel',aa)
        set(get(rts,'xlabel'),'String','Saffir-Simpson','fontweight','bold')
        colormap(ucolor)
    end    
end    
set(gca,'fontweight','bold')
%
% Add button to access documentation of script
%
if strcmp(Docdisp,'y')
    tempname=mfilename;
    mname = strcat('graphics_docs/html/',tempname,'_doc.html');
    comstr=['open ',mname,';delete(hButton)'];
    labelStr = '<html><center><a href="">Doc';
    hButton = uicontrol('string',labelStr,'pos',[10,10,50,20],'callback',comstr); 
end
%
% Clean up
%
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 