%
%  This script creates maps of the power dissipation index. The quantity 
%  mres 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 pdimap 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; 
%
%xmin2=xmin-mres;
%xmax2=xmax+mres;
%ymin2=ymin-mres;
%ymax2=ymax+mres;
%
%y=ymin2:mres:ymax2;
%x=xmin2:mres:xmax2; 
%
%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);
if strcmp(bas,'MT') || max(latbest(:,1)) >= 55   %  For Mediterranean basin, or polar lows rectify longitudes
    longbest(longbest>200)=longbest(longbest>200)-360.001;
end   
%
if strcmp(bas,'AL')
    longbest(longbest<100)=longbest(longbest<100)+360;
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);
%
figure
m_proj(projection,'long',[xmin xmax],'lat',[ymin ymax]);
m_gridm('fontname',axisfont,'fontsize',axisfontsize,'fontweight',axisfontweight, ...
    'linestyle',gridline,'backgroundcolor',oceancolor)
%
hold on
if strcmp(bas,'MT')
     m_gshhs_l('patch',landcolor,'edgecolor','none');
else     
     m_coast('patch',landcolor,'edgecolor','none');
end   
if strcmp(pstates,'y')
    m_states(xmin, xmax, ymin, ymax,'color',stcolor)
end  
if verLessThan('matlab', '8.5.0') && strcmp(wfill,'y')  %  Apply this to pre 2015 (R2014b) matlab versions
    [C,h]=m_contourf(transpose(x(kmin:kmax)),transpose(y(imin:imax)),transpose(z2(kmin:kmax,imin:imax)),15);
    thismap = get(gcf, 'Colormap');
    thismap(1,:)=[1 1 1];
    colormap(thismap)
    hold off
    alpha(0.99)
    hPatch = findobj(h,'Type','patch');
    set(hPatch,'facealpha',wtrans); 
    set(findobj(h,'FaceColor',[1 1 1]), 'FaceAlpha',0)
else
    [C1,h1]=m_contour(transpose(x(kmin:kmax)),transpose(y(imin:imax)),transpose(z2(kmin:kmax,imin:imax)),15);
    colormap(wmap)
    if strcmp(wfill,'y')
        h=m_pcolor(transpose(x(kmin:kmax)),transpose(y(imin:imax)),transpose(z2(kmin:kmax,imin:imax)));
        a=colormap(wmap);
        for i=1:4
            a(i,:)=oceancolor;
        end
        colormap(a)
        set(h,'EdgeAlpha',0,'FaceAlpha',wtrans);
        set(h,'FaceColor','interp');
    end
    hold off
end    
%
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')
%
% 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 