%
% 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.
% November 2008
% Copyright 2009 WindRiskTech L.L.C.
% Last Modified July, 2020 to account for larger size of colormaps
%-------------------------------------------------------------------------
%
params    %  Load parameters
%
if longmin < 0
  longmin=longmin+360;
end
if longmax < 0
  longmax=longmax+360;
end   
%clear daymod v vtrans vpot shear u850store v850store yearstore
if strcmp(swap,'y')
   dummy=landcolor;
   landcolor=oceancolor;
   oceancolor=dummy;
end   
lc=tcolor;
ld=landcolor;
colormap('default')
% 
if exist('shearstore','var') == 0
    shape='circ';
    load('temp.mat')
else
    load('temp.mat','vt')
    load('temp.mat','ut')
    warning ('off','all');
    load('temp.mat','yearstore')
end    
warning ('on','all');
if exist('nsortc','var') == 0
    fs3='sorted';
    load(fs3)
end    
if strcmp(bas,'GB')
    projection=gproject;
end    
%
if isempty(storm)
    prompt = {'Order by maximum at POI (1), maximum near POI (2), lifetime maximum (3), or by track number (4)?', ...
        'Track number'};
    dlg_title = 'Track Parameters';
    num_lines = 1;
    answer = inputdlg(prompt,dlg_title,num_lines);
    if isempty(answer) || isempty(answer{1})
        return
    end    
    nswitch=str2double(answer{1});
    if (strcmp(shape,'poly') || city_radius >= 2000) && nswitch == 1
        warndlg('Warning: Winds at POI not meaningful for this set; ordering by maximum near POI')
        nswitch=2;
    end        
    nttemp=str2double(answer{2});
else    
   nswitch=3;
   prompt={'Track Number'};
   dlg_title='Track Specification';
   answer = inputdlg(prompt,dlg_title,num_lines);
   if isempty(answer) || isempty(answer{1})
       return
   end
   nttemp=str2double(answer{1});
end   
%
if nswitch == 1
   nt=nsortc(nttemp);
elseif nswitch == 2
   nt=nsortp(nttemp);
elseif nswitch == 3   
   nt=nsort(nttemp);
else
   nt=nttemp;
end   
%
clf
pifac=acos(-1)/180;
ntstring=num2str(nt);
xmin=longmin;
xmax=longmax;
ymin=latmin;
ymax=latmax;
[ndum,m8max]=size(vstore);
[~,jmax]=min(vstore,[],2);
jmax(vstore(:,end)~=0)=m8max;
jmax=jmax-1;
iq=jmax(nt);
iqo=iq;
%
clear vtrans shear vpot monthm daym hour pres x y
vtrans=wunitsfac*(sqrt(ut(nt,1:iq).^2+vt(nt,1:iq).^2));
v=zeros(1,iq);   
v(1:iq)=wunitsfac*vnet(nt,1:iq);
pres(1:iq)=pstore(nt,1:iq);
x(1:iq)=longstore(nt,1:iq);
y(1:iq)=latstore(nt,1:iq);
for i=1:iq
    if x(1) > 200 && x(i) < 50
        x(i)=x(i)+360;
    end
end  
vpot(1:iq)=wunitsfac*vpstore(nt,1:iq);
shear(1:iq)=wunitsfac*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 monthm(i) == 12 && daym(i) > 31
      monthm(i)=1;
      daym(i)=daym(i)-31;
   end   
end
amod=mod(hour(1),2);
if amod ~= 0
    hour=hour-amod;
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;
[qcolor,~]=size(vcolor);
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 
daytime=daym+hour./24;     
daymod=daytime;
delday=max((daytime(2)-daytime(1)),(daytime(3)-daytime(2)));
iswitch=0;
daymax=floor(daytime(1));
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  
n=0;
m=0;
iswitch=0;
for i=1:iq
    if hour(i) <= 0.5
        m=m+1;
    end
end
if m == 0
    msgbox('This track too short to display','Error','error')
    return
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) <= 1.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
bd=max(size(xd));
ixlab=months(month(bd));
if month(bd) > month(1)
      ixlab=[char(months(month(1))) char(' - ') char(months(month(bd)))];
elseif month(bd) < month(1)
      ixlab='December - January';
