%
% This script creates maps of multiple 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. "Prep" must be run before this can be used, 
% and the m_map routines must be installed.  November 2008
% Copyright 2009 WindRiskTech L.L.C.
% Modified October, 2020 to account for MATLAB's new 256-member colormap
% Modified August, 2022, to speed up generation of multi-color tracks
%-------------------------------------------------------------------------
%
params    %  Load parameters
monthtext={'January','February','March','April','May','June','July','August', ...
    'September','October','November','December'};
%
clf('reset')
%
if exist('shearstore','var') == 0
    shape='circ';
    load('temp.mat')
end    
fs3='sorted';
load(fs3)
[nsize,jsize]=size(latstore);
%
if strcmp(bas,'GB')
       projection=gproject;
end
if strcmp(bas,'MT')
    for i=1:nsize
        for j=1:jsize
            if longstore (i,j) > 300
                longstore(i,j)=longstore(i,j)-360; %#ok<SAGROW>
            end
        end
    end
end    
if strcmp(bas,'SA')
    for i=1:nsize
        for j=1:jsize
            if longstore (i,j) < 100
                longstore(i,j)=longstore(i,j)+360; 
            end
        end
    end
end  
%
if 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    
%
[~,jmax]=min(vstore,[],2);
jmax(vstore(:,end)~=0)=jsize;
jmax=jmax-1;
%
hold off
%
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;
%
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)=nsort(ntt);
   %
   elseif nswitch == 2
   %
   % Ordered by most intense near POI
   %
   nt3(ntt)=nsortp(ntt);
   %
   elseif nswitch == 1
   %
   % Ordered by most intense at POI
   %
   nt3(ntt)=nsortc(ntt);
   %
   elseif nswitch == 4
   %
   % Randomly ordered
   %
   nt3(ntt)=1+floor((max(nsort)-1).*rand);
   %
   end
%
   for i=1:jsize
      if latstore(nt3(ntt),i) ~= 0
         y2=latstore(nt3(ntt),i);
         %
         % Preprocess longstore to handle greenwich crossings
         %
         if longstore(nt3(ntt),1) < 30 && longstore(nt3(ntt),i) > 200
             longstore(nt3(ntt),i)=longstore(nt3(ntt),i)-360; 
         end
         if longstore(nt3(ntt),1) > 200 && longstore(nt3(ntt),i) < 100
             longstore(nt3(ntt),i)=longstore(nt3(ntt),i)+360; 
         end
         %
         x2=longstore(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
%
mfactorx=(xmax-xmin)/30;
mfactory=(ymax-ymin)/30;
%
%  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;   
%
vcmax=0;vcmin=300;
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=latstore(nt3(ntt),1:jmax(nt3(ntt)));
   x=longstore(nt3(ntt),1:jmax(nt3(ntt)));
   v=vnet(nt3(ntt),1:jmax(nt3(ntt)));
   vcmax=max(vcmax,max(v));
   vcmin=min(vcmin,min(nonzeros(v)));
   %
   % Color and plot tracks
   %
   if colorscheme == 0
        s=char(col(ntcolor));
   elseif colorscheme == 1
        s=tcolor;
   else
        s=vcolor(1,:);
   end 
   sar=s;
   %
   [ncolors,~]=size(vcolor);
   colorfrac=ncolors/64;
   if colorscheme > 1
       if colorscheme == 2
           vms=v*1852/3600;
            vss=vms;
            vss(vms<33)=32;
            vss(vms>=33 & vms<45)=39;
            vss(vms>=45 & vms<54.4)=50;
            vss(vms>=54.4 & vms<62.4)=60;
            vss(vms>=62.4 & vms<69.5)=65;
            vss(vms>=69.5)=75;
            v=vss*3600/1852;
            colormap(jet(6))
       end    
      [X,Y]=m_ll2xy(x,y);
      p = patch([X NaN], [Y NaN], [v NaN], [v NaN], 'edgecolor', 'interp','linewidth',twidth);
   else
     m_plot(x,y,'color',s,'linewidth',twidth) 
   end 
   %
% Put arrowheads on tracks
%
   mmax=jmax(nt3(ntt));
   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 is
       %  desired
       %nt3s=['  ',char(monthtext(monthstore(nt3(ntt),1))),' ',char(num2str(yearstore(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 isempty(storm)
    datea(1)=upper(datea(1));
    tit1=datea;
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',12,'fontweight','bold','interpreter','none')
%
% Add colorbar
%
if strcmpi(cbar,'y') && colorscheme > 1
    vcmin=wunitsfac*vcmin;
    vcmax=wunitsfac*vcmax;
    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
        set(get(rts,'xlabel'),'String',wunitslabel,'fontweight','bold')
    end    
    if colorscheme == 2
        aa={'TS','Cat 1','Cat 2','Cat3','Cat4','Cat 5'};
        xttemp=get(rts,'xtick');
        xti=(max(xttemp)-min(xttemp))/5;
        xt=min(xttemp):xti:max(xttemp);
        set(rts,'xtick',xt,'xticklabel',aa)
        set(get(rts,'xlabel'),'String','Saffir-Simpson','fontsize',12,'fontweight','bold')
    end    
end  
set(gca,'fontweight','bold')
%
load temp longstore  %  Reset longstore
%
% 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 