%
% This script creates a map of a tropical cyclone track, specified 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 by track number. "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:
%
% Fraction of translation speed to ad to maximum circular wind
    transfac=0.7;
% Fit cubic spline to track (y or n)?
    cspline='n';
% Time resolution in cubic spline fit (hours; maximum = 6)
    trackres=0.5;
% Line thickness (integer)
	thick=3;
% Line color (standard MATLAB options; enforce for colorcoded=0 below)
	linecolor='b';
% Line color coding by intensity (0 for none, 1 for modified Safir-Simpson,
% 2 for continuous). Modified Safir-Simpson has five categories and sub-hurricane
% category; lower limits of categories 1 and 5 are same as normal Safir-Simspon,
% but category boundaries are in equal increments of v^2. For continuous curve,
% color range is linear in v, from minimum of 15 m/s to maximum of 80 m/s.
	colorcoded=1;
% Include colorbar for intensity scale (y or n)?
    cbar='y';
% Colormap used for colorcoding of track line
	trackcolormap='jet';
% Circles plotted at 00GMT positions (y or n)?
	cplot='y';   
% Color of circles (MATLAB options)
	ccol=[1 0.1 0.7];
%   ccol=[0.2 0.2 0.2];
% Circle size (Matlab MarkerSize)
    csize=16;
% Dates plotted next to 00 GMT positions (y or n)?
	dplot='y';
% Color of date labels (MATLAB options)
	textcolor='y';
% Fontsize of date labels
    textfontsize=8;
% Factor controling distance of text from track curve (0-3)
	textdfactor=1.2;
% Name of background map image
    %imagename='earth-map-mod.jpg';  
% Displacement of east (right) edge of map eastward from Greenwich Meridian
    displacement=30;    
% Factor governing size of boundary around map
    bfac=0.1;
% Lat-long increment in map background grid
    incr=10;              
% Background grid linestyle
    linestyle='none';       
% Background grid line width
    linewidth=0.5;          
% Background grid color
    color='g';                
% Axis label font
	axisfont='arial';
% Axis label font size
	axisfontsize=12;
% Axis label font weight
    axisfontweight='bold';
%---------------------------------------------------------------------------------
%
clear daymod v vtrans vpot shear
%
lc=linecolor;
colormap('default')
%
clear x* y* u* v* date offset pres daym hour
%
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   
nttemp=input('Input track number');
%
[ut,vt,jmax] = utrans(latstore,longstore);
%
if nswitch == 1
   nt=nsortp(nttemp);
elseif nswitch == 2   
   nt=nsort(nttemp);
else
   nt=nttemp;
end   
%
ntstring=num2str(nt);
clf
%A=imread(imagename);
load satimage
%
image(A)
[ay ax dum]=size(A);
aspect=ax./ay;
[ndum,m8max]=size(vstore);
iq=jmax(nt);
iqo=iq;
%
vtrans=sqrt(ut(nt,1:iq).^2+vt(nt,1:iq).^2);
v(1:iq)=vstore(nt,1:iq)+transfac.*vtrans;
pres(1:iq)=pstore(nt,1:iq);
x(1:iq)=ax.*(longstore(nt,1:iq)+360-displacement)/360;
y(1:iq)=ay-ay.*(latstore(nt,1:iq)+90)./180;
vpot(1:iq)=vpstore(nt,1:iq);
shear(1:iq)=shearstore(nt,1:iq);
monthm(1:iq)=monthstore(nt,1:iq);
daym(1:iq)=daystore(nt,1:iq);
hour(1:iq)=hourstore(nt,1:iq);
%
ntrack=1;
mini=iq;
np2=0;
for i=1:iq,
    if x(i) > ax
        ntrack=2;
        mini=min(i,mini);
        np2=np2+1;
        x2(np2)=x(i)-ax;
        y2(np2)=y(i);
    end
end
if ntrack == 2
    dumx(1:mini-1)=x(1:mini-1);
    dumy(1:mini-1)=y(1:mini-1);
    clear x y
    x=dumx;
    y=dumy;
    clear dumx dumy
    iq=mini-1;
end    
%
for i=1:iqo,
   if monthm(i) == 12 & daym(i) > 31
      monthm(i)=1;
      daym(i)=daym(i)-31;
   end   
   %hour(i)=hourstore(nt,i);
