%
% 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. This script requires installation of the "m_map" routines.
% K. Emanuel, November 2008
% Copyright 2009 WindRiskTech L.L.C.
% Last Modified January, 2011
%
%-------------------------------------------------------------------------
% Options:
%
% Fraction of translation speed to add to circular winds
    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='k';
% Fontsize of date labels
    textfontsize=8;
% Factor controling distance of text from track curve (0-3)
	textdfactor=1.2;
% Axis label font
	axisfont='arial';
% Axis label font size
	axisfontsize=12;
% Axis label font weight
    axisfontweight='bold';
% Grid line style ('-',':','-.','--', or 'none')
    gridline='none';
% Color of land (only effective when no topography contoured)
    landcolor=[0.45 0.6 0.45];
   %landcolor=[0.6 0.6 0.6];
% Ocean color 
    oceancolor=[0.95 0.95 1];
  %oceancolor=[0.99 0.99 1];
% Ocean color - land cover swap (y or n). This is sometimes necessary when 
% program produces wrong coloring.
	swap='n';
% Plot ../chip/data/states (y or n). Need to download political boundaries from 
%  http://www.maproom.psu.edu/cgi-bin/ian/points/index.cgi
%  to use this
    pstates='y';
% Color of state boundaries (MATLAB color options)
    stcolor=[.4 .4 .4];
% Coastline resolution(0 for none, 1 for very low, 2 for low, 3 for medium,  
% 4 for high. Note that higher resolution takes more memory and CPU and 
% slows down the plot. You must separately download high resolution data bases
% for choices 3 and 4.)
	cres=2;
% Color of shoreline ('none' or MATLAB color options)
	shorecolor='none';
% Resolution of topography/bathymetry (1 for low, 2 for high. Note that you
% must separately download high reoslution data base for choice 2.)
	tres=1;
% Contour bathymetry (y or n)?
    cbath='n';
% Filled bathymetry (y or n)?
	bfill='n';
% Colormap for bathymetry
	bathcolor='jet';
% Bathymetry contour levels 
	blevels=[-8000:500:-500];
% Contour topography (y or n)?
    ctopo='n';
% Filled topography (y or n)?
    tfill='y';
% Colormap for topography
	topocolor='summer';
% Topography contour levels
    tlevels=[250 500 750 1000];
% Map projection (Options include Mercator, Miller cylindrical, Equidistant
%  cylindrical, Transverse Mercator, Sinusoidal, Gall-Peters, 
%  Lambert conformal conic, Albers equal-area conic, Hammer-Aitoff,
%  and Mollweide. Other options are possible, but the call statement then
%  has to be changed; see M_Map documentation.)
    projection='equidistant cylindrical';
% Mode of determining lat-long bounds of map ('auto' or 'manu')
    mapmode='auto';
% Width of latitude and longitude buffers relative to storm track limits (degrees);
% affects 'auto' mode only
	dellat=5;
    dellong=15;
% Longitude and latitude bounds ('manu' mode only, longitude 0 to 360, 
% latitude -90 to 90)
   longmin=-80;
   longmax=-60;
   latmin=36;
   latmax=56;
   %
   if longmin < 0
      longmin=longmin+360;
   end
   if longmax < 0
      longmax=longmax+360;
   end   
%---------------------------------------------------------------------------------
%
clear daymod v vtrans vpot shear
if swap == 'y'
   dummy=landcolor;
   landcolor=oceancolor;
   oceancolor=dummy;
end   
lc=linecolor;
ld=landcolor;
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 track number (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
xmin=longmin;
xmax=longmax;
ymin=latmin;
ymax=latmax;
[ndum,m8max]=size(vstore);
iq=jmax(nt);
%
vtrans(1:iq)=transfac.*sqrt(ut(nt,1:iq).^2+vt(nt,1:iq).^2);
v(1:iq)=vstore(nt,1:iq)+vtrans;
pres(1:iq)=pstore(nt,1:iq);
x(1:iq)=longstore(nt,1:iq);
y(1:iq)=latstore(nt,1:iq);
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);
for i=1:iq,
    if x(1) > 200 & x(i) < 50
        x(i)=x(i)+360;
    end
