%
%  This script calculates and graphs the distribution of 
%  wind speed for a given storm at a specified time.
% Contours are displayed in Google Earth.
%
%       K. Emanuel  1/20/2011
%       Copyright 2011 WindRiskTech L.L.C.
%       Last modified January, 2019
%
%------------------------------------------------------------------------
trans=0.5;   % Transparency of filled contours
%------------------------------------------------------------------------
%
try
    kwindfield=kml('getempfile');
catch err
    errordlg('The free KML toolbox must be installed into MATLAB to use this script')
    return
end    
%
params    %  Load parameters
%
if strcmp(randfac,'y')
    if exist('tmr2','var')
        stop(tmr2)
        delete(tmr2)
        clear tmr2
    end   
    tmr2=timer;
    set(tmr2,'startdelay',3);
    set(tmr2,'timerfcn','try; delete(hms); catch err; end')  
    hms=msgbox('Caution: Random draw for wind radii turned on; results will differ with each call');
    start(tmr2)
end    
%
if exist('shearstore','var') == 0
    shape='circ';
    load('temp.mat')
end    
[nrm,mrm]=size(rmstore);
rfac=magfac.*(1+zeros(nrm,mrm));
%  
% Randomize radii of maximum wind
%
if strcmp(randfac,'y')
    rfac=magfac.*outer(nrm,mrm);
end
fs3='sorted';
load(fs3)
%
if exist('nt','var') == 0
    if isempty(storm)
        prompt = {'Order by maximum at POI (1), maximum near POI (2), lifetime maximum (3), or track number (4)?', ...
        'Track number'};
        dlg_title = 'Track Parameters';
        num_lines = 1;
        answer = inputdlg(prompt,dlg_title,num_lines);
        if isempty(answer) || isempty(answer{1})
            return
        end    
        nswitch=str2double(answer{1});
        if (strcmp(shape,'poly') || city_radius >= 2000) && nswitch == 1
            warndlg('Warning: Winds at POI not meaningful for this set; ordering by maximum near POI')
            nswitch=2;
        end                
        nttemp=str2double(answer{2});
    else    
       nswitch=3;
       prompt={'Track Number'};
       dlg_title='Track Specification';
       num_lines = 1;
       answer = inputdlg(prompt,dlg_title,num_lines);
       if isempty(answer) || isempty(answer{1})
            return
       end
       nttemp=str2double(answer{1});
    end  
    if nswitch == 1
        nt=nsortc(nttemp);
    elseif nswitch == 2
       nt=nsortp(nttemp);
    elseif nswitch == 3   
       nt=nsort(nttemp);
    else
       nt=nttemp;
    end   
end
%
prompt = {'Month (1-12)','Day of month', ...
    'GMT hour (may include decimal digits)'};
dlg_title = 'Date & Time Specification';
num_lines = 1;
answer = inputdlg(prompt,dlg_title,num_lines);
if isempty(answer) || isempty(answer{1})
    return
end    
monthplot=str2double(answer{1});
dayplot=str2double(answer{2});
hourplot=str2double(answer{3});
%
ntstring=num2str(nt);
%
[~,jmaxd]=min(vstore,[],2);
jmaxd=jmaxd-1;
%
[n,m]=size(vstore);
%
clear timev V
dum=zeros(1,jmaxd(nt));
V=[dum+2000; monthstore(nt,1:jmaxd(nt)); daystore(nt,1:jmaxd(nt));...
    hourstore(nt,1:jmaxd(nt)); dum; dum];
V=transpose(V);
timev=datenum(V);
timeplot=datenum([2000 monthplot dayplot hourplot 0 0]);
[~,jplot]=min(max(timeplot-timev,0));
jplot=max(jplot,2);
weight=(timeplot-timev(jplot-1))./(timev(jplot)-timev(jplot-1));
%
latstorm=weight.*latstore(nt,jplot)+(1-weight).*latstore(nt,jplot-1);
longstorm=weight.*longstore(nt,jplot)+(1-weight).*longstore(nt,jplot-1);
vstorm=weight.*vstore(nt,jplot)+(1-weight).*vstore(nt,jplot-1);
rmstorm=(weight.*rfac(nt,jplot).*rmstore(nt,jplot)+(1-weight).* ...
    rfac(nt,jplot-1).*rmstore(nt,jplot-1));
if exist('vsestore','var')
    vsestorm=weight.*vsestore(nt,jplot)+(1-weight).*vsestore(nt,jplot-1);
    rmsestorm=(weight.*rfac(nt,jplot).*rmsestore(nt,jplot)+(1-weight).* ...
        rfac(nt,jplot-1).*rmsestore(nt,jplot-1));