end
vtrans=smooth(vtrans);
shear=smooth(shear);
vpot=smooth(vpot);
months={'January','February','March','April','May','June','July','August','September','October','November','December'};
colormap('default')
colormap(trackcolormap)
vcolor=colormap;
%
hold on
%
if iq > 1
n=0;
m=0;
iswitch=0;
for i=1:iq,
    if hour(i) <= 0.5
        m=m+1;
    end
end
xd=zeros(1,m);
yd=zeros(1,m);
day=zeros(1,m);
month=zeros(1,m);
date=cellstr(num2str(daym(1:m)));
m=0;
for i=1:iq,
            if hour(i) <= 0.5
               m=m+1;
               xd(m)=x(i);
               yd(m)=y(i);
               day(m)=daym(i);
               month(m)=monthm(i);
               daystring=num2str(day(m));
               date(m)=cellstr(daystring);
            end      
end
%
if cspline == 'y'
   a=size(x);
   b=a(1,2);
   bx=[1:b];
   tint=floor(6./trackres);
   tintin=1./tint;
   bxm=[1:tintin:b];
   x=spline(bx,x,bxm);
   y=spline(bx,y,bxm);
   v=spline(bx,v,bxm);
   pres=spline(bx,pres,bxm);
   vpot=spline(bx,vpot,bxm);
   shear=spline(bx,shear,bxm);
   monthm=spline(bx,monthm,bxm);
   daym=spline(bx,daym,bxm);
   hour=spline(bx,hour,bxm);
   vtrans=spline(bx,vtrans,bxm);
end   
a=size(x);
b=a(1,2);
if colorcoded > 0
    xt=zeros(1,2);
    yt=zeros(1,2);
   for i=2:b,
      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 colorcoded == 1
         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',thick)
      lc=vtcolor;
   end   
else
   plot(x,y,lc,'linewidth',thick)
end
bd=size(xd);
bd=bd(1,2);
if dplot == 'y'  
    offset=zeros(1,bd);
if bd > 1   
    ul=zeros(1,bd); 
    vl=zeros(1,bd);
   for kl=2:bd-1;
      ul(kl)=0.5.*(xd(kl+1)-xd(kl-1));
      vl(kl)=0.5.*(yd(kl+1)-yd(kl-1));
   end
   ul(1)=4.*(x(2)-x(1));
   vl(1)=4.*(y(2)-y(1));
	ul(bd)=4.*(x(b)-x(b-1));
	vl(bd)=4.*(y(b)-y(b-1));
else
   ul(1)=x(2)-x(1);
   vl(1)=y(2)-y(1);
end   
for kl=1:bd,
   offset(kl)=-0.2;
   if day(kl)>9
      offset(kl)=-0.5;
   end
end   
vabs=0.001+sqrt(ul.^2+vl.^2);
dx=15.*(ax./3200).*textdfactor.*((vl./vabs)+offset);
dy=15.*(ay./1600).*textdfactor.*(-ul./vabs);
xdb=xd+dx;
ydb=yd+dy;
g=text(xdb,ydb,date);
set(g,'fontname','arial','fontweight','bold','fontsize',textfontsize,'color',textcolor)
end
ua=x(a(1,2))-x(a(1,2)-1);
va=y(a(1,2))-y(a(1,2)-1);
uabs=sqrt(ua.^2+va.^2)+0.001;
ua=10.*(ax./3600).*ua./uabs;
va=10.*(ay./1600).*va./uabs;
xq=x(a(1,2));
yq=y(a(1,2));
g1=quiver(xq,yq,ua,va,10);
set(g1,'maxheadsize',5);
set(g1,'showarrowhead','off')
set(g1,'showarrowhead','on')
set(g1,'linewidth',thick,'color',lc);
if cplot == 'y'
   plot(xd,yd,'.','color',ccol,'MarkerSize',csize)   
end   
%
end
%
if ntrack == 2
    iq=np2;
    clear x y
    x=x2;
    y=y2;
n=0;
m=0;
iswitch=0;
for i=1:iq,
    if hour(i) <= 0.5
        m=m+1;
    end
end
xd=zeros(1,m);
yd=zeros(1,m);
day=zeros(1,m);
month=zeros(1,m);
date=cellstr(num2str(daym(1:m)));
m=0;
for i=1:iq,
            if hour(i) <= 0.5
               m=m+1;
               xd(m)=x(i);
               yd(m)=y(i);
               day(m)=daym(i);
               month(m)=monthm(i);
               daystring=num2str(day(m));
               date(m)=cellstr(daystring);
            end      
