%
%  This script calculates and graphs the distribution of accumulated
%  precipitation, in mm, for a given individual storm.
%       January, 2009
%       Copyright 2009 WindRiskTech L.L.C.
%       Revised March, 2020
%       Revised December 2020 to remove cos(latitude) factor from hy and to
%       allow for higher resolution topography
%       Revised November, 2021 to handle storms that cross the Greenwich
%       meridion
%       Revised 26 December 2023 to allow for filled contours only land only
%-------------------------------------------------------------------------
%
params    %  Load parameters
%
if strcmp(randfac,'y')
    if exist('tmr7','var')
        stop(tmr7)
        delete(tmr7)
        clear tmr7
    end   
    tmr7=timer;
    set(tmr7,'startdelay',2);
    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 lat long v rm vse rmse ush vsh jmax y x h hx hy
clf('reset')
if exist('shearstore','var') == 0
    shape='circ';
    clear T600store u850store v850store
    load('temp.mat')
else
    load('temp.mat','vt')
    load('temp.mat','ut')    
    if strcmp(shape,'poly')
        load temp polytype passe xint yint jint kint kfrac 
    end    
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;
%
clear rain v
%
pifac=acos(-1)/180;
dfac=60.*1.852;
knotfac=1852./3600;
sfac=1./(topores.*60.*1852);
c1=0.622.*2.5e6./1e3;
c2=2.5e6.^2./(1e3.*461);
m_to_mm=1000;
rowa_over_rowl=0.00117;
%
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);
if bxmin < 0
    bxmin=bxmin+360;
end
if bxmax < 0
    bxmax=bxmax+360;
end    
%
[n,m]=size(vstore);
%
clear utd vtd
latsize=length(nonzeros(latstore(nt,:)));
utd(1,:)=ut(nt,1:latsize);
vtd(1,:)=vt(nt,1:latsize);
ush=0.*utd;
vsh=ush;
vdrift=1.5.*3600./1852;
vdrift=vdrift.*latstore(1,1)./(abs(latstore(1,1))+1e-8);
if exist('u850store','var') == 1
    ush(1,:)=5.*knotfac.*(utd-u850store(nt,1:latsize));
    vsh(1,:)=5.*knotfac.*(vtd-vdrift.*cos(pifac.*latstore(nt,1:latsize))- ...
        v850store(nt,1:latsize));
end    
lat(1,:)=latstore(nt,1:latsize);
long(1,:)=longstore(nt,1:latsize);
  
v(1,:)=vstore(nt,1:latsize);
vse(1,:)=vsestore(nt,1:latsize);
%  
% Scale and randomize radii of maximum wind
%
[nrm,mrm]=size(rmstore);
[~, jmaxd]=min(vstore,[],2);
%jmaxd(vstore(:,end)~=0)=m;
jmaxd(vstore(:,end)>3 & vstore(:,end) ~= vstore(:,end-1))=m;
jmaxd=jmaxd-1;
jmaxd=transpose(jmaxd);
rfac=zeros(nrm,mrm)+1;
if strcmp(randfac,'y')
    rfac=outer(nrm,mrm);
end  
%
rm(1,:)=transpose(nonzeros(magfac.*rfac(nt,:).*rmstore(nt,:)));
nrm=max(size(rm));
rmse(1,:)=magfac.*rfac(nt,1:nrm).*rmsestore(nt,1:nrm);
if max(long) - min(long) > 330  && mean(nonzeros(longstore)) < 50 % These lines pertain to Mediterranean storms
    long=wrapTo180(long);
    utd(long<0.5)=[];
    vtd(long<0.5)=[];
    lat(long<0.5)=[];
    v(long<0.5)=[];
    rm(long<0.5)=[];
    vse(long<0.5)=[];
    rmse(long<0.5)=[];
    long(long<0.5)=[];
end  
jmax=max(size(long));
for i=1:jmaxd(nt)
    if long(1) > 200 && long(i) < 50
        long(i)=long(i)+360;
    end
