%
% 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.
% K. Emanuel (First written June, 2009)
% Copyright 2009 WindRiskTech L.L.C.
% Last modified October, 2010
%
%------------------------------------------------------------------------
%  Options
%
method='globe';  % "local" for filtering by max wind at POI; "globe" for 
                 %  filtering by storm lifetime maximum wind
peakv=35;        %  Maximum wind speed (knots) that must be exceeded for
                 %     event to be included
startv=30;       %  First point wind speed (knots) must exceed this value
                 %
transfac=0.6;    %  Factor of translation speed to add to circular winds
                 %
mres=2.5;        %  Map resolution (degrees)
%
%------------------------------------------------------------------------                
clf
clear latscat longscat alat along npoints ut
shape='circ';
load('temp.mat')
load('sorted')
[nn,m]=size(vstore);
%
pi=acos(-1);
pifac=pi./180;
%
%  Get translations speeds
%
[ut,vt,jmax] = utrans(latstore,longstore);
vtemp=vstore+transfac.*sqrt(ut.^2+vt.^2);
%
if method == 'local'
    [vs, dir, dayk] = pointseries(latstore,longstore,vstore,rmstore, ...
      monthstore,daystore,hourstore,ut,vt,jmax,clat,clong);
    vmax=max(vs,[],2);
else
    [vmax,imax]=max(vtemp,[],2);
end   
pdi=freq.*sum(sum(vtemp.^3))./nn;
%
[dum,jmin]=min(max((startv-vtemp),0),[],2);
if bas == 'AL' 
   m_proj('gall-peters','long',[260 350], ...
      'lat',[0 60]);
   alat=[0:mres:60];
   along=[250:mres:360];
elseif bas == 'EP'
   m_proj('gall-peters','long',[160 285], ...
      'lat',[0 40]);
   alat=[0:mres:40];
   along=[160:mres:285];
elseif bas == 'WP'
   m_proj('gall-peters','long',[100 220], ...
      'lat',[0 55]);
   alat=[0:mres:55];
   along=[100:mres:220];
elseif bas == 'IO'
   m_proj('gall-peters','long',[40 100], ...
      'lat',[0 40]);
   alat=[0:mres:40];
   along=[40:mres:100];
elseif bas == 'SH'
   m_proj('gall-peters','long',[30 240], ...
      'lat',[-40 0]);
   alat=[-40:mres:0];
   along=[30:mres:240];
elseif bas == 'BM'
   m_proj('gall-peters','long',[100 180], ...
      'lat',[-40 0]);
   alat=[-40:mres:0];
   along=[100:mres:180];  
elseif bas == 'SA'
   m_proj('gall-peters','long',[290 390], ...
      'lat',[-50 0]);
   alat=[-50:mres:0];
   along=[290:mres:390];
elseif bas == 'GB'
   m_proj('gall-peters','long',[0 360], ...
      'lat',[-70 70]);
   alat=[-70:mres:70];
   along=[0:mres:360];
end   
an=max(size(alat));
am=max(size(along));
npoints(1:an,1:am)=0.0;
for j=2:an-1,
   cosfac=1./cos(pifac.*alat(j));
   for k=2:am-1,
      for ii=1:nn,
          if vmax(ii) >= peakv
              for jj =jmin(ii):jmax(ii),   
                 if latstore(ii,jj) >= alat(j-1) & latstore(ii,jj) < alat(j+1) & longstore(ii,jj) >= along(k-1) & longstore(ii,jj) < along(k+1)
                    npoints(j,k)=npoints(j,k)+cosfac;
                    break
                 end
              end   
          end
      end
   end
end   
npoints=freq.*npoints./nn;
axes('color',[0.99 0.99 1])
m_contourf(along,alat,npoints,10);
m_coast('patch',[0.45 0.6 0.45],'edgecolor','none');
m_gridm('fontname','ariel','fontsize',9,'fontweight','bold','linestyle','none')
cb=colorbar;
set(cb,'fontweight','bold')
if isempty(storm)
    tit1=datea;
else
    tit1=[storm(1:2) ' ' storm(3:4) ' ' datea(1:4) ' ' datea(5:6) ' ' datea(7:8)...
    ' ' datea(9:10) ' GMT'];
end    
title(tit1,'fontsize',12,'fontweight','bold','interpreter','none')