%
%  This script creates maps of the power dissipation index. The quantity 
%  mres specifies the resolution of the maps. The results are displayed in
%  Google Earth.
%  November, 2008
%  Copyright 2009 WindRiskTech L.L.C. 
%  Last modified April, 2024 to take advantage of "ismembertol"
%------------------------------------------------------------------------
trans=0.5;   % Transparency of filled contours
%------------------------------------------------------------------------
%
try
    kpdimap=kml('getempfile');
catch err
    errordlg('The free KML toolbox must be installed into MATLAB to use this script')
    return
end
%
params    %  Load parameters
%
if exist('shearstore','var') == 0 || exist('yearstore','var') == 0
    shape='circ';
    load temp
end
load sorted
clear z
%
[n,m]=size(vnet);
%
[~,jmax]=min(vnet,[],2);
jmax=jmax-1;
if exist('yearstore','var') == 0
    yearstore=2000+zeros(1,n);
    freqyear=freq;
end
%
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:n
        for j=1:m
            if longstore(i,j) >= 200
                longstore(i,j)=longstore(i,j)-360.001; %#ok<SAGROW>
            end
        end    
    end
end    
%
pi=acos(-1);
pifac=pi/180;
fac=2*3600*(1852/3600)^3; % Factor to convert power dissipation rate to m^3/s^2
fac=fac/mres^2;
%
if strcmp(bas,'GB')
   projection=gproject;
end   
%
if strcmp(mapmode,'auto')
    mask=ceil(min(max(vnet-30,0),1));  % Added to further restrict area over which computations performed
    xmin=min(nonzeros(mask.*longstore));
    xmax=max(nonzeros(mask.*longstore));
    ymin=min(nonzeros(mask.*latstore));
    ymax=max(nonzeros(mask.*latstore));
    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,'GB')
    xmin=0;
    xmax=360;
    ymin=-87;
    ymax=87;
end 
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));
%
yearmin=min(yearstore);
yearmax=max(yearstore);
qy=yearmax-yearmin+1;
year=yearmin:yearmax;
%
dum=nonzeros(latmask.*(longstore+0.01))-0.01;
[longl,ax]=ismembertol(dum,x,1,'DataScale',0.5*mres);
clear dum
%
dum=nonzeros(latmask.*latstore);
[latl,ay]=ismembertol(dum,y,1,'DataScale',0.5*mres);

%
vmr=nonzeros(latmask.*(vnet+0.01))-0.01;
dum=repmat(max(yearstore',1),[1,m]); 
yr=nonzeros(double(dum).*latmask);
if max(size(freqyear)) > 1
    freqmask=freqyear(yr+1-min(yr));
else
    freqmask=freq;
end    
clear dum
%
[yearl,an]=ismember(yr,year);
ax=max(ax,1);
ay=max(ay,1);
an=max(an,1);
%
z3=zeros(nx,ny,qy);
N=max(size(latl));
pdi=vmr.^3.*freqmask'.*latl.*longl.*yearl./cos(pifac*latl);
for i=1:N
    z3(ax(i),ay(i),an(i))=z3(ax(i),ay(i),an(i))+pdi(i);
end
% 
z(:,:)=sum(z3,3);
z=1e-4.*z.*fac./n;
%
% Calculate new map limits
%
imin=1;
imax=ny;
kmin=1;
kmax=nx;
if strcmp(mapmode,'auto')
    xmin1=xmax;
    xmax1=xmin;
    ymin1=ymax;
    ymax1=ymin;
    kmin=nx;
    kmax=1;
    imin=ny;
    imax=1;
    zcrit=0.01*max(max(z));
    for j=2:ny-1
        for k=2:nx-1
            if z(k,j) > 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')
%    xmin=max(xmin,220+mod(220,mres));
%end
%
cmax=(max(max(z(kmin:kmax,imin:imax))));
cint=cmax/25;
%
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';
%
kpdimap.contourf(xm,ym,z(kmin:kmax,imin:imax)','numberOfLevels',ceil(cmax/cint),'transparency',trans);
kpdimap.contour(xm,ym,z(kmin:kmax,imin:imax)','numberOfLevels',ceil(cmax/cint));
%
a=jet;
kpdimap.colorbar([0,max(max(z))],a);
%
kpdimap.run
%
if strcmp(bas,'MT') %  For Mediterranean basin, rest longitudes
    load temp longstore
end    
%
% Clean up
%
delete('Colorbar.png')
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 