%
%  This script calculates and graphs the distribution of 
%  rain rate, in mm/hr, for a given storm at a specified time.
%
%       December 2011
%       Copyright 2011 WindRiskTech L.L.C.
%       Revised December 2020 to remove cos(latitude) factor from hy and to
%       allow for higher resolution topography
%
%------------------------------------------------------------------------
%
params    %  Load parameters
%
if strcmp(randfac,'y')
    if exist('tmr7','var')
        stop(tmr7)
        delete(tmr7)
        clear tmr7
    end   
    tmr7=timer;
    set(tmr7,'startdelay',6);
    set(tmr7,'timerfcn','try; delete(hms); catch err; end')
    hms=msgbox('Caution: Random draw for wind radii turned on; results will differ with each call');
    start(tmr7)
end 
%
clear rain
%
clf('reset')
if exist('shearstore','var') == 0
    shape='circ';
    clear T600store u850store v850store
    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 strcmpi(topores,'high')
    load bathymetry_high
else
    load bathymetry
end    
[ntopo,~]=size(bathy);
topores=360/ntopo;
toporesi=1/topores;
%
pifac=acos(-1)/180;
dfac=60.*1.852;
knotfac=1852./3600;
sfac=1./(topores.*60.*1852);
%
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
ntstring=num2str(nt);
%
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});
%
[n,m]=size(vstore);
[~, jmaxd]=min(vstore,[],2);
jmaxd=jmaxd-1;
jmaxd=transpose(jmaxd);
%
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]);
[dum2,jplot]=min(max(timeplot-timev,0));
jplot=max(jplot,2);
jplot=min(jplot,399);
%
latstorm=latstore(nt,jplot-1:jplot+1);
longstorm=longstore(nt,jplot-1:jplot+1);
vstorm=vstore(nt,jplot-1:jplot+1);
rmstorm=rfac(nt,jplot-1:jplot+1).*(rmstore(nt,jplot-1:jplot+1));
vsestorm=vsestore(nt,jplot-1:jplot+1);
rmsestorm=rfac(nt,jplot-1:jplot+1).*(rmsestore(nt,jplot-1:jplot+1));
utstorm=ut(nt,jplot-1:jplot+1);
vtstorm=vt(nt,jplot-1:jplot+1);
ush=zeros(1,3);
vsh=ush;
if exist('u850store','var')
    vdrift=1.5.*3600./1852;
    vdrift=vdrift.*latstore(1,1)./(abs(latstore(1,1))+1e-8);
    u850storm=u850store(nt,jplot-1:jplot+1);
    v850storm=v850store(nt,jplot-1:jplot+1);
    ush(1,:)=5.*knotfac.*(utstorm-u850storm);
    vsh(1,:)=5.*knotfac.*(vtstorm-vdrift.*cos(pifac.*latstorm)-v850storm);
end 
%
if char(bound) == char('auto')
    bxmin=floor(longstorm(1,2)-deltax);
    bxmax=ceil(longstorm(1,2)+deltax);
    bymin=floor(latstorm(1,2)-deltay);
    bymax=ceil(latstorm(1,2)+deltay);
end
%
x=bxmin:dellatlong:bxmax;
y=bymin:dellatlong:bymax;
sx=max(size(x));
sy=max(size(y));
%  This statement added because in pointwshort it is assumed that if sx=sy
%  POI's are a vector rather than a grid.
if sx == sy
    sx=sx+1;
    x=[x bxmax+dellatlong];
