%
% This script contours the number density of all 6 hour best track points,  
% filtered (if desired) to include only those events whose winds speeds 
% exceed some threshold, peakv. It requires installation of the
% "m_map" routines. Tracks interpolated to 2 hours. 
% Copyright 2009 WindRiskTech L.L.C.
% Last modified April, 2024 to take advantage of "ismembertol"
% Modified April, 2026 to divide z2 by mresbest instead of mresbest^2
%------------------------------------------------------------------------
%
if exist('shape','var') == 0 || exist('xmax','var') == 0
    f=msgbox('Need to run trackdensity before running this script');
    return
end
%
params    %  Load parameters
%               
figure
x=xmin:mresbest:xmax;
y=ymin:mresbest:ymax;
%
pi=acos(-1);
pifac=pi./180;
%
bestproc   %  Process best tracks
%
yearbeg=min(bestyears);
yearend=max(bestyears);
%
if strcmp(bas,'GB')
    projection=gproject;
end 
%
[n,m]=size(vbest);
%
yearsize=max(size(bestyears));
%
lat=latbest;long=longbest;
lat(vbest < minv)=NaN;  % Negate portions of tracks where vnet < minv
long(vbest < minv)=NaN;
%
mask = (vmaxb' < peakv);  % logical vector (nCurves x 1)

lat(mask, :)  = NaN; % Throw away whole track if vmax < peakv
long(mask, :) = NaN;
%
z2=countCurvesInLatLonBoxes_interp(lat, long, x, y, 0.5, true);
z2=z2./(yearsize*mresbest);
x=x+0.5*(x(2)-x(1));
x(end)=[];
y=y+0.5*(y(2)-y(1));
y(end)=[];
an=length(y);
am=length(x);
%
clf('reset')
m_proj(projection,'long',[xmin xmax],'lat',[ymin ymax]);
m_gridm('fontname',axisfont,'fontsize',axisfontsize,'fontweight',axisfontweight, ...
    'linestyle',gridline,'backgroundcolor',oceancolor)
%
if strcmp(bas,'MT')
     m_gshhs_l('patch',landcolor,'edgecolor','none');
else     
     m_coast('patch',landcolor,'edgecolor','none');
end   
if strcmp(pstates,'y')
    m_states(xmin, xmax, ymin, ymax,'color',stcolor)
end  
hold on
cmax=(max(max(z2)));
cint=cmax/15;
if verLessThan('matlab', '8.5.0')  %  Apply this to pre 2015 (R2014b) matlab versions
    [C,h]=m_contourf(x,y,z2,cint:cint:cmax);
    thismap = get(gcf, 'Colormap');
    thismap(1,:)=[1 1 1];
    colormap(thismap)
    clear alpha
    alpha(0.99)
    hPatch = findobj(h,'Type','patch');
    set(hPatch,'facealpha',wtrans); 
    set(findobj(h,'FaceColor',[1 1 1]), 'FaceAlpha',0)
else
    [C1,h1]=m_contour(x,y,z2,cint:cint:cmax);
    if strcmp(wfill,'y')
        %z(z<0.08*max(max(z)))=NaN;
        h=m_pcolor(x,y,z2);
        a=colormap(wmap);
        for i=1:3
            a(i,:)=oceancolor;
        end
        colormap(a)
        set(h,'EdgeAlpha',0,'FaceAlpha',wtrans);
        set(h,'FaceColor','interp');
    end
end    
%
jj=colorbar;
set(jj,'fontname',axisfont,'fontweight','bold')
cbIm = findobj(jj,'Type','image');
alpha(cbIm,wtrans)
yvec=get(jj,'ylim');
yvec(1,1)=0.1*yvec(1,2);
set(jj,'ylim',yvec)
hold off
%
titbest=['Best track density from ' num2str(yearbeg) ' to ' num2str(yearend)];
title(titbest,'fontsize',12,'fontweight','bold','interpreter','none')
%
% 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 