%
%  This script creates maps of the power dissipation index. The quantity 
%  mresbest specifies the resolution of the maps. Comment out line 
%  vmax=min(vmax,1) below if trackdensity desired. The routines
%  "m_map" must be installed. 
%  Copyright 2009 WindRiskTech L.L.C.; 
%  Last modified April, 2024 to take advantage of "ismembertol"
%-------------------------------------------------------------------------
%
if exist('bas','var') == 0 || exist('x','var') == 0
    f=msgbox('Need to run pdimapmt before running this script');
    return
end   
params    %  Load parameters
%
bestproc   %  Process best tracks
%
if strcmp(bas,'GB')
   projection=gproject;
end  
%
yearbeg=min(bestyears);
yearend=max(bestyears);
yearsize=max(size(bestyears));
fbest=1/yearsize; 
%
nx=max(size(x));
ny=max(size(y));
%
pi=acos(-1);
pifac=pi/180;
fac=2*3600*(1852/3600)^3; % Factor to convert power dissipation rate to m^3/s^3
fac=fac/mres^2;
%
[n,m]=size(vbest);
%
latmask=min(ceil(abs(latbest)),1);
longbest=min(longbest,359.99);
if strcmp(bas,'MT') || strcmp(bas,'AL') %  For Mediterranean and Atlantic basins, rectify longitudes
    longbest = wrapTo180(longbest);
end 
%
[~,jmax]=min(vbest,[],2);
jmax=jmax-1;
%
%
dum=nonzeros(latmask.*(longbest+0.01))-0.01;
[longl,ax]=ismembertol(dum,x,1,'DataScale',0.5*mres);
clear dum
%
dum=nonzeros(latmask.*latbest);
[latl,ay]=ismembertol(dum,y,1,'DataScale',0.5*mres);
%
vmr=nonzeros(latmask.*(vbest+0.01))-0.01;
ax=max(ax,1);
ay=max(ay,1);
%
z2=zeros(nx,ny);
N=max(size(latl));
pdi=vmr.^3.*latl.*longl./cos(pifac*latl);
for i=1:N
    z2(ax(i),ay(i))=z2(ax(i),ay(i))+pdi(i);
end
% 
z2=1e-4.*z2.*fac/(1+yearend-yearbeg);
z2=z2';
%
gh=figure('Visible','off');
axesm('MapProjection',projection,'MapLatLimit',[ymin ymax],'MapLonLimit',[xmin xmax], ...
     'frame','on','ffacecolor',oceancolor,'fontname',axisfont,'fontsize',axisfontsize, ...
     'fontweight',axisfontweight,'labelrotation','on')
delx=floor((xmax-xmin)/5);
dely=floor((ymax-ymin)/5);
gridm('mlinelocation',delx,'plinelocation',dely,'mlabellocation',delx,'plabellocation',dely,'mlabelparallel','south')
if strcmp(gridline,'none')
    gridm('off')
end
plabel('fontweight','bold'); mlabel('fontweight','bold');
axis off
R = georasterref('RasterSize', size(z2(imin:imax,kmin:kmax)), ...
   'Latlim', [ymin ymax], 'Lonlim', [xmin xmax]);
%
geoback(xmin, xmax, ymin, ymax)
mlabel('south')
%
colormap(wmap)
if strcmp(wfill,'y')
    a=colormap(wmap);
    for i=1:8
        a(i,:)=oceancolor;
    end
    colormap(a)
    [C,h]=contourm(z2(imin:imax,kmin:kmax),R,12,'fill','on');
    hPatch = findobj(h,'Type','patch');
    set(hPatch,'facealpha',wtrans); 
    qh=max(size(hPatch));
    set(hPatch(qh,1),'facealpha',0)
else
    [C1,h1]=contourm(z2(imin:imax,kmin:kmax),R,15);
end    
%
set(gh,'visible','on')
%
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)
titbest=['Best Track Power Dissipation Density from ' num2str(yearbeg) ' to ' num2str(yearend)];
title(titbest,'fontsize',12,'fontweight','bold','interpreter','none')
tightmap
%
% 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 