%
% Plots the tracks and intensities of events from best track data sets.
% The latter must be available in matlab binary format. 
% Copyright 2009 WindRiskTech L.L.C.
% Modified October, 2020 to account for MATLAB's change in colormaps
% Last modified September, 2021
%------------------------------------------------------------------------
%
params    %  Load parameters
%
if longmin < 0
   longmin=longmin+360;
end
if longmax < 0
   longmax=longmax+360;
end
%
if strcmpi(swap,'y')
   dummy=landcolor;
   landcolor=oceancolor;
   oceancolor=dummy;
end   
lc=tcolor;
ld=landcolor;
%
% Specify storm basin (at, ep, wp, io, sh, or bm)
%
	basin=input('Input basin (al, ep, wp, io, or sh)','s');
%
basin2=basin;
if strcmp(basin,'al')
   basin2='at';
end   
if strcmp(basin,'bm')
    basin2='bom';
end    
file=strcat(trkdir,'/',basin2,'tracks');
clf
clear x* y* v* u* date offset
xmin=longmin;
xmax=longmax;
ymin=latmin;
ymax=latmax;
load(file)
%
% Convert wind speed units and define S-S category boundaries
%
vsm=vsm*(3600/1852)*wunitsfac;
dc(1)=1;dc(2)=64;dc(3)=83;dc(4)=96;dc(5)=113;dc(6)=137;
dc=dc*wunitsfac;
vcmax=10*round(0.1*wunitsfac*190);
vcmin=10*round(0.1*wunitsfac*20);
%
[iq,mq]=size(daym);
year=input('Year of storm');
clear temp
temp=input('Storm number or name','s');
clear ns
ns=str2num(temp); %#ok<ST2NM>
[nname,mname]=size(stname);
if mname > nname
    stname=stname';
end    
if isempty(ns)
    for i=1:iq
        if strcmpi(temp,strtrim(stname(i,:))) && yeari(i) == year
            ns=nsi(i);
        end
    end
    if isempty(ns)
        disp('Storm not found, quitting')
        return
    end
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('jet(264)')
colormap(trackcolormap)
vcolor=colormap;
[cosize,~]=size(vcolor);
dcolor=1+(cosize-1)*(dc-dc(1))/(dc(end)-dc(1));
dcolor=round(dcolor);
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)
      inew=i;
      clear stormname
      if min(size(stname)) == 1
          stormname=stname(i);
      else    
        stormname=strtrim(stname(i,:)); 
      end  
      x=zeros(1,mq);
      y=x;
      v=x;
      xd=x;
      yd=x;
      day=x;
      month=x;
      date=cell(mq);
      for j=1:mq
         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);
               day(m)=daym(i,j);
               month(m)=monthm(i,j);
               daystring=num2str(day(m));
               date(m)=cellstr(daystring);
            end      
         end
      end
      x(n+1:end)=[];
      y(n+1:end)=[];
      v(n+1:end)=[];
      xd(m+1:end)=[];
      yd(m+1:end)=[];
      day(m+1:end)=[];
      month(m+1:end)=[];
      date(m+1:end)=[];
   end
end
xstorm=x;ystorm=y;
if iswitch == 0
   message='No data found for this storm, quitting';
   disp(message)
	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]);
m_gridm('fontname',axisfont,'fontsize',axisfontsize,'fontweight',axisfontweight, ...
   'linestyle',gridline,'backcolor',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   
%
if (xmax-xmin)*(ymax-ymin) > 1e3 % (Usually not enough memory for higher resolution coasts)
    cres=1;
end
%
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);
[ncolors,~]=size(vcolor);
hold on
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);
      vta=v(i-1);
      if colorcoded == 1
         if vta < dc(2)
            vtcolor=vcolor(dcolor(1),:);
         elseif vta >= dc(2) && vta < dc(3)
            vtcolor=vcolor(dcolor(2),:);
         elseif vta >= dc(3) && vta < dc(4)
            vtcolor=vcolor(dcolor(3),:);
         elseif vta >= dc(4) && vta < dc(5)
            vtcolor=vcolor(dcolor(4),:);
         elseif vta >= dc(5) && vta < dc(6)
            vtcolor=vcolor(dcolor(5),:);
         else
            vtcolor=vcolor(dcolor(6),:);
         end
      else 
         vtr=round(cosize*(vta-vcmin)/(vcmax-vcmin));
         %vtr=round(0.8583*(vta-30*1852/3600));
         %vtr=round(vtr*cosize/64);
         vtr=max(vtr,1);
         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 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
   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',gfont,'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 ' 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 max(size(atit))~= 8
     if isempty(atit) || (strcmp(atit(1),'N') && strcmp(atit(2),'O'))
         titstring=['Storm number ' nstring ',  ' ixlab ' ' ystring];
     end    
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   
title(titstring,'fontsize',axisfontsize,'fontweight','bold','interpreter','none')
if colorcoded > 0 && strcmp(cbar,'y')
   if colorcoded == 1    
      ucolor=[vcolor(dcolor(1),:);vcolor(dcolor(2),:);vcolor(dcolor(3),:); ...
          vcolor(dcolor(4),:);vcolor(dcolor(5),:);vcolor(dcolor(6),:)];
      %aa=0:6;
      aa={'TS','Cat 1','Cat 2','Cat3','Cat4','Cat 5',' '};
      bb=1/6;
      tstring='Category';      
   else
      ucolor=vcolor;
      if strcmpi(wunits,'kts')
          aa=40:20:180;
          tstring='knots';
      elseif strcmpi(wunits,'kmh')
          aa=70:35:315;
          tstring='km/h';
      else
          aa=40:20:180;
          tstring='mph';
      end
      bb=1/7;  
   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)
   set(rts,'xaxislocation','bottom')
   set(rts,'fontname',axisfont,'fontweight','bold')
   set(rts,'xtick',0:bb:1,'xticklabel',aa)
   set(gca,'position',posit)
   set(get(rts,'xlabel'),'String',tstring,'fontweight','bold')
   
end
%
% Add button to access documentation of script
%
if strcmp(Docdisp,'y')
    tempname=mfilename;
    mname = strcat('graphics_docs/html/',tempname,'_doc.html');
    comstr=['open ',mname,';delete(hButton)'];
    labelStr = '<html><center><a href="">Doc';
    hButton = uicontrol('string',labelStr,'pos',[10,10,50,20],'callback',comstr); 
end