end
%
if cspline == 'y'
   a=size(x);
   b=a(1,2);
   bx=[1:b];
   tint=floor(6./trackres);
   tintin=1./tint;
   bxm=[1:tintin:b];
   x=spline(bx,x,bxm);
   y=spline(bx,y,bxm);
   v=spline(bx,v,bxm);
   pres=spline(bx,pres,bxm);
   vpot=spline(bx,vpot,bxm);
   shear=spline(bx,shear,bxm);
   monthm=spline(bx,monthm,bxm);
   daym=spline(bx,daym,bxm);
   hour=spline(bx,hour,bxm);
   vtrans=spline(bx,vtrans,bxm);
end   
a=size(x);
b=a(1,2);
if colorcoded > 0
    xt=zeros(1,2);
    yt=zeros(1,2);
   for i=2:b,
      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 colorcoded == 1
         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',thick)
      lc=vtcolor;
   end   
else
   plot(x,y,lc,'linewidth',thick)
end
bd=size(xd);
bd=bd(1,2);
if dplot == 'y'  
    offset=zeros(1,bd);
if bd > 1   
    ul=zeros(1,bd); 
    vl=zeros(1,bd);
   for kl=2:bd-1;
      ul(kl)=0.5.*(xd(kl+1)-xd(kl-1));
      vl(kl)=0.5.*(yd(kl+1)-yd(kl-1));
   end
   ul(1)=4.*(x(2)-x(1));
   vl(1)=4.*(y(2)-y(1));
	ul(bd)=4.*(x(b)-x(b-1));
	vl(bd)=4.*(y(b)-y(b-1));
else
   ul(1)=x(2)-x(1);
   vl(1)=y(2)-y(1);
end   
for kl=1:bd,
   offset(kl)=-0.2;
   if day(kl)>9
      offset(kl)=-0.5;
   end
end   
vabs=0.001+sqrt(ul.^2+vl.^2);
dx=15.*(ax./3200).*textdfactor.*((vl./vabs)+offset);
dy=15.*(ay./1600).*textdfactor.*(-ul./vabs);
xdb=xd+dx;
ydb=yd+dy;
g=text(xdb,ydb,date);
set(g,'fontname','arial','fontweight','bold','fontsize',textfontsize,'color',textcolor)
end
%
ua=x(a(1,2))-x(a(1,2)-1);
va=y(a(1,2))-y(a(1,2)-1);
uabs=sqrt(ua.^2+va.^2)+0.001;
ua=10.*(ax./3600).*ua./uabs;
va=10.*(ay./1600).*va./uabs;
xq=x(a(1,2));
yq=y(a(1,2));
g1=quiver(xq,yq,ua,va,10);
set(g1,'maxheadsize',5);
set(g1,'showarrowhead','off')
set(g1,'showarrowhead','on')
set(g1,'linewidth',thick,'color',lc);
if cplot == 'y'
   plot(xd,yd,'.','color',ccol,'MarkerSize',csize)   
end   
%
if colorcoded > 0 & cbar == 'y'
   if colorcoded == 1
      ucolor=[vcolor(1,:);vcolor(13,:);vcolor(30,:);vcolor(39,:);vcolor(49,:);vcolor(61,:)];
      aa={'TS';'1';'2';'3';'4';'5'};
      bb=1/6;
      tstring='Category';
   else
      ucolor=vcolor;
      aa=[20:10:80];
      bb=1/6;
      tstring='m/s';
   end   
   colormap(ucolor);
   rts=colorbar('horiz');
   posit=get(gca,'position');
   rr=get(rts,'position');
   xl=get(gca,'xlim');
   yl=get(gca,'ylim');
   yfrac=(yl(2)-yl(1))./(xl(2)-xl(1));
   rr2=rr;
   rr2(3)=0.75.*rr(3);
   rr2(3)=max(rr2(3),0.4);
   rr2(1)=posit(1)+0.5.*posit(3)-0.5*rr2(3);
   rr2(4)=0.015;
   %
   set(rts,'position',rr2)
   set(rts,'ticklength',[0 0])
   set(rts,'fontweight','bold')
   set(rts,'xticklabel',aa)
   set(gca,'position',posit)
   set(get(rts,'xlabel'),'String',tstring,'fontweight','bold')
end
%
end
ixlab=months(month(bd));
if month(bd) > month(1)
   ixlab=strcat(months(month(1)), ' - ', months(month(bd)));
