%
% 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. 
%
% May, 2016. Last modified September, 2021
% Copyright 2016 WindRiskTech L.L.C.
%
%------------------------------------------------------------------------
%
params    %  Load parameters
load temp eset yearstore shape vcrit vstore
%
if longmin < 0
    longmin=longmin+360;
end
if longmax < 0
    longmax=longmax+360;
end    
clf('reset')
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);
%
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;
%
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    
%
%  Read in image
%
imagename=strcat('map_',imageres);
try
    load(imagename)
    %A=imread('eocene_cont_map.jpg');
    %load eomap2  % These two lines for eocene continents background
    %displacement=0;
catch err
    errordlg('Error:  Image not available. May need to run initialize.m.')
    return
end 
%   
im=image(A);
set(im,'visible','off')
set(gcf,'visible','off')
hold on
[ay, ax, dum]=size(A);
aspect=ax/ay;
%
longmin=ax*(longmin-displacement)/360;
longmax=ax*(longmax-displacement)/360;
latmin=ay-ay*(90+latmin)/180;
latmax=ay-ay*(90+latmax)/180;
%
% 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
%
end
%
% Label for tracks
%
mmax=0;
if strcmp(mapmode,'auto')
    longmin=2.*ax;
    longmax=1;
    latmin=ay;
    latmax=1;
