%
% This script contours the number density of all 2 hour 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.
% Copyright 2009 WindRiskTech L.L.C.
% Last modified August, 2015
% Modified April, 2026 to divide z by mres instead of mres^2
%------------------------------------------------------------------------
%
params    %  Load parameters
%
clf('reset')
clear latscat longscat y x z vmax2
if exist('shearstore','var') == 0
    shape='circ';
    load('temp.mat','vnet','latstore','longstore','yearstore','freqyear', ...
        'freq','vmax','bas')
end
load sorted
[nn,m]=size(vnet);
if strcmp(bas,'GB')
    projection=gproject;
end    
if exist('yearstore','var') == 0
    yearstore=2000+zeros(1,nn);
    freqyear=freq;
end
qy=length(freqyear);
%
pi=acos(-1);
pifac=pi/180;
%
if strcmp(mapmode,'auto')
    xmin=min(nonzeros(longstore))-dellong;
    xmax=max(nonzeros(longstore))+dellong;
    ymin=min(nonzeros(latstore))-dellat;
    ymax=max(nonzeros(latstore))+dellat;
    xmin=mres*floor(xmin/mres);
    xmax=mres*ceil(xmax/mres);
    ymin=mres*floor(ymin/mres);
    ymax=mres*ceil(ymax/mres);
else
    xmin=longmin;
    xmax=longmax;
    ymin=latmin;
    ymax=latmax;
end 
if ((strcmp(bas,'AL') || strcmp(bas,'AP')) && max(latstore(:,1)) < 55) || strcmp(bas,'EW')
    xmin=max(xmin,180);
    xmax=min(xmax,360);
    ymin=max(ymin,0);
    %ymax=80;
%elseif strcmp(bas,'GB')
%    xmin=0;
%    xmax=360;
%    ymin=-84;
%    ymax=84;
end    
x=xmin:mres:xmax;
y=ymin:mres:ymax;
%
lat=latstore;long=longstore;
lat(vnet < minv)=NaN;  % Negate portions of tracks where vnet < minv
long(vnet < minv)=NaN;
%
mask = (vmax' < peakv);  % logical vector (nCurves x 1)

lat(mask, :)  = NaN; % Throw away whole track if vmax < peakv
long(mask, :) = NaN;
%
z=countCurvesInLatLonBoxes_interp(lat, long, x, y, 0.5, true);
z=freq*z/(nn*mres);
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);
%
% Calculate new map limits
%
imin=1;
imax=an;
kmin=1;
kmax=am;
if strcmp(mapmode,'auto')
    xmin1=xmax;
    xmax1=xmin;
    ymin1=ymax;
    ymax1=ymin;
    kmin=am;
    kmax=1;
    imin=an;
    imax=1;
    zcrit=0.01*max(max(z));
    for j=1:an
        for k=1:am
            if z(j,k) > zcrit
                xmin1=min(xmin1,x(k));
                xmax1=max(xmax1,x(k));
                ymin1=min(ymin1,y(j));
                ymax1=max(ymax1,y(j));
                kmin=min(kmin,k);
                kmax=max(kmax,k);
                imin=min(imin,j);
                imax=max(imax,j);
            end    
        end
    end    
    xmin=xmin1;
    ymin=ymin1;
    xmax=xmax1;
    ymax=ymax1;
end
%if strcmp(bas,'AL') || strcmp(bas,'EW')
%    xmin=max(xmin,220);
%    xmax=min(xmax,360);
%    ymin=max(ymin,0);
%    ymax=80;
%end
%
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
%[C,h]=m_contourf(x(kmin:kmax),y(imin:imax),z(imin:imax,kmin:kmax),15);
cmax=(max(max(z(imin:imax,kmin:kmax))));
cint=cmax/15;
if verLessThan('matlab', '8.5.0') && strcmp(wfill,'y') %  Apply this to pre 2015 (R2014b) matlab versions
    [C,h]=m_contourf(x(kmin:kmax),y(imin:imax),z(imin:imax,kmin:kmax),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
    %cmin=(min(min(z(imin:imax,kmin:kmax))));
    %[C,h]=m_contour(x(kmin:kmax),y(imin:imax),z(imin:imax,kmin:kmax),cint:cint:cmax) % This was changed on 2/8/2017
    %
    %h1=m_contour(x(kmin:kmax),y(imin:imax),z(imin:imax,kmin:kmax),cint:cint:cmax);
    h1=m_contour(x(kmin:kmax),y(imin:imax),z(imin:imax,kmin:kmax));
    if strcmp(wfill,'y')
        %z(z<0.08*max(max(z)))=NaN;
        h=m_pcolor(x(kmin:kmax),y(imin:imax),z(imin:imax,kmin:kmax));
        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
%
if isempty(storm)
    datea(1)=upper(datea(1));
    tit1=datea;
else
    tit1=[storm(1:2) ' ' storm(3:4) ' ' datea(1:4) ' ' datea(5:6) ' ' datea(7:8)...
    ' ' datea(9:10) ' GMT'];
end   
tita='Track Density';
tittot={tita;tit1};
title(tittot,'fontsize',12,'fontweight','bold','interpreter','none')
%
if strcmp(bas,'MT') || max(latstore(:,1)) >= 55  %  For Mediterranean basin, and polar lows reset longitudes
    load temp longstore
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
%
% Clean up
%
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 