elseif month(bd) < month(1)
      ixlab=['December - January'];
end   
ixlab=char(ixlab);
titstring=['Track number ' ntstring ',   ' ixlab ];
titstring2=['Track number ' ntstring  ];
if isempty(storm)
     tit1=[datea ' ' titstring2];
else
    tit1=[storm(1:2) storm(3:4) ' ' datea(1:4) ' ' datea(5:6) ' ' datea(7:8)...
    ' ' datea(9:10) ' GMT'];
    tit2=titstring2;
    tit1=[tit1 '  ' tit2];
end    
tit3={tit1;ixlab};
title(tit3,'fontsize',12,'fontweight','bold','interpreter','none')
%
% 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
%
xmin=floor(min(x));
xmax=ceil(max(x));
ymin=floor(min(y));
ymax=ceil(max(y));
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);
ymax=ceil(ymax);
set(gca,'xlim',[xmin xmax],'ylim', [ymin ymax]) 
%
pause
%
iq=a(1,2);
clf
daytime=daym+hour./24;
daymod=daytime;
delday=max((daytime(2)-daytime(1)),(daytime(3)-daytime(2)));
iswitch=0;
for i=2:iqo,
   if daytime(i) < daymod(i-1)
      daymod(i)=daymod(i-1)+delday;
      if iswitch == 0
         iswitch=1;
         daymax=floor(daytime(i-1));
      end   
   end
end   
%vtrans=smooth(vtrans);
%vpot=smooth(vpot);
shear=smooth(shear);
h=plot(daymod,v,daymod,vtrans,daymod,vpot,daymod,shear);
set(h,'linewidth',3)
set(gca,'fontweight','bold','fontsize',12)
legend('V_{max}','V_{translation}','V_{pot}','Shear',0)
if monthm(1) ~= monthm (iq)
  jg=get(gca,'xticklabel');
  hj=str2num(char(jg));
  %set(gca,'xtickmode','auto');
  %set(gca,'xtickmode','manual');
  asiz=size(hj);
  asize=asiz(1,1);
  iswitch=0;
  for iss=1:asize,
     if hj(iss) >daymax
        hj(iss)=hj(iss)-daymax;
     end
  end
  hjs=cellstr(num2str(hj));
  set(gca,'xticklabel',hjs);
  set(gca,'fontweight','bold','fontsize',12)
end   
xlabel(ixlab,'fontweight','bold','fontsize',12)
ylabel('Maximum surface wind speed (knots)','fontweight','bold','fontsize',12)
title(tit1,'fontsize',12,'fontweight','bold','interpreter','none')
%
pause
%
clf
cla
pp=1;
while pp == 1
h=plot(daymod,pres);
set(h,'linewidth',3)
set(gca,'fontweight','bold','fontsize',12)
if monthm(1) ~= monthm (iq)
  jg=get(gca,'xticklabel');
  hj=str2num(char(jg));
  %set(gca,'xtickmode','auto');
  %set(gca,'xtickmode','manual');
  asiz=size(hj);
  asize=asiz(1,1);
  iswitch=0;
  for iss=1:asize,
     if hj(iss) >daymax
        hj(iss)=hj(iss)-daymax;
     end
  end
  hjs=cellstr(num2str(hj));
  set(gca,'xticklabel',hjs);
  set(gca,'fontweight','bold','fontsize',12)
end   
xlabel(ixlab,'fontweight','bold','fontsize',12)
ylabel('Central pressure (hPa)','fontweight','bold','fontsize',12)
title(tit1,'fontsize',12,'fontweight','bold','interpreter','none')
pause 
break
end
h=plot(daymod,shear);
set(h,'linewidth',3)
set(gca,'fontweight','bold','fontsize',12)
if monthm(1) ~= monthm (iq)
  jg=get(gca,'xticklabel');
  hj=str2num(char(jg));
  asiz=size(hj);
  asize=asiz(1,1);
  iswitch=0;
  for iss=1:asize,
     if hj(iss) >daymax
        hj(iss)=hj(iss)-daymax;
     end
  end
  hjs=cellstr(num2str(hj));
  set(gca,'xticklabel',hjs);
  set(gca,'fontweight','bold','fontsize',12)
end   
xlabel(ixlab,'fontweight','bold','fontsize',12)
ylabel('850 - 250 hPa shear (knots)','fontweight','bold','fontsize',12)
title(tit1,'fontsize',12,'fontweight','bold','interpreter','none')