end    
%
for i=1:iq,
   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'};
transcolor=0.5.*(landcolor+oceancolor);
colormap('default')
colormap(trackcolormap)
vcolor=colormap;
if cbath == 'n' & ctopo == 'y'
   colormap(topocolor)
else
   colormap(bathcolor)
end   
if cres == 0 & ctopo == 'y' & tfill == 'y'
   uc=colormap;
   topocolor=[landcolor;uc];
   colormap(topocolor); 
end   
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 mapmode == 'auto'
   delcrit=15;
   if (max(max(x))-min(min(x)))<delcrit
      dellong=dellong+delcrit-max(max(x))+min(min(x));
   end   
   xmin=min(min(x))-dellong;
   xmax=max(max(x))+dellong;
   ymin=min(min(y))-dellat;
   ymax=max(max(y))+dellat;
end   
m_proj(projection,'long',[xmin xmax], ...
   'lat',[ymin ymax]);
set(gca,'color',oceancolor)
if cbath == 'y'
   if bfill == 'y'
      if tres == 2
         m_tbase('contourf',blevels);
	   else
         m_elev('contourf',blevels);
      end   
   else
      if tres == 2
         m_tbase('contour',blevels);
	   else
         m_elev('contour',blevels);
      end   
   end  
end   
hold on
if cres == 1
      m_coast('patch',ld,'edgecolor',shorecolor);
elseif cres == 2  
      m_gshhs_l('patch',ld,'edgecolor',shorecolor);
elseif cres == 3  
      m_gshhs_i('patch',ld,'edgecolor',shorecolor);
elseif cres == 4   
      m_gshhs_h('patch',ld,'edgecolor',shorecolor);
end   
if pstates == 'y'
   m_states(xmin,xmax,ymin,ymax,'color',stcolor)
end
if ctopo == 'y'
   if tfill == 'y'
      if tres == 2
         m_tbase('contourf',tlevels);
      else   
         m_elev('contourf',tlevels);
      end   
   else
      if tres == 2
         m_tbase('contour',tlevels);
      else   
         m_elev('contour',tlevels);
      end   
   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   
      m_plot(xt,yt,'color',vtcolor,'linewidth',thick)
      lc=vtcolor;
   end   
else
   m_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);
mfactorx=(xmax-xmin)./30;
mfactory=(ymax-ymin)./30;
dx=mfactorx.*textdfactor.*((vl./vabs)+offset);
dy=mfactory.*textdfactor.*(-ul./vabs);
xdb=xd+dx;
ydb=yd+dy;
g=m_text(xdb,ydb,date);
set(g,'fontname','arial','fontweight','bold','fontsize',textfontsize,'color',textcolor)
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  ];
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;
scalef=sqrt(mfactorx.^2+mfactory.^2);
ua=2.*ua.*scalef./uabs;
va=2.*va.*scalef./uabs;
xq=x(a(1,2));
yq=y(a(1,2));
g1=m_quiver(xq,yq,ua,va,0);
set(g1,'maxheadsize',4.5);
set(g1,'showarrowhead','off')
set(g1,'showarrowhead','on')
set(g1,'linewidth',thick,'color',lc);
if cplot == 'y'
   m_plot(xd,yd,'.','color',ccol,'MarkerSize',csize)   
end   
m_gridm('fontname',axisfont,'fontsize',axisfontsize,'fontweight',axisfontweight,'linestyle',gridline)
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')
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,'xaxislocation','bottom')
   set(rts,'fontweight','bold')
   set(rts,'xtick',[0:bb:1],'xticklabel',aa)
   set(gca,'position',posit)
   set(get(rts,'xlabel'),'String',tstring,'fontweight','bold')
end
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:iq,
   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));
  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));
  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')