end    
%
if strcmp(bound,'auto')
    if (strcmp(shape,'circ') && city_radius > 2000) || (strcmp(shape,'poly') && strcmp(polytype,'closed'))
        bxmin=min(nonzeros(long))-deltax;
        bxmin=max(bxmin,0.5);
        bxmax=max(nonzeros(long))+deltax;
        bymin=min(nonzeros(lat))-deltay;
        bymax=max(nonzeros(lat))+deltay;
    else    
        if strcmp(bound,'auto')
            if strcmp(shape,'circ')
                bxmin=clong-deltax;
                bxmax=clong+deltax;
                bymin=clat-deltay;
                bymax=clat+deltay;
            else
                [~,xsize]=size(xint);
                vtemp=0;
                ixmax=1; 
                for ix=1:xsize
                    if kint(nt,ix) > 0 && vnet(nt,kint(nt,ix)) > vtemp
                        vtemp=vnet(nt,kint(nt,ix));
                        ixmax=ix;
                    end
                end    
                xm=xint(nt,ixmax);
                ym=yint(nt,ixmax);
                bxmin=xm-deltax;
                bxmin=max(bxmin,0.5);
                bxmax=xm+deltax;
                bymin=ym-deltay;
                bymax=ym+deltay;
            end
        end  
    end
end
%
x=(bxmin:dellatlongs:bxmax);
y=(bymin:dellatlongs: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+dellatlongs];
end    
%
h=zeros(sx,sy);
hx=zeros(sx,sy);
hy=zeros(sx,sy);
%
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,:), vstore(nt,:) );
    wq = pointwshortnq(lat,long,q900,v,rm,vse,rmse,utd,vtd,ush,vsh,y,x,h,hx,hy,timeres,wrad);
else
    w = pointwshortn(lat,long,v,rm,vse,rmse,utd,vtd,ush,vsh,y,x,h,hx,hy,timeres);
    wq=max(w-wrad,0).*q900;
    clear w
end    
temp=eprecip.*m_to_mm.*timeres.*3600.*rowa_over_rowl.*sum(wq,2);
rain(:,:)=temp;
z=rain;
clear temp
rain=transpose(rain);
%
if strcmpi(landfill,'y')
    rain(h'<1)=NaN;     % No contours over ocean
    rain(rain<10)=NaN;  % No contours where rainfall is less than 10 mm
end    
%
clf
m_proj(projection,'long',[bxmin bxmax], ...
   'lat',[bymin bymax]);
ld=landcolor;
hold on
if cres == 1
      m_coast('patch',ld,'edgecolor',shorecolor);
     % m_gshhs_l('color',transcolor);
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 pstates == 'y'
   m_states(bxmin,bxmax,bymin,bymax,'color',stcolor)
end
%
if strcmpi(wfill,'n') 
    [C1,h1]=m_contour(x,y,runitsfac*rain,runitsfac*raincntrs,'linewidth',1);
    colormap(wmap)
elseif strcmpi(wfill,'y') && strcmpi(sfill,'n')
    [C1,h1]=m_contour(x,y,runitsfac*rain,30,'linewidth',1);
    colormap(wmap)
    hp=m_pcolor(x,y,runitsfac*rain);
    a=colormap(wmap);
    for i=1:2
        a(i,:)=oceancolor;
    end
    colormap(a)
    set(hp,'EdgeAlpha',0,'FaceAlpha',wtrans);
    set(hp,'FaceColor','interp');
elseif strcmpi(wfill,'y') && strcmpi(sfill,'y')    
    [C1,h1]=m_contourf(x,y,runitsfac*rain,25);
    jeti=colormap('hsv');
    colormap(jeti(30:end,:));
end
m_grid('fontname',axisfont,'fontsize',axisfontsize,'fontweight',axisfontweight, ...
    'linestyle',gridline,'backgroundcolor',oceancolor)
if plottrack == 'y'
    m_plot(long(1:jmax),lat(1:jmax),'k','linewidth',3)
    %   This code places dots on track at 12 hour intervals
    %
    m_plot(long(1:6:jmax),lat(1:6:jmax),'.r','markersize',24)
    %
end    
if strcmpi(shape,'circ')
    m_plot(clong,clat,'.w','markersize',30)
end    
hold off
cb=colorbar;
set(cb,'fontname',axisfont,'fontweight','bold')
set(get(cb,'ylabel'),'String',runitslabel,'fontweight','bold','rotation',0,'horizontalalignment','left');
if isempty(storm)
    datea(1)=upper(datea(1));
    tit1=datea;
    if exist('yearstore','var') == 1 
        %if exist('tit3','var')
            %tit1=tit3;
        %else    
            tit1={tit1; ['Track number ' ntstring ', year ' char(num2str(yearstore(nt)))]};
        %end
    else
        tit1={tit1;['Track number ' ntstring]};
    end    
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')
%title('Accumulated Precipitation (mm)','fontweight','bold')
if exist('tmr7','var')
    stop(tmr7)
    delete(tmr7)
    clear tmr7
end
%
% 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 