%
% 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.  K. Emanuel, November 2008
%  Copyright 2009 WindRiskTech L.L.C.
% Last modified January, 2011
%
%-------------------------------------------------------------------------
% Options:
%
transfac=0.7;           % Fraction of translation speed to add to maximum
                        % circular wind speed
                        %
vthresh=0;              % Threshold liftime peak wind speed if tracks are
                        % ordered by lifetime maximum wind speed
%
colorscheme=3;          % Set to zero for alternating track colors; 1 for 
                        % single color, 2 for modified Saffir-Simpson
                        % scale, or 3 for continuous color coded by wind
                        % speed
%
tcolor='b';             % Track color, if single color used
%
twidth=0.9;             % Linewidth of tracks
%
afactor=2;              % Size of arrowheads
%
textlabel=0;            % Set to 1 if track number labels are desired
%
projection='Miller Cylindrical';  % Map projection
%
pstates='y';            % Plot political boundaries
% 
stcolor=[.4 .4 .4];     % Color of political boundaries (MATLAB color options)
%
%  Delx and dely are the increments of longitude and latitude the specify
%  how far the maps extend beyond the outer most limits of the tracks
%
delx=15.0;
dely=5.0;
%
%------------------------------------------------------------------------
clf
%
shape='circ';
load('temp.mat')
fs3=['sorted'];
load(fs3)
%
if isempty(storm)
   nswitch=input('Order by maximum at POI (1), lifetime maximum (2), or random (3)?');
else    
   nswitch=2;
end   
nt=input('Starting track number');
net=input('Ending track number');
%
[ut,vt,jmax] = utrans(latstore,longstore);
vnet=vstore+transfac.*sqrt(ut.^2+vt.^2);
vmax=max(vnet,[],2);
%
hold off
%
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('default');
vcolor=colormap;
%
ymin=90;
ymax=-90;
xmin=360;
xmax=0;
%
% Determine ordering of tracks
%
for ntt=nt:net,
   if nswitch == 2
   %
   % Ordered by most intense over lifetime
   %
   nt3(ntt)=nsort(ntt);
   %
   elseif nswitch == 1
   %
   % Ordered by most intense at city
   %
   nt3(ntt)=nsortp(ntt);
   %
   elseif nswitch == 3
   %
   % Randomly ordered
   %
   nt3(ntt)=1+floor((max(nsort)-1).*rand);
   %
   end
%
   for i=1:400,
      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) > 300
             longstore(nt3(ntt),i)=longstore(nt3(ntt),i)-360;
         end
         if longstore(nt3(ntt),1) > 200 & longstore(nt3(ntt),i) < 60
             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+dely;
ymin=ymin-dely;
xmax=xmax+delx;
xmin=xmin-delx;
%
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',[.45 .6 .45],'edgecolor','none');
%
%  Higher resolution coastlines:
%
m_gshhs_l('patch',[.45 .6 .45],'edgecolor','none');
%
if char(pstates) == char('y')
    m_states(xmin, xmax, ymin, ymax,'color',stcolor)
end    
%
set(gca,'color',[.95 .95 1])
%
hold on
mmax=0;
for ntt=net:-1:nt,
    if nswitch == 2
        if vmax(nt3(ntt)) < vthresh
            break
        end
    end    
   nt3s=num2str(nt3(ntt));
   ntcolor=ntt-nt+1;
   clear x y s
   x=zeros(1,jmax(nt3(ntt)));
   y=zeros(1,jmax(nt3(ntt)));
   v=zeros(1,jmax(nt3(ntt)));
   %
   % Find tracks
   %
   for i=1:jmax(nt3(ntt)),
      if latstore(nt3(ntt),i) ~= 0 
         y(i)=latstore(nt3(ntt),i);
         x(i)=longstore(nt3(ntt),i);
         v(i)=vstore(nt3(ntt),i)+transfac.*sqrt(ut(nt3(ntt),i).^2+ ...
             ut(nt3(ntt),i).^2);
      end
   end   
   %
   % Color and plot tracks
   %
   if colorscheme == 0
        s=char(col(ntcolor));
   elseif colorscheme == 1
        s=tcolor;
   else
        s=vcolor(1,:);
   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);
      %vt=0.5.*(v(i-1)+v(i));
      vt=(1852./3600).*v(i-1);
      if colorscheme == 2
         if vt < 33
            vtcolor=vcolor(1,:);
         elseif vt >= 33 & vt < 45
            vtcolor=vcolor(13,:);
         elseif vt >= 45 & vt < 54.4
            vtcolor=vcolor(30,:);
         elseif vt >= 54.4 & vt < 62.4
            vtcolor=vcolor(39,:);
         elseif vt >= 62.4 & vt < 69.5
            vtcolor=vcolor(49,:);
         else
            vtcolor=vcolor(61,:);
         end
      else   
         vtr=round(vt)-15;
         if vtr > 64
            vtr=64;
         end
         if vtr < 1
            vtr=1;
         end   
         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));
   ua=x(mmax)-x(mmax-1);
   va=y(mmax)-y(mmax-1);
   xq=x(mmax-1);
   yq=y(mmax-1);
   scale=1./sqrt(ua.^2+va.^2);
   ua=ua.*5.*scale;
   va=va.*5.*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);
   g2=m_plot(xb,yb,'.','color',s,'markersize',6);
   if textlabel == 1
       g2=m_text(xb,yb,nt3s);
       set(g2,'fontname','arial','fontweight','bold','fontsize',9,'color',s)
   end
end
hold off
m_gridm('fontname','arial','fontsize',11,'fontweight','bold','linestyle','none')
if isempty(storm)
    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')