end   
if exist('yearstore','var') == 1 
    ixlab=[char(ixlab), ' ', char(num2str(yearstore(nt)))];
end    
titstring=['track number ' ntstring ',   ' ixlab ];
titstring2=['track number ' ntstring  ];
if isempty(storm)
    datea(1)=upper(datea(1));
    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;char(ixlab)};
save templot nt ntstring tit3 date month day xd yd vcolor months ...
    vpot shear vtrans v x y monthm daym hour jmax iq daymod ...
    tit1 pres daymax ixlab lc ld bd
%
nplot=2;
while nplot > 1
    nplot=menu('Track plots','Stop','Track on map background', 'Track on satellite background', ...
        'Track in Google Earth','Max winds along track','Vmax including secondary eyewalls', ...
        'Wind radii','Minimum pressure along track', ...
        'Shear along track','Maximum wind swath','Surface windfield at selected time', ...
        'Surface wind arrows at selected time',  ...
        'Surface pressure field',  ...
        'Wind time series at selected point','Accumulated rain swath', ...
        'Rain field at selected time','Time series of rain at selected point', ...
        'Print as .png figure');
    if nplot == 1
        %  Delete temporary file
        delete templot.mat
        %
        break
    elseif nplot == 2   
        load templot
        params
        clf('reset')
        colormap('default')
        colormap(trackcolormap)
        vcolor=colormap;
        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   
        
        if strcmp(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]);
        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);
           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);
              %vta=(1852/3600)*v(i-1);
              vta=(1852/3600)*v(i);
              if colorcoded == 1
                 if vta < 33
                    vtcolor=vcolor(1,:);
                 elseif vta >= 33 && vta < 45
                    vtcolor=vcolor(round(0.2*qcolor),:);
                 elseif vta >= 45 && vta < 54.4
                    vtcolor=vcolor(round(0.47*qcolor),:);
                 elseif vta >= 54.4 && vta < 62.4
                    vtcolor=vcolor(round(0.6*qcolor),:);
                 elseif vta >= 62.4 && vta < 69.5
                    vtcolor=vcolor(round(.77*qcolor),:);
                 else
                    vtcolor=vcolor(round(0.95*qcolor),:);
                 end
              else   
                 vtr=round(0.4267.*(qcolor/64)*(3600.*vta./1852-30));
                 vtr=min(vtr,qcolor);
                 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
        if strcmp(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
        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,'bac',oceancolor)
        title(tit3,'fontsize',11,'fontweight','bold','interpreter','none')
        if colorcoded > 0 && strcmp(cbar,'y')
           if colorcoded == 1
              ucolor=[vcolor(1,:);vcolor(round(0.2*qcolor),:);vcolor(round(0.47*qcolor),:); ...
                  vcolor(round(0.6*qcolor),:);vcolor(round(0.77*qcolor),:);vcolor(round(0.95*qcolor),:)];
              aa={'TS';'1';'2';'3';'4';'5';' '};
              bb=1/6;
              tstring='Category';
           else
              ucolor=vcolor;
              aa=40:10:180;
              bb=1/14;
              tstring=wunitslabel;
           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')
            try
                delete(hButton)
            catch
            end    
            tempname='trplot';
            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             
    %
    elseif nplot == 3  
        load templot
        clf('reset')
        %clear x y
        params
        %
        %  Read in image
        %
        imagename=strcat('map_',imageres);
        try
            load(imagename);
            %A=imread('eocene_map.jpg');
        catch err
            errordlg('Cannot find image file; may neeed to run initialize.m')
            break
        end    
        %
        [ay, ax, dum]=size(A);
        aspect=ax/ay;
        %       
        image(A)
        %
        hold on
        x(1:iq)=ax.*(longstore(nt,1:iq)+360-displacement)/360;
        y(1:iq)=ay-ay.*(latstore(nt,1:iq)+90)./180;
        %
        if iq > 1
        %
        ntrack=1;
        mini=iq;
        np2=0;
        x2=zeros(1,iq);
        y2=x2;
        v2=x2;pres2=x2;vpot2=x2;shear2=x2;monthm2=x2;daym2=x2;hour2=x2;
        vtrans2=x2;
        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);
                v2(np2)=v(i);
                pres2(np2)=pres(i);
                vpot2(np2)=vpot(i);
                shear2(np2)=shear(i);
                monthm2(np2)=monthm(i);
                daym2(np2)=daym(i);
                hour2(np2)=hour(i);
                vtrans2(np2)=vtrans(i);
            end
        end
        x2(np2+1:end)=[];
        y2(np2+1: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);
              vta=(1852/3600)*v(i-1);
              if colorcoded == 1
                 if vta < 33
                    vtcolor=vcolor(1,:);
                 elseif vta >= 33 && vta < 45
                    vtcolor=vcolor(round(0.2*qcolor),:);
                 elseif vta >= 45 && vta < 54.4
                    vtcolor=vcolor(round(0.47*qcolor),:);
                 elseif vta >= 54.4 && vta < 62.4
                    vtcolor=vcolor(round(0.6*qcolor),:);
                 elseif vta >= 62.4 && vta < 69.5
                    vtcolor=vcolor(round(.77*qcolor),:);
                 else
                    vtcolor=vcolor(round(0.95*qcolor),:);
                 end
              else   
                 vtr=round(0.4267.*(qcolor/64)*(3600.*vta./1852-30));
                 %vtr=round(vta)-15;
                 %vtr=round(vtr*(qcolor/64));
                 vtr=min(vtr,qcolor);
                 vtr=max(vtr,1);
                 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 strcmp(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).*textdfacsat.*((vl./vabs)+offset);
        dy=15.*(ay/1600).*textdfacsat.*(-ul./vabs);
        xdb=xd+dx;
        ydb=yd+dy;
        g=text(xdb,ydb,date);
        set(g,'fontname',gfont,'fontweight',gfontweight,'fontsize',textfsizesat,'color',textcolorsat)
        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 strcmp(cplot,'y')
           plot(xd,yd,'.','color',ccol,'MarkerSize',csize)   
        end   
        %
        end
        %
        if ntrack == 2
            iq=np2;
            clear x y
            x=x2;
            y=y2;
            v=v2;
            pres=pres2;
            vpot=vpot2;
            shear=shear2;
            monthm=monthm2;
            daym=daym2;
            hour=hour2;
            vtrans=vtrans2;
        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 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);
           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);
              vta=(1852/3600)*v(i-1);
              if colorcoded == 1
                 if vta < 33
                    vtcolor=vcolor(1,:);
                 elseif vta >= 33 && vta < 45
                    vtcolor=vcolor(round(0.2*qcolor),:);
                 elseif vta >= 45 && vta < 54.4
                    vtcolor=vcolor(round(0.47*qcolor),:);
                 elseif vta >= 54.4 && vta < 62.4
                    vtcolor=vcolor(round(0.6*qcolor),:);
                 elseif vta >= 62.4 && vta < 69.5
                    vtcolor=vcolor(round(.77*qcolor),:);
                 else
                    vtcolor=vcolor(round(0.95*qcolor),:);
                 end
              else   
                 vtr=round(0.4267.*(qcolor/64)*(3600.*vta./1852-30)); 
                 %vtr=round(vta)-15;
                 %vtr=round(vtr*(qcolor/64));
                 vtr=min(vtr,qcolor);
                 vtr=max(vtr,1);                 
                 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).*textdfacsat.*((vl./vabs)+offset);
        dy=15.*(ay./1600).*textdfacsat.*(-ul./vabs);
        xdb=xd+dx;
        ydb=yd+dy;
        g=text(xdb,ydb,date);
        set(g,'fontname',gfont,'fontweight',gfontweight,'fontsize',textfsizesat,'color',textcolorsat)
        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(round(0.2*qcolor),:);vcolor(round(0.47*qcolor),:); ...
                  vcolor(round(0.6*qcolor),:);vcolor(round(0.77*qcolor),:);vcolor(round(0.95*qcolor),:)];
              aa={'TS';'1';'2';'3';'4';'5'};
              bb=1/6;
              tstring='Category';
           else
              ucolor=vcolor;
              aa=40:10:180;
              bb=1/14;
              tstring=wunitslabel;
           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,'fontname',gfont,'fontweight','bold')
           set(rts,'xtick',0:bb:1,'xticklabel',aa)
           set(gca,'position',posit)
           set(get(rts,'xlabel'),'String',tstring,'fontweight','bold')
        end
        %
        end
        title(tit3,'fontsize',11,'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,'fontname',gfont,'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',gridlstyle,'linewidth',glinewidth,'color',gridcolor)
        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',gridlstyle,'linewidth',glinewidth,'color',gridcolor)
        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);
        xmax=min(xmax,ax);
        set(gca,'xlim',[xmin xmax],'ylim', [ymin ymax]) 
        %
        % Add button to access documentation of script
        %
        if strcmp(Docdisp,'y')
            try
                delete(hButton)
            catch
            end    
            tempname='trmap';
            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                     
    elseif nplot == 4
        %------------------------------------------------------------------
        twidthge=4;         %  Width of track
        arrowangle=45;      %  Angle of arrowheads
        arrowlength=1;      %  Length of arrowheads
        dotsize=1;          %  Size of dots at date positions
        dotcolor=[0 1 1];   %  Color of date dots
        textoffset=1;       %  Latitude and Longitude offset of text (degrees)
        getextsize=2;       %  Scale factor of text size
        %------------------------------------------------------------------       
        try
            ktrplot=kml('getempfile');
        catch err
            errordlg('The free KML toolbox must be installed into MATLAB to use this script')
            return
        end    
        load templot
        params
        %close
        vcolor=jet;
        n=0;
        m=0;
        iswitch=0;
        for i=1:iqo
            if hour(i) <= 0.5
                m=m+1;
            end
        end
        xd=zeros(1,m);
        yd=zeros(1,m);
        day=zeros(1,m);
        date=cellstr(num2str(daym(1:m)));
        month=date;
        m=0;
        for i=1:iq
            if hour(i) <= 0.5
               m=m+1;
               if x(i) <= 180
                   xd(m)=x(i);
               else
                   xd(m)=x(i)-360;
               end    
               yd(m)=y(i);
               day(m)=daym(i);
               monthstring=num2str(monthm(i));
               month(m)=cellstr(monthstring);
               daystring=num2str(day(m));
               date(m)=cellstr(daystring);
            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);
           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);
              vta=(1852/3600)*v(i-1);
              if colorcoded == 1
                if vta < 33
                    vtcolor=vcolor(1,:);
                 elseif vta >= 33 && vta < 45
                    vtcolor=vcolor(round(0.2*qcolor),:);
                 elseif vta >= 45 && vta < 54.4
                    vtcolor=vcolor(round(0.47*qcolor),:);
                 elseif vta >= 54.4 && vta < 62.4
                    vtcolor=vcolor(round(0.6*qcolor),:);
                 elseif vta >= 62.4 && vta < 69.5
                    vtcolor=vcolor(round(.77*qcolor),:);
                 else
                    vtcolor=vcolor(round(0.95*qcolor),:);
                 end
              else   
                 vtr=round(vta)-15;
                 vtr=round(vtr*(qcolor/64));
                 vtr=min(vtr,qcolor);
                 vtr=max(vtr,1);
                 vtcolor=vcolor(vtr,:);
              end   
              s=(strcat('FF',rgbconv(fliplr(vtcolor))));
              ktrplot.plot(xt,yt,'linecolor',s,'linewidth',twidthge);
           end                  
        else
           ss=rgb(tcolor);
           s=(strcat('FF',rgbconv(fliplr(ss)))); 
           ktrplot.plot(x,y,'linecolor',s,'linewidth',twidthge);
           xt(1)=x(end-1);
           xt(2)=x(end);
           yt(1)=y(end-1);
           yt(2)=y(end);
        end
        if strcmp(cplot,'y')
           sd=(strcat('FF',rgbconv(fliplr(ccol))));
           ktrplot.scatter(xd,yd,'iconcolor',sd,'iconscale',dotsize);
        end
        if strcmp(dplot,'y')  
            qd=max(size(xd));
            for kk=1:qd
                ud=strcat(month(kk),'/',date(kk));
                ktrplot.text(xd(kk)+textoffset,yd(kk)+textoffset,0,ud,'labelscale',getextsize);
            end
        end
       xquiv=xt(2);
       yquiv=yt(2);
       uquiv=(xt(2)-xt(1));
       vquiv=(yt(2)-yt(1));    
       di=1./sqrt(uquiv^2+vquiv^2)+1e-8;
       ct=uquiv.*di;
       st=vquiv.*di;
       angle=pifac.*(180-arrowangle);
       san=sin(angle);
       can=cos(angle);
       dx1=arrowlength*(ct*can+st*san);
       dx2=arrowlength*(ct*can-st*san);
       dy1=arrowlength*(st*can-ct*san);
       dy2=arrowlength*(st*can+ct*san);
       xa=[xt(2) xt(2)+dx1];
       ya=[yt(2) yt(2)+dy1];
       xb=[xt(2) xt(2)+dx2];
       yb=[yt(2) yt(2)+dy2];
       ktrplot.plot(xa,ya,'linecolor',s,'linewidth',twidthge); 
       ktrplot.plot(xb,yb,'linecolor',s,'linewidth',twidthge);  
       ktrplot.run
    elseif nplot == 5
        load templot
        params
        clf 
        shear=smooth(shear);
        h=plot(daymod,wunitsfac*v,daymod,wunitsfac*vtrans,daymod,wunitsfac*vpot,daymod,wunitsfac*shear);
        set(h,'linewidth',3)
        set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        legend('V_{max}','V_{translation}','V_{pot}','Shear','location','best')
        if monthm(1) ~= monthm (iq)
          jg=get(gca,'xticklabel');
          hj=str2num(char(jg)); %#ok<ST2NM>
          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,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        end   
        xlabel(ixlab,'fontweight','bold','fontsize',12)
        ylabl=char(strcat('Maximum surface wind speed (',wunitslabel,')'));
        ylabel(ylabl,'fontweight','bold','fontsize',12)
        title(tit1,'fontsize',12,'fontweight','bold','interpreter','none')
        %
        % Add button to access documentation of script
        %
        if strcmp(Docdisp,'y')
            try
                delete(hButton)
            catch
            end    
            tempname='trspeeds';
            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                     
    %
    elseif nplot == 6
        vdoubletime
    elseif nplot == 7
        rdoubletime
    elseif nplot == 8
    %
        load templot
        params
        clf
        cla
        h=plot(daymod,pres);
        set(h,'linewidth',3)
        set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        if monthm(1) ~= monthm (iq)
          jg=get(gca,'xticklabel');
          hj=str2num(char(jg)); %#ok<ST2NM>
          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,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        end   
        xlabel(ixlab,'fontweight','bold','fontsize',12)
        ylabel('Central pressure (hPa)','fontweight','bold','fontsize',12)
        title(tit1,'fontsize',12,'fontweight','bold','interpreter','none')
        %
        % Add button to access documentation of script
        %
        if strcmp(Docdisp,'y')
            try
                delete(hButton)
            catch
            end    
            tempname='trpress';
            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                             
    %
    elseif nplot == 9
    %
        load templot   
        params
        clf
        h=plot(daymod,wunitsfac*shear);
        set(h,'linewidth',3)
        set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        if monthm(1) ~= monthm (iq)
          jg=get(gca,'xticklabel');
          hj=str2num(char(jg)); %#ok<ST2NM>
          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,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        end   
        xlabel(ixlab,'fontweight','bold','fontsize',12)
        ylabl=char(strcat('850 - 250 hPa shear (',wunitslabel,')'));
        ylabel(ylabl,'fontweight','bold','fontsize',12)
        title(tit1,'fontsize',12,'fontweight','bold','interpreter','none')
        %
        % Add button to access documentation of script
        %
        if strcmp(Docdisp,'y')
            try
                delete(hButton)
            catch
            end    
            tempname='trshear';
            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                                     
    elseif nplot == 10
        load templot
        swath
    elseif nplot == 11
        load templot
        windfield
    elseif nplot == 12
        load templot
        windfield_quiver    
    elseif nplot == 13
        load templot
        pressurefield    
    elseif nplot == 14
        load templot
        pointplot
    elseif nplot == 15
        load templot
        rainswath
    elseif nplot == 16
        load templot
        rainfield
    elseif nplot == 17
        load templot
        rainseries
    elseif nplot == 18
        prompt = {'Enter Figure Name'};
        dlg_title = 'Figure name';
        num_lines = 1;
        answer = inputdlg(prompt,dlg_title,num_lines);
        if isempty(answer) || isempty(answer{1})
            return
        end    
        figurename=answer{1};
        print(figurename,'-dpng','-r300')
    end
end    
%
% Clean up
%
if exist('Clearvlbs','var') && strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 