end  
%
h=zeros(sx,sy);
hx=zeros(sx,sy);
hy=zeros(sx,sy);
m_to_mm=1000;
rowa_over_rowl=0.00117;
%
bathy=max(bathy,-1);
dhdx=sfac*(circshift(bathy,-1,1)-circshift(bathy,0,1));
dhdy=sfac*(circshift(bathy,-1,2)-circshift(bathy,0,2));
%
for i=1:sx
    plong=x(i);
    for j=1:sy
        plat=y(j);
        %
        ib=1+floor(toporesi.*plong);
        ibp=ib+1;
        if ibp > ntopo
            ibp=ibp-ntopo;
        end    
        ibs=1+floor(toporesi.*plong-0.5);
        ibsp=ibs+1;
        plongs=plong;
        if ibs < 0
            ibs=ntopo;
            plongs=plong+360;
        end    
        if ibsp > ntopo
            ibsp=ibsp-ntopo;
        end
        jb=1+floor(toporesi.*(plat+90));
        jbs=1+floor(toporesi.*(plat+90)-0.5);
        b1=bathy(ib,jb);
        b2=bathy(ib,jb+1);
        b3=bathy(ibp,jb);
        b4=bathy(ibp,jb+1);
        dely=1+toporesi.*(plat+90)-jb;
        delx=1+toporesi.*plong-ib;
        d1=(1.-delx).*(1.-dely);
        d2=dely.*(1.-delx);
        d3=delx.*(1.-dely);
        d4=delx.*dely;
        h(i,j)=exp(d1.*log(b1+11)+d2.*log(b2+11)+d3.*log(b3+11)+d4.*log(b4+11))-11;
        %
        b1=dhdx(ibs,jbs);
        b2=dhdx(ibs,jbs+1);
        b3=dhdx(ibsp,jbs);
        b4=dhdx(ibsp,jbs+1);
        dely=0.5+toporesi.*(plat+90)-jbs;
        delx=0.5+toporesi.*plongs-ibs;
        d1=(1.-delx).*(1.-dely);
        d2=dely.*(1.-delx);
        d3=delx.*(1.-dely);
        d4=delx.*dely;
        hx(i,j)=(b1.*d1+b2.*d2+b3.*d3+b4.*d4)./cos(pifac.*plat);
        %
        b1=dhdy(ibs,jbs);
        b2=dhdy(ibs,jbs+1);
        b3=dhdy(ibsp,jbs);
        b4=dhdy(ibsp,jbs+1);
        hy(i,j)=(b1.*d1+b2.*d2+b3.*d3+b4.*d4);    
        %
    end
end  
if exist('T600store','var')
    clear q900 q600
    [ q900, q600 ] = qs900b( T600store(nt,jplot), vstore(nt,jplot) );
end    
w = pointwfield(latstorm,longstorm,vstorm,rmstorm,vsestorm,rmsestorm,utstorm,vtstorm,ush,vsh, ...
    y,x,h,hx,hy);
temp=eprecip.*m_to_mm.*3600.*rowa_over_rowl.*q900.*max(w(1,2,:,:)-wrad,0);
rain(:,:)=temp(1,1,:,:);
rain=transpose(rain);
clear temp
%
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(rain), ...
   'Latlim', [bymin bymax], 'Lonlim', [bxmin bxmax]);
%
geoback(bxmin, bxmax, bymin, bymax)
%
[C,h]=contourm(runitsfac*rain,R,runitsfac*raincntrs,'linewidth',2);
%thismap = get(gcf, 'Colormap');
%thismap(1,:)=[1 1 1];
%colormap(thismap)
%
if plottrack == 'y'
    g1=plotm(latstore(nt,1:jplot),longstore(nt,1:jplot));
    set(g1,'color','k')
    set(g1,'linewidth',3)
    g2=plotm(latstore(nt,jplot),longstore(nt,jplot),'marker','.','markersize',18,'markeredgecolor','r');
end    
%alpha(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')
ystr=char(strcat(runitslabel,'/hr'));
set(get(cb,'ylabel'),'String', ystr,'fontweight','bold','rotation',0,'horizontalalignment','left');
yvec=get(cb,'ylim');
yvec(1,1)=0.1*yvec(1,2);
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')
if exist('tmr7','var')
    stop(tmr7)
    delete(tmr7)
    clear tmr7
end
clear hms
tightmap
set(gcf,'Renderer','zbuffer')
z=runitsfac*rain;
%
% 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
%
% Clean up
%
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 