%
%  This script calculates and graphs the distribution of 
%  wind speed for a given storm at a specified time.
%
%       K. Emanuel  1/20/2011
%       Copyright 2011 WindRiskTech L.L.C.
%       Last modified January, 2019
%
%------------------------------------------------------------------------
%
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    
%
clf('reset')
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    
%
clf
axesm('MapProjection',projection,'MapLatLimit',[bymin bymax],'MapLonLimit',[bxmin bxmax], ...
     'frame','on','ffacecolor',oceancolor,'fontname',axisfont,'fontsize',axisfontsize, ...
     'fontweight',axisfontweight,'labelrotation','on')
delx=floor((bxmax-bxmin)/5);
dely=floor((bymax-bymin)/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(vspeed), ...
   'Latlim', [bymin bymax], 'Lonlim', [bxmin bxmax]);
%
geoback(bxmin, bxmax, bymin, bymax)
%
z=vspeed;
[C,h]=contourm(vspeed,R,20:vinc2:vvmm,'linewidth',2);
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;
    g1=geoshow(latplot,longplot);
    set(g1,'color','k')
    set(g1,'linewidth',3)
    g2=geoshow(latplot(1,jplot),longplot(1,jplot),'displaytype','point');
    set(g2,'marker','.','markersize',18,'markeredgecolor','r')
end  
%thismap = get(gcf, 'Colormap');
%thismap(1,:)=[1 1 1];
%colormap(thismap)
 
%hPatch = findobj(h,'Type','patch');
%set(hPatch,'facealpha',wtrans); 
%set(findobj(h,'FaceColor',[1 1 1]), 'FaceAlpha',0)
cb=colorbar;
%cbIm = findobj(cb,'Type','image');
%alpha(cbIm,wtrans)
set(cb,'fontname',axisfont,'fontweight','bold')
set(get(cb,'ylabel'),'String',wunitslabel,'fontweight','bold','rotation',0,'horizontalalignment','left');
yvec=get(cb,'ylim');
yvec(1,1)=20;
set(cb,'ylim',yvec)
if isempty(storm)
    datea(1)=upper(datea(1));
    tit1=datea;
    adate=datestr(timeplot,'mmmm dd, HH:MM');
    tit1={tit1;['Track number ' ntstring ', ' adate ' GMT']};
else
    tit1=[storm(1:2) ' ' storm(3:4) ' ' datea(1:4) ' ' datea(5:6) ' ' datea(7:8)...
    ' ' datea(9:10) ' GMT'];
    titstring2=['Track number ' ntstring  ];
    tit2=titstring2;
    tit1=[tit1 '  ' tit2];
end    
title(tit1,'fontsize',12,'fontweight','bold','interpreter','none')
tightmap
set(gcf,'Renderer','zbuffer')
%
% Add button to access documentation of script
%
if strcmp(Docdisp,'y')
    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
clear cd
%
% Clean up
%
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 