% Plots the tracks and intensities of events from best track data sets.
% The latter must be available in matlab binary format. 
% Last modified August, 2021
%--------------------------------------------------------------------------
% Options:
%
% Velocities in 'knots', 'ms' (meters per second), or 'kmh' (kilometers/hr)
   units='knots';      
% Track line thickness (integer)
   thick=3;
% Fit cubic spline to track (y or n)?
   cspline='y';
% Time resolution in cubic spline fit (hours; maximum = 6)
   trackres=0.5;
% Track line color (standard MATLAB options; enforce for colorcoded=0 below)
   linecolor='b';
% Track 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=2;
% Colormap used for colorcoding of track line
   trackcolormap='jet';
% Include colorbar for intensity scale (y or n)?
   cbar='y';
% Circles plotted at 00GMT positions (y or n)?
   cplot='y';   
% Color of circles (MATLAB options)
   ccol=[1 0.1 0.7];
% 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';
% Font for date labels
   dfont='arial';
% Fontsize of date labels
   textfontsize=8;
% Factor controlling 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 
   landcolor=[0.45 0.6 0.45];
% Ocean color 
   oceancolor=[0.8 0.8 0.9];
% Ocean color - land cover swap (y or n). This is sometimes necessary when 
%         program produces wrong coloring.
   swap='n';
% Plot states (y or n)
   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=3;
% 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='y';
% Colormap for bathymetry
	bathcolor='gray';
% 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 2000 3000];
% 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='lambert';
% 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=2;
   dellong=5;
% Longitude and latitude bounds ('manu' mode only, longitude 0 to 360, 
% latitude -90 to 90)
   longmin=250;
   longmax=325;
   latmin=14;
   latmax=48;
%---------------------------------------------------------------------------------
%
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;
%
% Specify storm basin (at, ep, wp, io, sh)
%
basin=input('Input basin (al, ep, wp, io, or sh)','s');
%
if strcmp(basin,'al')
   basin='at';
end   
file=strcat('best_tracks/',basin,'tracks');
atemp=units;
clear x* y* u* v* date offset
units=atemp;
clear atemp
xmin=longmin;
xmax=longmax;
ymin=latmin;
ymax=latmax;
load(file)
q=size(daym);
iq=q(1,1);
year=input('Year of storm');
clear temp
temp=input('Storm number or name','s');
clear ns
clf
ns=str2num(temp); %#ok<ST2NM>
ns0=0;
if isempty(ns)
    for i=1:iq
        if strcmpi(temp,strtrim(stname(i,:))) && yeari(i) == year
            ns=nsi(i);
            ns0=ns;
        end
    end
end    
if isempty(ns) 
    disp(' ')
    disp('Storm not found')
    return
end
ystring=num2str(year);
ystringp=num2str(year+1);
ystringm=num2str(year-1);
nstring=num2str(ns);
months={'January','February','March','April','May','June','July','August','September','October','November','December'};
transcolor=0.5.*(landcolor+oceancolor);
colormap('default')
colormap(trackcolormap)
vcolor=colormap;
[cosize,~]=size(vcolor);
dc1=1;dc2=round(cosize*13/64);dc3=round(cosize*30/64);dc4=round(cosize*39/64);dc5=round(cosize*49/64);
dc6=round(cosize*61/64);
if strcmp(cbath,'n') && strcmp(ctopo,'y')
   colormap(topocolor)
else
   colormap(bathcolor)
end   
if cres == 0 && strcmp(ctopo,'y') && strcmp(tfill,'y')
   uc=colormap;
   topocolor=[landcolor;uc];
   colormap(topocolor); 
end   
n=0;
m=0;
iswitch=0;
for i=1:iq
   if year == yeari(i) && ns == nsi(i)
      istorm=i;
      clear stormname
      if min(size(stname)) == 1
          stormname=stname(i);
      else    
        stormname=strtrim(stname(i,:)); 
      end  
      x=zeros(1,120);
      y=zeros(1,120);
      v=zeros(1,120);
      xd=zeros(1,120);
      yd=zeros(1,120);
      hour=zeros(1,120);
      day=zeros(1,120);
      month=zeros(1,120);
      date=cell(120);
      for j=1:120
         if latm(i,j)~=0
            iswitch=1;
            n=n+1;
            x(n)=longm(i,j);
            if strcmp(basin,'at') && x(n) < 100
               x(n)=x(n)+360;
            end   
            y(n)=latm(i,j);
            v(n)=vsm(i,j);
            if hourm(i,j)<2
               m=m+1;
               xd(m)=x(n);
               yd(m)=y(n);
               hour(m)=hourm(i,j);
               day(m)=daym(i,j);
               month(m)=monthm(i,j);
               daystring=num2str(day(m));
               date(m)=cellstr(daystring); 
            end      
         end
      end
   end
end
x(n+1:end)=[];
y(n+1:end)=[];
v(n+1:end)=[];
xd(m+1:end)=[];
yd(m+1:end)=[];
hour(m+1:end)=[];
day(m+1:end)=[];
month(m+1:end)=[];
date(m+1:end)=[];
if iswitch == 0
   disp('No data found for this storm')
	return   