end
utstorm=weight.*uinc(nt,jplot)+(1-weight).*uinc(nt,jplot-1);
vtstorm=weight.*vinc(nt,jplot)+(1-weight).*vinc(nt,jplot-1);
%
pifac=acos(-1)/180;
dfac=60.*1.852;
%
if strcmp(bound,'auto')
    bxmin=floor(longstorm-deltax);
    bxmax=ceil(longstorm+deltax);
    bymin=floor(latstorm-deltay);
    bymax=ceil(latstorm+deltay);
end
%
x=bxmin:dellatlong:bxmax;
y=bymin:dellatlong:bymax;
sx=max(size(x));
sy=max(size(y));
%
if strcmp(wreduc,'y')
    logfac=log(wheight/500)/0.35; % Factor to Reduce gradient wind to 10 m wind
    %
    % Load neutral drag coefficients
    %
    load C_Drag500 cd
    load bathymetry bathy
    mincd=min(min(cd));
    cd(bathy<0)=mincd; % This statement added to set ocean point drag coefficients to 
                       % their proper value; this treats a few bogus points
                       % in the drag coefficient data set resulting from
                       % ECMWF's interpolation from a 1 x 1 grid
    clear bathy
    rat=1./(1+sqrt(mincd)*logfac); % Factor to ensure correct wind speeds over water
    %
    % Interpolate drag coefficient to POI
    %
    cdrag=zeros(sx,sy);
    for i=1:sx
        for j=1:sy
            ib=1+floor(4.*x(i));
            ibp=ib+1;
            if ibp > 1440
                ibp=1;
            end    
            jb=1+floor(4.*(y(j)+90));
            b1=cd(ib,jb);
            b2=cd(ib,jb+1);
            b3=cd(ibp,jb);
            b4=cd(ibp,jb+1);
            %
            %  This block interpolates drag coefficient to plong, plat
            %
            dely=1+4.*(y(j)+90)-jb;
            delx=1+4.*x(i)-ib;
            d1=(1.-delx).*(1.-dely);
            d2=dely.*(1.-delx);
            d3=delx.*(1.-dely);
            d4=delx.*dely;
            cdrag(i,j)=1./(d1./b1+d2./b2+d3./b3+d4./b4); % Conservative estimate; weights small roughness
            %
        end
    end
end 
%
% Calculate distance of each POI from track
%
dx=zeros(sx,sy);
dy=dx;
for i=1:sx
    for j=1:sy
        dx(i,j)=dfac.*cos(pifac.*latstorm).*(x(i)-longstorm);
        dy(i,j)=dfac.*(y(j)-latstorm);
    end
end    
r=sqrt(dx.^2+dy.^2);
if exist('vsestore','var')
    vcirc=windprofilem(vstorm,rmstorm,vsestorm,rmsestorm,r,wprofile);
else
    vcirc=windprofiles(vstorm,rmstorm,r,wprofile);
end 
vcirc=vcirc.*latstorm./(abs(latstorm)+1e-8);
un=-vcirc.*dy./max(r,1e-8)+utstorm;
vn=vcirc.*dx./max(r,1e-8)+vtstorm;
%
rfac=r./sqrt(r.^2+4.*radcity.^2);
vspeed=sqrt(un.^2+vn.^2)-rfac.*sqrt(utstorm.^2+vtstorm.^2);
%
if strcmp(wreduc,'y')
    cdfac=max(1+sqrt(cdrag).*logfac,0);
    vspeed=vspeed.*rat.*cdfac;
end
vspeed=max(vspeed,0);
vspeed=wunitsfac*transpose(vspeed);
%
vvm=max(max(vspeed));
vvmm=10*ceil(0.1*vvm);
vinc2=10;
if vvmm <= 100
    vinc2=5;
end    
z=vspeed;
%
xm=repmat(x,max(size(y)),1);
if min(x) > 200
    xm=xm-360;
end
ym=repmat(y,max(size(x)),1);
ym=ym';
%
kwindfield.contourf(xm,ym,z,'numberOfLevels',20:vinc2:vvmm,'transparency',trans);
kwindfield.contour(xm,ym,z,'numberOfLevels',20:vinc2:vvmm,'lineWidth',2);
%
a=jet;
kwindfield.colorbar([0,max(max(z))],a);
%
if strcmp(plottrack,'y')
    longplot=zeros(1,jplot);
    longplot(1,1:jplot-1)=longstore(nt,1:jplot-1);
    longplot(1,jplot)=longstorm;
    latplot=zeros(1,jplot);
    latplot(1,1:jplot-1)=latstore(nt,1:jplot-1);
    latplot(1,jplot)=latstorm;
    kwindfield.plot(longplot,latplot,'lineWidth',4,'lineColor','FF000000');
    sd=(strcat('FF',rgbconv(fliplr(ccol))));
    xplot=longplot(1,jplot);
    if longplot(1,jplot) > 180
        xplot=xplot-360;
    end    
    kwindfield.scatter(xplot,latplot(1,jplot),'iconcolor',sd,'iconscale',1);
end    
clear cd
%
kwindfield.run
%
% Clean up
%
delete('Colorbar.png')
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 