%
%  This script calculates and graphs the distribution of 
%  wind speed for a given storm at a specified time and graphs the result
%  as wind quivers; otherwise, it is like windfield.m
%
%       K. Emanuel  6/28/2022
%       Copyright 2022 WindRiskTech L.L.C.
%-------------------------------------------------------------------------
quiverlinewidth=1.2;  % Governs size of arrows
%-------------------------------------------------------------------------
%
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);
%
u850storm=weight.*u850store(nt,jplot)+(1-weight).*u850store(nt,jplot-1);
v850storm=weight.*v850store(nt,jplot)+(1-weight).*v850store(nt,jplot-1);
%
pifac=acos(-1)/180;
dfac=60.*1.852;
%
if strcmp(bound,'auto')
    bxmin=floor(longstorm-0.75*deltax); % Factor of 0.75 added to zoom in a bit
    bxmax=ceil(longstorm+0.75*deltax);
    bymin=floor(latstorm-0.75*deltay);
    bymax=ceil(latstorm+0.75*deltay);
end
%
x=bxmin:dellatlongq:bxmax;
y=bymin:dellatlongq: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 
%
if strcmpi(plotisobars,'y')
    [psurface,rw] = psurf(vstorm,rmstorm,latstorm,wprofile,vsestorm,rmsestorm);
    pcirc=interp1(rw,psurface,r);
    %
    % Add background pressure field
    %
    f=(pi/(6*3600))*sin(latstorm*pi/180);
    coef=f*1852/(287*300*3600); % f/(R_d*T_s) with factor to convert knots to m/s
    ycoef=60*1852;
    xcoef=ycoef*cos(latstorm*pi/180);
    pcirc=pcirc.*exp(coef*(v850storm*xcoef*repmat(x'-mean(x),[1,sy])- ...
        u850storm*ycoef*repmat(y-mean(y),[sx,1])));
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);
un=-vcirc.*dy./max(r,1e-8)+(1-rfac)*utstorm+rfac*u850storm;
vn=vcirc.*dx./max(r,1e-8)+(1-rfac)*vtstorm+rfac*v850storm;
vspeed=sqrt(un.^2+vn.^2);
%
if strcmp(wreduc,'y')
    cdfac=max(1+sqrt(cdrag).*logfac,0);
    vspeed=vspeed.*rat.*cdfac;
end
vspeed=max(vspeed,0);
%        
vspeed=transpose(vspeed);
%
clf
m_proj(projection,'long',[bxmin bxmax], ...
   'lat',[bymin bymax]);
ld=landcolor;
hold on
if cres == 1
      m_coast('patch',ld,'edgecolor',shorecolor);
elseif cres == 2  
      m_gshhs_l('patch',ld,'edgecolor',shorecolor);
elseif cres == 3  
      m_gshhs_i('patch',ld,'edgecolor',shorecolor);
elseif cres == 4   
      m_gshhs_h('patch',ld,'edgecolor',shorecolor);
end   
%
if strcmp(pstates,'y')
   m_states(bxmin,bxmax,bymin,bymax,'color',stcolor)
end
%
vvm=max(max(vspeed));
vvmm=10*ceil(0.1*vvm);
vinc2=10;
if vvmm <= 100
    vinc2=5;
end    
z=vspeed;
[xn,yn]=meshgrid(x,y);
qm=m_quiver(xn,yn,un',vn');
m_gridm('fontname',axisfont,'fontsize',axisfontsize,'fontweight',axisfontweight, ...
    'linestyle',gridline,'backgroundcolor',oceancolor)
mags = sqrt(sum(cat(2, qm.UData(:), qm.VData(:) ...
).^2, 2));
currentColormap = colormap(jet);
[~, ~, ind] = histcounts(mags, size(currentColormap, 1));
cmap = uint8(ind2rgb(ind(:), currentColormap) * 255);
cmap(:,:,4) = 255;
cmap = permute(repmat(cmap, [1 3 1]), [2 1 3]);
set(qm.Tail, ...
'ColorBinding', 'interpolated', ...
'ColorData', reshape(cmap(1:2,:,:), [], 4).');
set(qm.Head, ...
'ColorBinding', 'interpolated', ...
'ColorData', reshape(cmap(1:3,:,:), [], 4).');   %'
set(qm,'linewidth',quiverlinewidth)
wmax=10*ceil(0.1*max(max(vspeed)));
caxis([0 wmax]);
C=colorbar;
set(C,'fontsize',12)
C.Label.FontSize=12;
C.Label.FontWeight='bold';
C.Label.FontName=gfont;
wunitslab=wunitslabel;
wunitslab(1)=upper(wunitslab(1));
C.Label.String = wunitslab;
C.Label.Rotation=0;
C.Label.HorizontalAlignment='left';
%
if strcmpi(plotisobars,'y')
    [~,CP]=m_contour(xn,yn,pcirc',20,'k');
end
%
plottrack='n';  %  This map does not look too nice with track ploted
%
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;
    m_plot(longplot,latplot,'k','linewidth',3)
    m_plot(longplot(1,jplot),latplot(1,jplot),'.r','markersize',18)
end    
hold off
%{
cb=colorbar;
set(cb,'fontname',gfont,'fontweight','bold')
set(get(cb,'ylabel'),'String',wunitslabel,'fontweight','bold','rotation',0,'horizontalalignment','left');
%}
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  
if strcmpi(plotisobars,'y')
    uu=CP.LevelList;
    titstring3=['Lowest isobar = ',num2str(round(uu(1))),' hPa'];
    tit1=[tit1; titstring3];
end 
title(tit1,'fontsize',12,'fontweight','bold','interpreter','none')
%
% 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 