%
% 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. Contours are displayed using Google Earth. 
% Copyright 2009 WindRiskTech L.L.C.
% Last modified April, 2024 to take advantage of "ismembertol"
%------------------------------------------------------------------------
trans=0.5;   % Transparency of filled contours
%------------------------------------------------------------------------
%
try
    ktrackdensity=kml('getempfile');
catch err
    errordlg('The free KML toolbox must be installed into MATLAB to use this script')
    return
end    
%
params    %  Load parameters
%
clear latscat longscat y x z vmax2
if exist('shearstore','var') == 0
    shape='circ';
    load('temp.mat')
end
load('sorted')
[nn,m]=size(vstore);   
if exist('yearstore','var') == 0
    yearstore=2000+zeros(1,nn);
    freqyear=freq;
end
%
pi=acos(-1);
pifac=pi/180;
%
latmask=min(ceil(abs(latstore)),1);
%if strcmp(bas,'MT') || max(latstore(:,1)) >= 55   %  For Mediterranean basin, or polar lows rectify longitudes
%    for i=1:nn,
%        for j=1:m,
%            if longstore(i,j) >= 200
%                longstore(i,j)=longstore(i,j)-360.001;  %#ok<SAGROW> The 0.001 prevents artificial zeros
%            end
%        end    
%    end
%end    
%
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') && max(latstore(:,1)) < 55) || strcmp(bas,'EW')
    xmin=max(xmin,200);
    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;
an=max(size(y));
am=max(size(x));
z=zeros(an,am);
%
vmax2(:,1)=ceil(min(max(vmax-peakv,0),1));  
vmax3=repmat(vmax2,[1,m]);
freqmask=freqyear(1+yearstore-min(yearstore));
freqmask2=repmat(freqmask,[1,m]);
dum=nonzeros(vmax3'.*latstore'.*latmask');
[latl,ay]=ismembertol(dum,y,1,'DataScale',0.5*mres);
dum=nonzeros(vmax3'.*(longstore+0.001)'.*latmask');
[longl,ax]=ismembertol(dum,x,1,'DataScale',0.5*mres);
clear dum
vl=ceil(min(max((nonzeros(vmax3'.*vnet'.*latmask')-minv),0),1));
%
ax=max(ax,1);
ay=max(ay,1);
%
N=max(size(vl));
den=freqmask2(1:N)'.*vl(1:N).*latl(1:N).*longl(1:N)./cos(pifac*latl(1:N));
%
for i=2:N
    if ay(i) ~= ay(i-1) || ax(i) ~= ax(i-1)
        z(ay(i),ax(i))=z(ay(i),ax(i))+den(i);
    end    
end
%
z=z./(nn*mres^2);
%
% 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
%
%[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;
xm=repmat(x(kmin:kmax),imax-imin+1,1);
if strcmp(bas,'AL') || strcmp(bas,'EP') || strcmp(bas,'SA')
    xm=xm-360;
end
ym=repmat(y(imin:imax),kmax-kmin+1,1);
ym=ym';
%
ktrackdensity.contourf(xm,ym,z(imin:imax,kmin:kmax),'numberOfLevels',ceil(cmax/cint),'transparency',trans);
ktrackdensity.contour(xm,ym,z(imin:imax,kmin:kmax),'numberOfLevels',ceil(cmax/cint));
%
a=jet;
ktrackdensity.colorbar([0,max(max(z))],a);
%
if strcmp(bas,'MT') || max(latstore(:,1)) >= 55  %  For Mediterranean basin, and polar lows reset longitudes
    load temp longstore
end  
%
ktrackdensity.run
%
% Clean up
%
delete('Colorbar.png')
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 