%
% This script creates maps of multiple tropical cyclone tracks, ordered 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 randomly. "Prep" must be run before this can be used.
% This verion plots the tracks in Google Earth and requires the
% installation into MATLAB of the free KML toolbox.
%
% Copyright 2012 WindRiskTech L.L.C.
% Last modified August, 2015
%
%--------------------------------------------------------------------------
twidthge=1.5;       %  Width of track
arrowangle=45;      %  Angle of arrowheads
arrowlength=1;      %  Length of arrowheads
startdotsize=1;     %  Size of dots at starting positions
textoffset=0.5;     %  Latitude and Longitude offset of text (degrees)
getextsize=1;       %  Scale factor of text size
%--------------------------------------------------------------------------
try
    kalltrmap=kml('getempfile');
catch err
    errordlg('The free KML toolbox must be installed into MATLAB to use this script')
    return
end    
%
params    %  Load parameters
%
if longmin < 0
    longmin=longmin+360;
end
if longmax < 0
    longmax=longmax+360;
end    
clear xquiv yquiv uquiv vquiv xbeg ybeg
if exist('shearstore','var') == 0
    shape='circ';
    load('temp.mat')
end
fs3='sorted';
load(fs3)
[nsize,jsize]=size(latstore);
%
clear answer
if isempty(storm)
    prompt = {'Order by maximum at POI (1), maximum near POI (2), lifetime maximum (3), or random (4)?', ...
        'Starting track number ', ...
        'Ending track number (hit return to use whole set)', ...
        'Order best tracks by maximum at POI (1), maximum near POI (2), lifetime maximum (3), or random (4)?', ...
        'Number of best tracks to plot (can be zero)'};
    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        
    nt=str2double(answer{2});
    net=str2double(answer{3});
    nbswitch=str2double(answer{4});
    nbtracks=str2double(answer{5});
    if (strcmp(shape,'poly') || city_radius >= 2000) && nbswitch == 1
        nbswitch=2;
    end    
else    
   nswitch=3;
   nbswitch=1;
   nbtracks=0;
end   
if isnan(net)
    net=nsize;
end  
%
[~,jmax]=min(vstore,[],2);
jmax=jmax-1;
%
pifac=acos(-1)/180;
clear col
col=cell(1,2000);
col(1)=cellstr('r');
col(2)=cellstr('y');
col(3)=cellstr('g');
col(4)=cellstr('c');
col(5)=cellstr('m');
%
for i=6:2000
   col(i)=col(i-5);
end   
vcolor=jet;
close
vcolor(1:20,3)=1;
vcolor(1:20,2)=0.6;
if colorscheme == 4
    vcolor=flipud(vcolor);
end    
uv=vcolor(20,:); 
%   
% Determine ordering of tracks
%
rng('shuffle')   %  Set seed of random number generator based on time
nt3=zeros(1,net);
for ntt=nt:net
   if nswitch == 3
   %
   % Ordered by most intense over lifetime
   %
   nt3(ntt)=nsort(ntt);
   %
   elseif nswitch == 2
   %
   % Ordered by most intense near POI
   %
   nt3(ntt)=nsortp(ntt);
   %
   elseif nswitch == 1
   %
   % Ordered by most intense at POI
   %
   nt3(ntt)=nsortc(ntt);
   %
   elseif nswitch == 4
   %
   % Randomly ordered
   %
   nt3(ntt)=1+floor((max(nsort)-1).*rand);
   %
   end
