%
% This script contours the number density of 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. It requires installation of the
% "m_map" routines. Units are events per 1 degree latitude square per year.
% 
% Copyright 2009 WindRiskTech L.L.C.
% Last modified April, 2024 to take advantage of "ismembertol"
%
%------------------------------------------------------------------------
%
params    %  Load parameters
%
clear latscat longscat y x z imin imax kmin kmax
clf('reset')
if exist('shearstore','var') == 0 || exist('vstore','var') == 0
    storm=[];
    shape='circ';
    load('temp.mat')
end   
load sorted
if strcmp(bas,'GB')
    projection=gproject;
end    
[nn,m]=size(vstore);
pifac=acos(-1)./180;
%
[~,jmax]=min(vstore,[],2);
jmax=jmax-1;
%
[~,jmin]=min(max((startv-vnet),0),[],2);
longscat=zeros(1,nn);
latscat=zeros(1,nn);
freqmask=zeros(1,nn);
if exist('yearstore','var') == 0
    yeartemp=2000+zeros(1,nn);
    freqtemp=freq;
else
    yeartemp=yearstore;
    freqtemp=freqyear;
    atemp=ismember((min(yearstore):max(yearstore)),yearstore);
    atemp=cast(atemp,'like',yearstore);
    yearset=nonzeros(atemp.*(min(yearstore):max(yearstore)));
end 
n=0;
for i=1:nn
   if vmax(i) >= peakv
      n=n+1; 
      longscat(n)=longstore(i,jmin(i));
      latscat(n)=latstore(i,jmin(i));
      if (strcmp(bas,'MT')|| max(latstore(:,1)) > 50) && longscat(n) > 200
          longscat(n)=longscat(n)-360;
      end    
      [~,nyy]=ismember(yeartemp(i),yearset);
      freqmask(n)=freqtemp(nyy);
   end   
end
longscat(n+1:nn)=[];
latscat(n+1:nn)=[];
freqmask(n+1:nn)=[];
if strcmp(mapmode,'auto')
    xmin=min(nonzeros(longscat))-dellong;
    xmax=max(longscat)+dellong;
    ymin=min(nonzeros(latscat))-dellat;
    ymax=max(latscat)+dellat;
    xmin=gres*floor(xmin/gres);
    xmax=gres*ceil(xmax/gres);
    ymin=gres*floor(ymin/gres);
    ymax=gres*ceil(ymax/gres);    
else
    xmin=longmin;
    xmax=longmax;
    ymin=latmin;
    ymax=latmax;
end 
if strcmp(bas,'AL') && max(latstore(:,1)) < 55
    xmin=max(xmin,250);
    xmax=min(xmax,360);
    ymin=max(ymin,0);
%elseif strcmp(bas,'GB')
%    xmin=0;
%    xmax=360;
end 
x=xmin:gres:xmax;
y=ymin:gres:ymax;
an=max(size(y));
am=max(size(x));
z=zeros(an,am);
%
[latl,ay]=ismembertol(latscat,y,1,'DataScale',0.5*gres);
[longl,ax]=ismembertol(longscat,x,1,'DataScale',0.5*gres);
ax=max(ax,1);
ay=max(ay,1);
N=max(size(latl));
ga=freqmask.*latl.*longl./cos(pifac*latl);  % Corrected March, 2016
for i=1:N
    z(ay(i),ax(i))=z(ay(i),ax(i))+ga(i);
end
z=z/(nn*gres^2);
zsyn=z;
%
% 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=2:an-1
        for k=2:am-1
            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
%
m_proj(projection,'long',[xmin xmax],'lat',[ymin ymax]);
cmax=(max(max(z(imin:imax,kmin:kmax))));
cint=cmax/20;
[C1,h1]=m_contour(x(kmin:kmax),y(imin:imax),z(imin:imax,kmin:kmax),cint:cint:cmax);
colormap(wmap)
%
if strcmp(wfill,'y')
    hold on
    h=m_pcolor(x(kmin:kmax),y(imin:imax),z(imin:imax,kmin:kmax));
    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)
if strcmp(bas,'MT')
     m_gshhs_l('patch',landcolor,'edgecolor','none');
else     
     m_coast('patch',landcolor,'edgecolor','none');
end     
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)
if isempty(storm)
    datea(1)=upper(datea(1));
    tit1=[num2str(nn) ' Tracks in ' datea];
else
    tit1=[storm(1:2) ' ' storm(3:4) ' ' datea(1:4) ' ' datea(5:6) ' ' datea(7:8)...
    ' ' datea(9:10) ' GMT'];
end    
tita='Genesis Density';
tittot={tita;tit1};
title(tittot,'fontsize',12,'fontweight','bold','interpreter','none')
%
% 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 