%
%  This script overlays best tracks on track maps created using alltrmap,
%  which must be run directly before this script is executed.
%  January, 2009.  
%  Copyright 2009 WindRiskTech L.L.C.
%  Last modified August, 2015
%--------------------------------------------------------------------------
%
if exist('shape','var') == 0 || exist('x','var') == 0
    f=msgbox('Need to run alltrmap before running this script');
    return
end   
%
params    %  Load parameters
%
prompt = {'Order by maximum at POI (1), maximum near POI (2), lifetime maximum (3), or random (4)?', ...
           'Number of tracks to plot'};
dlg_title = 'Best 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    
ntracks=str2double(answer{2});
%
bestproc   %  Process best tracks
%
yearbeg=min(bestyears);
yearend=max(bestyears);
rng('shuffle')   %  Set seed of random number generator based on time
%
% Find tracks, allowing for one potential wrap in longitude
%
hold on
nact=min(ntracks,max(size(nsortb)));
for k=1:nact
   if nswitch == 3
   %
   % Ordered by most intense over lifetime
   %
   i=nsortb(k);
   %
   elseif nswitch == 2
   %
   % Ordered by most intense near POI
   %
   i=nsortpb(k);
   %
   elseif nswitch == 1
   %
   % Ordered by most intense at POI
   %
   i=nsortcb(k);
   %
   elseif nswitch == 4
   %
   % Randomly ordered
   %
   i=1+floor((max(nsortb)-1).*rand);
   %
   end 
    np1=0;
    np2=0;
    ntrack=1;
    [~, jmax]=min(vbest(i,:),[],2);
    jmax=jmax-1;
    jmax=max(jmax,1);
    clear x y x2 y2 v v2
    y=zeros(1,jmax);
    x=y;
    v=y;
    x2=x;
    y2=y;
    v2=v;
    for j=1:jmax
          np1=np1+1;
          y(np1)=ay-ay.*(latbest(i,j)+90)./180;
          x(np1)=ax.*(longbest(i,j)+360-displacement)/360;
          v(np1)=vbest(i,j);
          if x(np1) <= ax
             xmax=max(xmax,x(np1));
             xmin=min(xmin,x(np1));
             ymax=max(ymax,y(np1));
             ymin=min(ymin,y(np1));
          end  
          if x(np1) > ax 
              ntrack=2;
              np2=np2+1;
              x2(np2)=x(np1)-ax;
              y2(np2)=y(np1);
              v2(np2)=v(np1);
              np1=np1-1;
              xmax=max(xmax,x2(np2));
              xmin=min(xmin,x2(np2));
              ymax=max(ymax,y2(np2));
              ymin=min(ymin,y2(np2));
          end   
    end    
    if np1 > 1
        plot(x(1:np1),y(1:np1),'color', besttckcolmp,'linewidth',besttckwidth) 
    end    
    if ntrack == 2
        plot(x2,y2,'color', besttckcolmp,'linewidth',besttckwidth)
    end           
%
% Put arrowheads on tracks
%
   if ntrack == 1 && np1 > 1
       ua=x(np1)-x(np1-1);
       va=y(np1)-y(np1-1);
       mag=sqrt(ua^2+va^2);
       ua=bestafac*(ax/3600)*ua/mag;
       va=bestafac*(ay/1800)*va/mag;
       xq=x(np1);
       yq=y(np1);
       g1=quiver(xq,yq,ua,va,10);
       set(g1,'maxheadsize',5);
       set(g1,'showarrowhead','off')
       set(g1,'showarrowhead','on')
       set(g1,'linewidth',besttckwidth,'color', besttckcolmp);
       xb=x(1);
       yb=y(1);
       g2=plot(xb,yb,'.','color', besttckcolmp,'markersize',floor(12.*besttckwidth));
   end    
   if  ntrack == 2 && np2 > 1
       ua=x2(np2)-x2(np2-1);
       va=y2(np2)-y2(np2-1);
       mag=sqrt(ua^2+va^2);
       ua=bestafac*(ax/3600)*ua/mag;
       va=bestafac*(ay/1800)*va/mag;
       xq=x2(np2);
       yq=y2(np2);
       g1=quiver(xq,yq,ua,va,10);
       set(g1,'maxheadsize',5);
       set(g1,'showarrowhead','off')
       set(g1,'showarrowhead','on')
       set(g1,'linewidth',besttckwidth,'color', besttckcolmp);
       xb=x2(1);
       yb=y2(1);
       g2=plot(xb,yb,'.','color', besttckcolmp,'markersize',floor(12.*besttckwidth));      
   end    
end   
hold off
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')
    delete(hButton)
    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
%
% Clean up
%
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end    