%
end
%
%Label for tracks
%
mmax=0; 
for ntt=net:-1:nt
   if nswitch==1
      nt3s=num2str(ntt);
   else
      for jj=1:max(nsort)
         if nsort(jj)==nt3(ntt)
            nt3s=num2str(jj);
         end
      end   
   end   
   ntcolor=ntt-nt+1;
   clear x y v s x2 y2 v2
   ntrack=1;
   np1=0;
   np2=0;
   xbeg(ntt-nt+1)=longstore(nt3(ntt),1);
   ybeg(ntt-nt+1)=latstore(nt3(ntt),1);
   if xbeg(ntt-nt+1) > 180
       xbeg(ntt-nt+1)=xbeg(ntt-nt+1)-360;
   end    
   %
   % Find tracks
   %
      for i=1:jmax(nt3(ntt))
         np1=np1+1;
         if colorscheme < 4
             v(np1)=vstore(nt3(ntt),i)+sqrt(uinc(nt3(ntt),i).^2+vinc(nt3(ntt),i).^2);
             v(np1)=v(np1).*1852./3600;
         else
             v(np1)=pstore(nt3(ntt),i);
         end     
      end 
   %
   % Color and plot tracks
   %
   if colorscheme == 0
        s=char(col(ntcolor));
        s=rgb(s);
   elseif colorscheme == 1
        s=tcolor;
        s=rgb(s);
   elseif colorscheme == 2 || colorscheme == 3
        s=vcolor(1,:);
   elseif colorscheme == 4
        s=vcolor(64,:);
   end 
   sar=s;
   scolorb=(strcat('FF',rgbconv(fliplr(s))));
   if colorscheme > 1
        xt=zeros(1,2);
        yt=zeros(1,2);
        for i=2:np1
          xt(1)=longstore(nt3(ntt),i-1);
          xt(2)=longstore(nt3(ntt),i);
          yt(1)=latstore(nt3(ntt),i-1);
          yt(2)=latstore(nt3(ntt),i);
          vta=v(i-1);
          if colorscheme == 2
             if vta < 33
                vtcolor=vcolor(1,:);
             elseif vta >= 33 && vta < 45
                vtcolor=vcolor(13,:);
             elseif vta >= 45 && vta < 54.4
                vtcolor=vcolor(30,:);
             elseif vta >= 54.4 && vta < 62.4
                vtcolor=vcolor(39,:);
             elseif vta >= 62.4 && vta < 69.5
                vtcolor=vcolor(49,:);
             else
                vtcolor=vcolor(61,:);
             end
          else   
             if colorscheme == 3 
                vtr=round(vta)-15;
             else
                vtr=round(0.5.*(vta-900));
             end   
             if vtr > 64
                vtr=64;
             end
             if vtr < 1
                vtr=1;
             end   
             vtcolor=vcolor(vtr,:);
          end   
          sar=vtcolor;
          scolor=(strcat('FF',rgbconv(fliplr(vtcolor))));
          kalltrmap.plot(xt,yt,'linecolor',scolor,'linewidth',twidthge);
        end  
       ii=net-ntt+1;
       xquiv(ii)=xt(2);
       yquiv(ii)=yt(2);
       uquiv(ii)=(xt(2)-xt(1));
       vquiv(ii)=(yt(2)-yt(1));    
       di=1./sqrt(uquiv(ii)^2+vquiv(ii)^2)+1e-8;
       ct=uquiv(ii).*di;
       st=vquiv(ii).*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];
       kalltrmap.plot(xa,ya,'linecolor',scolor,'linewidth',twidthge); 
       kalltrmap.plot(xb,yb,'linecolor',scolorb,'linewidth',twidthge); 
       %
       % Add labels to start positions
       %
       if strcmp(textlabel,'y')
           g2=kalltrmap.text(xbeg(ntt-nt+1)+textoffset,ybeg(ntt-nt+1)+textoffset,0,nt3s,'labelscale',getextsize);
       end                  
   %
   else 
       ii=ntt-nt+1;
       xquiv(ii)=longstore(nt3(ntt),np1);
       yquiv(ii)=latstore(nt3(ntt),np1);
       uquiv(ii)=longstore(nt3(ntt),np1)-longstore(nt3(ntt),np1-1);
       vquiv(ii)=latstore(nt3(ntt),np1)-latstore(nt3(ntt),np1-1);    
       kalltrmap.plot(longstore(nt3(ntt),1:np1),latstore(nt3(ntt),1:np1),'linecolor',scolorb,'linewidth',twidthge); 
       scolor=scolorb;
       %
       % Add dots to starting positions
       %
       if strcmp(startdot,'y')
           xd=longstore(nt3(ntt),1);
           if xd > 180
               xd(m)=xd-360;
           end               
           kalltrmap.scatter(xd,latstore(nt3(ntt),1),'iconcolor',scolor,'iconscale',startdotsize);
       end
       %
       % Add labels to start positions
       %
       if strcmp(textlabel,'y')
           xd=longstore(nt3(ntt),1);
           if xd > 180
               xd(m)=xd-360;
           end                          
           g2=kalltrmap.text(xd+textoffset,latstore(nt3(ntt),1)+textoffset,0,nt3s,'labelscale',getextsize);
       end         
       %
       %  Add arrowheads
       %
       di=1./sqrt(uquiv(ii)^2+vquiv(ii)^2)+1e-8;
       ct=uquiv(ii).*di;
       st=vquiv(ii).*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=[longstore(nt3(ntt),np1) longstore(nt3(ntt),np1)+dx1];
       ya=[latstore(nt3(ntt),np1) latstore(nt3(ntt),np1)+dy1];
       xb=[longstore(nt3(ntt),np1) longstore(nt3(ntt),np1)+dx2];
       yb=[latstore(nt3(ntt),np1) latstore(nt3(ntt),np1)+dy2];
       kalltrmap.plot(xa,ya,'linecolor',scolorb,'linewidth',twidthge); 
       kalltrmap.plot(xb,yb,'linecolor',scolorb,'linewidth',twidthge); 
   end