end
if strcmp(mapmode,'auto')
   delcritlong=15;
   delcritlat=10;
   if (max(max(x))-min(min(x)))<delcritlong
      dellong=dellong+delcritlong-max(max(x))+min(min(x));
   end   
   if (max(max(y))-min(min(y)))<delcritlat
      dellat=dellat+delcritlat-max(max(y))+min(min(y));
   end   
   xmin=min(min(x))-dellong;
   xmax=max(max(x))+dellong;
   ymin=min(min(y))-dellat;
   ymax=max(max(y))+dellat+0.5;
end   
m_proj(projection,'long',[xmin xmax], ...
   'lat',[ymin ymax]);
axes('color',oceancolor)
if strcmp(cbath,'y')
   if strcmp(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 strcmp(pstates,'y')
    m_states(xmin,xmax,ymin,ymax,'color',stcolor)      
end
if strcmp(ctopo,'y')
   if strcmp(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 strcmp(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);
end   
a=size(x);
b=a(1,2);
if colorcoded > 0
   for i=2:b
      xt(1)=x(i-1);
      yt(1)=y(i-1);
      xt(2)=x(i);
      yt(2)=y(i);
      vt=v(i-1);
      if colorcoded == 1
         if vt < 33
            vtcolor=vcolor(dc1,:);
         elseif vt >= 33 && vt < 45
            vtcolor=vcolor(dc2,:);
         elseif vt >= 45 && vt < 54.4
            vtcolor=vcolor(dc3,:);
         elseif vt >= 54.4 && vt < 62.4
            vtcolor=vcolor(dc4,:);
         elseif vt >= 62.4 && vt < 69.5
            vtcolor=vcolor(dc5,:);
         else
            vtcolor=vcolor(dc6,:);
         end
      else   
         vtr=round(0.4267.*(3600.*vt./1852-30));
         vtr=round(vtr*cosize/64);
         %vtr=min(vtr,64);
         vtr=max(vtr,1);
         vtcolor=vcolor(vtr,:);
      end   
      gg=m_plot(xt,yt,'color',vtcolor,'linewidth',thick);
      lc=vtcolor;
   end   
else
   gg=m_plot(x,y,lc,'linewidth',thick);
end
bd=size(xd);
bd=bd(1,2);
if strcmp(dplot,'y')
    if bd > 1   
        ul=zeros(1,bd);
        vl=ul;
       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   
offset=zeros(1,bd)-0.2;
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',dfont,'fontweight','bold','fontsize',textfontsize,'color',textcolor)
end
ixlab=months(month(bd));
if month(bd) > month(1)
   ixlab=[char(months(month(1))) ' - ' char(months(month(bd)))];
elseif month(bd) < month(1)
   ixlab=['December ' ystring ' - January'];
   ystring=ystringp;   
end   
%{
if strcmp(basin,'sh') && month(1) > 6 && month(1) <= month(bd)
   ystring=ystringm;
end
%}
ixlab=char(ixlab);
titstring=[char(stormname) ',  ' ixlab ' ' ystring];
atit=char(stormname);   
 if strcmp(atit, 'NOTNAMED') || strcmp(atit, 'UNNAMED')  
     titstring=['Storm Number ' nstring ',  ' ixlab ' ' ystring];
 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;
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 strcmp(cplot,'y')
   m_plot(xd,yd,'.','color',ccol,'MarkerSize',csize)   
end   
m_gridm('fontname',axisfont,'fontsize',axisfontsize,'fontweight',axisfontweight, ...
    'linestyle',gridline,'backgroundcolor',oceancolor)
title(titstring,'fontsize',axisfontsize+1,'fontweight','bold','interpreter','none')
if colorcoded > 0 && strcmp(cbar,'y')
   if colorcoded == 1
      ucolor=[vcolor(dc1,:);vcolor(dc2,:);vcolor(dc3,:);vcolor(dc4,:);vcolor(dc5,:);vcolor(dc6,:)];
      aa=0:6;
      bb=1/6;
      tstring='Category';
   else
      ucolor=vcolor;
      if strcmpi(units,'knots')
          aa=40:20:175;
          tstring='knots';
      elseif strcmpi(units,'ms')
          aa=20:10:90;
          tstring='m/s';    
      elseif strcmpi(units,'kmh')
          aa=70:35:315;
          tstring='km/h';
      else
          aa=40:20:180;
          tstring='mph';
      end
      bb=1/7;
   end  
   ab=cellstr(num2str(aa'));
   if colorcoded == 1
       ab{7}='';
   end    
   rts=colorbar('horiz');
   colormap(rts,ucolor);
   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)=rr(1)+0.125.*rr(3);
   rr2(1)=posit(1)+0.5*posit(3)-0.5*rr2(3);
   rr2(4)=0.015;
   %
   set(rts,'position',rr2)
   set(rts,'ticklength',0)
   set(rts,'xaxislocation','bottom')
   set(rts,'fontname',axisfont,'fontweight','bold')
   uz=get(rts,'limits');
   set(rts,'ticks',uz(1):(uz(2)-uz(1))*bb:uz(2))
   set(rts,'ticklabels',ab)
   set(gca,'position',posit)
   set(get(rts,'xlabel'),'String',tstring,'fontweight','bold')
end