%
% 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.
% K. Emanuel, November 2008
% Copyright 2009 WindRiskTech L.L.C.
% Modified for use by AER, October, 2010
%
%------------------------------------------------------------------------
% Options:
%
transfac=0.7;           % Fraction of translation speed to add to maximum
                        % circular 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=[0.8 0.8 0.0];   % Track color, if single color used
%
twidth=0.9;             % Linewidth of tracks
%
afactor=8;              % Size of arrowheads
%
textlabel=0;            % Set to 1 if track number labels are desired
%
%imagename='equirectangular-projection7.jpg';
%
displacement=30;        % Displacement of east edge of image eastward from
                        %   the Greenwich Meridian
                        %
bfac=0.1;               % Factor governing size of boundary around map
%
incr=15;                % Lat-long increment in map background grid
%
linestyle='none';       % Background grid linestyle
%
linewidth=0.5;          % Background grid line width
%
color='g';              % Background grid color
%
%--------------------------------------------------------------------------
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');
%
%  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
%
%------------------------------------------------------------------------
reduction=0.7;  % Specified fraction of translation speed to add to 
                % circular wind speed.
%------------------------------------------------------------------------
[ut,vt,jmax] = utrans(latstore,longstore);
[nsize,jsize]=size(latstore);
net=min(net,nsize);
%
delx=15.0;
dely=5.0;
%
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;
vcolor(1:20,3)=1;
vcolor(1:20,2)=0.6;
uv=vcolor(20,:);
%
%  Read in image
%
%A=imread(imagename);  
%load(imagename)
load satimage
%
image(A)
[ay ax dum]=size(A);
aspect=ax./ay;
%
% 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 point of interest
   %
   nt3(ntt)=nsortp(ntt);
   %
   elseif nswitch == 3
   %
   % Randomly ordered
   %
   nt3(ntt)=1+floor((max(nsort)-1).*rand);
   %
   end
%
end
hold on
%
%Label for tracks
%
mmax=0;
xmin=2.*ax;
xmax=1;
ymin=ay;
ymax=1;
for ntt=net:-1:nt,
   if nswitch==1
      nt3s=num2str(ntt);
   else
      for jj=1:max(nsort),
         if nsort(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.*(latstore(nt3(ntt),i)+90)./180;
         x(np1)=ax.*(longstore(nt3(ntt),i)+360-displacement)/360;
         v(np1)=vstore(nt3(ntt),i)+transfac.*sqrt(ut(nt3(ntt),i).^2+ ...
             ut(nt3(ntt),i).^2);
         v(np1)=v(np1).*1852./3600;
         %v(np1)=pstore(nt3(ntt),i);
         if x(np1) <= ax
            xmax=max(xmax,x(np1));
            xmin=min(xmin,x(np1));
            ymax=max(ymax,y(np1));
            ymin=min(ymin,y(np1));
         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;
             xmax=max(xmax,x2(np2));
             xmin=min(xmin,x2(np2));
             ymax=max(ymax,y2(np2));
             ymin=min(ymin,y2(np2));
         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:np1,
      xt(1)=x(i-1);
      xt(2)=x(i);
      yt(1)=y(i-1);
      yt(2)=y(i);
      vt=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;
      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);
      vt=v2(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   
      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,'headwidth',3,'headheight',3);
       xb=x(1);
       yb=y(1);
       g2=plot(xb,yb,'.','color',s,'markersize',6);
       if textlabel == 1
           g2=text(xb,yb,nt3s);
           set(g2,'fontname','arial','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,'headwidth',3,'headheight',3);
       xb=x2(1);
       yb=y2(1);
       g2=plot(xb,yb,'.','color',s,'markersize',6);
       if textlabel == 1
           g2=text(xb,yb,nt3s);
           set(g2,'fontname','arial','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,'fontweight','bold')
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',linestyle,'linewidth',linewidth,'color',color)
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',linestyle,'linewidth',linewidth,'color',color)
hold off
%
delx=xmax-xmin;
dely=ymax-ymin;
xmin=xmin-bfac.*delx;
xmin=max(xmin,1);
xmax=xmax+bfac.*delx;
xmax=min(xmax,ax);
ymin=ymin-bfac.*dely;
ymin=max(ymin,1);
ymax=ymax+bfac.*dely;
ymax=min(ymax,ay);
delx=xmax-xmin;
dely=ymax-ymin;
delx2=aspect.*dely;
dely2=delx./aspect;
if delx2 > delx
    delxf=delx2;
    delyf=dely;
else
    delxf=delx;
    delyf=dely2;
end    
xmaxn=0.5.*(xmin+xmax+delxf);
xminn=0.5.*(xmin+xmax-delxf);
xmax=xmaxn;
xmin=xminn;
ymaxn=0.5.*(ymin+ymax+delyf);
yminn=0.5.*(ymin+ymax-delyf);
ymax=ymaxn;
ymin=yminn;
xmin=floor(xmin);
xmax=ceil(xmax);
ymin=floor(ymin);
ymin=max(ymin,0);
ymax=ceil(ymax);
xmax=min(xmax,ax);
set(gca,'xlim',[xmin xmax],'ylim', [ymin ymax]) 
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')
xminold=xmin;
xmaxold=xmax;
yminold=ymin;
ymaxold=ymax;