%
% This script contours the number density of best track genesis points,  
% defined as the first point of each track where wind speeds exceed startv,
% filtered (if desired) to include only those events whose winds speeds 
% exceed some threshold, peakv, and filtered in the same way as the 
% synthetic tracks. Gendensity must be run first. This requires 
% installation of the "m_map" routines.
% Copyright 2009 WindRiskTech L.L.C.
% Last modified April, 2024 to take advantage of "ismembertol"
%------------------------------------------------------------------------
%
if exist('shape','var') == 0 || exist('xmax','var') == 0
    f=msgbox('Need to run gendensity before running this script');
    return
end   
%
params    %  Load parameters
%
figure
clear latscat longscat z2
%
pifac=acos(-1)/180;
x=xmin:gresbest:xmax;
y=ymin:gresbest:ymax;
coslati=1./(0.001+cos(pifac.*y));
%
bestproc   %  Process best tracks
%
[nn,m]=size(vbest);
%
[~,jmax]=min(vbest,[],2);
jmax=jmax-1;
%
[~,jmin]=min(max((startv-vbest),0),[],2);
%
yearbeg=min(bestyears);
yearend=max(bestyears);
%
n=0;
latscat=zeros(1,nn);
longscat=zeros(1,nn);
for i=1:nn
    if vmaxb(i) >= peakv 
      n=n+1;
      longscat(n)=longbest(i,jmin(i));
      latscat(n)=latbest(i,jmin(i));
    end   
end
longscat(n+1:nn)=[];
latscat(n+1:nn)=[];
%
an=max(size(y));
am=max(size(x));
z2=zeros(an,am);
%
[latl,ay]=ismembertol(latscat,y,1,'DataScale',0.5*gresbest);
[longl,ax]=ismembertol(longscat,x,1,'DataScale',0.5*gresbest);
ax=max(ax,1);
ay=max(ay,1);
N=max(size(latl));
ga=latl.*longl./(0.001+cos(pifac*latl));  
for i=1:N
    z2(ay(i),ax(i))=z2(ay(i),ax(i))+ga(i);
end
qp=max(size(bestyears));
z2=z2/(qp*gresbest^2);
%
%cmax=(max(max(z2(imin:imax,kmin:kmax))));
%cint=cmax/15;
if strcmp(bas,'GB')
    projection=gproject;
end    
m_proj(projection,'long',[xmin xmax],'lat',[ymin ymax]);
hold on
%[C1,h1]=m_contour(x(kmin:kmax),y(imin:imax),z2(imin:imax,kmin:kmax),cint:cint:cmax);
[C1,h1]=m_contour(x,y,z2,30);
colormap(wmap)
if strcmp(wfill,'y')
    hold on
    %h=m_pcolor(x(kmin:kmax),y(imin:imax),z2(imin:imax,kmin:kmax));
    h=m_pcolor(x,y,z2);
    a=colormap(wmap);
    for i=1:8
        a(i,:)=oceancolor;
    end
    colormap(a)
    set(h,'EdgeAlpha',0,'FaceAlpha',wtrans);
    set(h,'FaceColor','interp');
    hold off
end
cb=colorbar;
set(cb,'fontweight',axisfontweight)
m_coast('patch',landcolor,'edgecolor','none');
%m_gshhs_l('patch',landcolor,'edgecolor','none');
if char(pstates) == char('y')
    m_states(xmin, xmax, ymin, ymax,'color',stcolor)
end 
m_gridm('fontname',axisfont,'fontsize',axisfontsize,'fontweight',axisfontweight, ...
    'linestyle',gridline,'backgroundcolor',oceancolor)
hold off
titbest=['Best track genesis 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')
    try
        delete(hButton)
    catch    
    end    
    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 