end    
for ntt=net:-1:nt
   if nswitch==1
      nt3s=num2str(ntt);
   else
      for jj=1:max(nsortb)
         if nsortb(jj)==nt3(ntt)
            nt3s=num2str(jj);
         end
      end   
   end   
   ntcolor=ntt-nt+1;
   clear x y v s x2 y2 v2
   ntrack=1;
   np1=0;
   np2=0;
   %
   % Find tracks, allowing for one potential wrap in longitude
   %
  for i=1:jmax(nt3(ntt))
     np1=np1+1;
     y(np1)=ay-ay.*(latbest(nt3(ntt),i)+90)./180;
     x(np1)=ax.*(longbest(nt3(ntt),i)+360-displacement)/360;
     if colorscheme < 4
         v(np1)=vbest(nt3(ntt),i);
         %vcmax=max(vcmax,v(np1));
         %vcmin=min(vcmin,v(np1));
     else
         v(np1)=pbest(nt3(ntt),i);
     end   
     if char(mapmode) == char('auto')
         if x(np1) <= ax
            longmax=max(longmax,x(np1));
            longmin=min(longmin,x(np1));
            latmax=max(latmax,y(np1));
            latmin=min(latmin,y(np1));
         end   
     end    
     if x(np1) > ax
         ntrack=2;
         np2=np2+1;
         x2(np2)=x(np1)-ax;
         y2(np2)=y(np1);
         v2(np2)=v(np1);
         np1=np1-1;
         if char(mapmode) == char('auto')
             longmax=max(longmax,x2(np2));
             longmin=min(longmin,x2(np2));
             latmax=max(latmax,y2(np2));
             latmin=min(latmin,y2(np2));
         end    
     end    
  end 
   %
   % 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:np1
      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;
      plot(xt,yt,'color',vtcolor,'linewidth',twidth)
    end  
    if ntrack == 2
    xt=zeros(1,2);
    yt=zeros(1,2);
    for i=2:np2
      xt(1)=x2(i-1);
      xt(2)=x2(i);
      yt(1)=y2(i-1);
      yt(2)=y2(i);
      vta=v2(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   
      plot(xt,yt,'color',vtcolor,'linewidth',twidth)
      sar=vtcolor;
    end  
    end
   else 
      plot(x(1:np1),y(1:np1),'color',s,'linewidth',twidth) 
      if ntrack == 2
          plot(x2,y2,'color',s,'linewidth',twidth)
      end       
   end
%
% Put arrowheads on tracks
%
   if ntrack == 1 && np1 > 1
       ua=x(np1)-x(np1-1);
       va=y(np1)-y(np1-1);
       mag=sqrt(ua.^2+va.^2);
       ua=afactor*(ax/3600)*ua/mag;
       va=afactor*(ay/1800)*va/mag;
       xq=x(np1);
       yq=y(np1);
       g1=plot_arrow(x(np1-1),y(np1-1),x(np1),y(np1),...
           'linewidth',twidth,'edgecolor',sar,'facecolor',sar,'color',sar,'headwidth',4,'headheight',4);
       xb=x(1);
       yb=y(1);
       if strcmp(startdot,'y')
        g2=plot(xb,yb,'.','color',s,'markersize',6);
       end
       if strcmp(textlabel,'y')
           g2=text(xb,yb,nt3s);
           set(g2,'fontname','gfont','fontweight','bold','fontsize',9,'color',s)
       end    
   end    
   if  ntrack == 2 && np2 > 1
       ua=x2(np2)-x2(np2-1);
       va=y2(np2)-y2(np2-1);
       mag=sqrt(ua.^2+va.^2);
       ua=afactor*(ax/3600)*ua/mag;
       va=afactor*(ay/1800)*va/mag;
       xq=x2(np2);
       yq=y2(np2);
       g1=plot_arrow(x2(np2-1),y2(np2-1),x2(np2),y2(np2),...
           'linewidth',twidth,'edgecolor',sar,'facecolor',sar,'color',sar,'headwidth',4,'headheight',4);
       xb=x2(1);
       yb=y2(1);
       if strcmp(startdot,'y')
        g2=plot(xb,yb,'.','color',s,'markersize',6);
       end 
       if strcmp(textlabel,'y')
           g2=text(xb,yb,nt3s);
           set(g2,'fontname',gfont,'fontweight','bold','fontsize',9,'color',s)
       end    
   end    
end
hold off
%
% Label map
%
xbeg=(incr-mod(displacement,incr))*ax/360;
xend=ax*(1-incr/360)+xbeg;
if xbeg == 0
    xend=ax;
end    
xtic=xbeg:incr.*ax/360:xend;
xtic(1)=xbeg+0.005.*ax;
set(gca,'xtick',xtic);
xinc=displacement+incr:incr:360+displacement;
xsize=max(size(xinc));
xtlab=cell(1,xsize);
for i=1:xsize
    if xinc(i) <= 180
        xtlab(i)={strcat(num2str(xinc(i)), 'E')};
    elseif xinc(i) > 180 && xinc(i) <360
        xtlab(i)={strcat(num2str(360-xinc(i)),'W')};
    else
        xtlab(i)={strcat(num2str(xinc(i)-360),'E')};
    end
end    
set(gca,'xticklabel',xtlab)
ytic=0:incr*ay/180:ay;
ytic(1)=0.005.*ay;
set(gca,'ytick',ytic)
yinc=90:-incr:-90;
ysize=max(size(yinc));
ytlab=cell(1,ysize);
for i=1:ysize
    if yinc(i) > 0
        ytlab(i)={strcat(num2str(yinc(i)), 'N')};
    elseif yinc(i) < 0
        ytlab(i)={strcat(num2str(-yinc(i)),'S')};
    else
        ytlab(i)={num2str(yinc(i))};
    end
end    
set(gca,'yticklabel',ytlab)
set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
set(gca,'tickdir','out')   
%
% Add grid lines to map
%
xline=zeros(max(size(ytic)),ax);
yline=zeros(max(size(ytic)),ax);
for i=1:max(size(ytic))
    xline(i,1:ax)=1:ax;
    yline(i,1:ax)=ytic(i);
end
xline=transpose(xline);
yline=transpose(yline);
hold on
line(xline,yline,'linestyle',gridlstyle,'linewidth',glinewidth,'color',gridcolor)
clear xline yline
xline=zeros(max(size(xtic)),ay);
yline=zeros(max(size(xtic)),ay);
hold off
hold on
for i=1:max(size(xtic))
    yline(i,1:ay)=1:ay;
    xline(i,1:ay)=xtic(i);
end
xline=transpose(xline);
yline=transpose(yline);
line(xline,yline,'linestyle',gridlstyle,'linewidth',glinewidth,'color',gridcolor)
hold off
%
delx=longmax-longmin;
dely=latmax-latmin;
longmin=longmin-bfac*delx;
longmin=max(longmin,1);
longmax=longmax+bfac*delx;
longmax=min(longmax,ax);
latmin=latmin-bfac*dely;
latmin=max(latmin,1);
latmax=latmax+bfac*dely;
latmax=min(latmax,ay);
delx=longmax-longmin;
dely=latmax-latmin;
delx2=aspect*dely;
dely2=delx/aspect;
if delx2 > delx
    delxf=delx2;
    delyf=dely;
else
    delxf=delx;
    delyf=dely2;
end    
longmaxn=0.5*(longmin+longmax+delxf);
longminn=0.5*(longmin+longmax-delxf);
longmax=longmaxn;
longmin=longminn;
latmaxn=0.5*(latmin+latmax+delyf);
latminn=0.5*(latmin+latmax-delyf);
latmax=latmaxn;
latmin=latminn;
longmin=floor(longmin);
longmax=ceil(longmax);
latmin=floor(latmin);
latmin=max(latmin,0);
latmax=floor(latmax);
latmax=min(latmax,0.95*ay);
longmax=min(longmax,ax); 
set(gca,'xlim',[longmin longmax],'ylim', [latmin latmax]) 
set(im,'visible','on')
set(gcf,'visible','on')
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')
xmin=longmin;
xmax=longmax;
ymin=latmin;
ymax=latmax;
%
% 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    
%
% 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 