end
if colorscheme > 0 && strcmp(startdot,'y')
    kalltrmap.scatter(xbeg,ybeg,'iconcolor',scolor,'iconscale',startdotsize);
end
%
%  Now plot best tracks
%
bestproc   %  Process best tracks
%
nmax=max(size(nsortb));
%
rng('shuffle')   %  Set seed of random number generator based on time
if nbtracks > 0
    nact=min(nbtracks,nmax);
    for k=1:nact
       if nbswitch == 3
       %
       % Ordered by most intense over lifetime
       %
       i=nsortb(k);
       %
       elseif nbswitch == 2
       %
       % Ordered by most intense near POI
       %
       i=nsortpb(k);
       %
       elseif nbswitch == 1
       %
       % Ordered by most intense at POI
       %
       i=nsortcb(k);
       %
       elseif nbswitch == 4
       %
       % Randomly ordered
       %
       i=1+floor((max(nsortb)-1)*rand);
       %
       end
        [~, jmax]=min(vbest(i,:),[],2);
        jmax=jmax-1;
        sb=rgb(besttckcolor);
        sbb=strcat('FF',fliplr(rgbconv(sb)));
        kalltrmap.plot(longbest(i,1:jmax),latbest(i,1:jmax),'linecolor',sbb,'linewidth',besttckwidth); 
       %
       % Add dots to starting positions
       %
       xd=longbest(i,1);
       if xd > 180
           xd=xd-360;
       end               
       kalltrmap.scatter(xd,latbest(i,1),'iconcolor',sbb,'iconscale',startdotsize);
       %
       %  Add arrowheads
       %
       xquiv=longbest(i,jmax);
       yquiv=latbest(i,jmax);
       uquiv=longbest(i,jmax)-longbest(i,jmax-1);
       vquiv=latbest(i,jmax)-latbest(i,jmax-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=[longbest(i,jmax) longbest(i,jmax)+dx1];
       ya=[latbest(i,jmax) latbest(i,jmax)+dy1];
       xb=[longbest(i,jmax) longbest(i,jmax)+dx2];
       yb=[latbest(i,jmax) latbest(i,jmax)+dy2];
       kalltrmap.plot(xa,ya,'linecolor',sbb,'linewidth',besttckwidth); 
       kalltrmap.plot(xb,yb,'linecolor',sbb,'linewidth',besttckwidth);  
    end   
end
kalltrmap.run
%
% Clean up
%
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 