%
%  This script calculates and graphs a time series of 
%  rain rate, in mm/day, for a given storm at a specified location.
%
%       Copyright 2012 WindRiskTech L.L.C.
%       Revised December 2020 to remove cos(latitude) factor from hy and to
%       allow for higher resolution topography
%       Revised 1/2/2024 to change dayk to dayk1 on line 256
%--------------------------------------------------------------------------
%
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 
%
clf('reset')
if exist('shearstore','var') == 0
    shape='circ';
    clear T600store u850store v850store
    load('temp.mat')
end
fs3='sorted';
load(fs3)
%
if exist('nt','var') == 0
    if isempty(storm)
        prompt = {'Order by maximum at POI (1), lifetime maximum (2), or track number (3)?', ...
            '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});
        nttemp=str2double(answer{2});
    else    
       nswitch=2;
       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  
    %nttemp=input('Input track number');
    if nswitch == 1
       nt=nsortp(nttemp);
    elseif nswitch == 2   
       nt=nsort(nttemp);
    else
       nt=nttemp;
    end  
end
ntstring=num2str(nt);
%
%
if strcmpi(topores,'high')
    load bathymetry_high
else
    load bathymetry
end    
[ntopo,~]=size(bathy);
topores=360/ntopo;
toporesi=1/topores;
%
clear rain rtemp rainrate rainratec plat
%
pifac=acos(-1)/180;
dfac=60.*1.852;
knotfac=1852./3600;
sfac=1./(topores.*60.*1852);
%
%if strcmp(shape,'poly') == 1 || city_radius >= 2000
%    prompt = {'Enter POI latitude','Enter POI longitude','Enter POI name'};
%    dlg_title = 'Enter Point of Interest Data';
%    num_lines = 1;
%    answer = inputdlg(prompt,dlg_title,num_lines);
%    if isempty(answer) || isempty(answer{1})
%        return
%    end            
%    plat=str2double(answer{1});
%    plong=str2double(answer{2});
%    tmp=answer{3};
%    poi=cellstr(tmp);
%else
%    plat=clat;
%    plong=clong;
%    clear poi
%    datea(1)=upper(datea(1));
%    poi=cellstr(datea);
%end
iswitch='y';
if exist('vs','var')
    vs1=vs(nt,:);
end
if strcmp(shape,'circ') && city_radius < 2000
    plat=clat;
    plong=clong;
    clear poi
    datea(1)=upper(datea(1));
    poi=cellstr(datea);
    prompt = {'Use new point of interest (POI)? (y or n)'};
    dlg_title = 'New Point of Interest';
    num_lines = 1;
    answer = inputdlg(prompt,dlg_title,num_lines);
    iswitch=answer{1};
    if isempty(answer) || isempty(answer{1})
        iswitch='n';
    end        
end
%
[~,mm]=size(vstore);
if strcmp(iswitch,'y')
    prompt = {'Enter POI latitude:','Enter POI Longitude:'};
    dlg_title = 'Input for POI position';
    num_lines = 1;
    answer = inputdlg(prompt,dlg_title,num_lines);
    if isempty(answer) || isempty(answer{1})
        return
    end        
    plat=str2double(answer{1});
    plong=str2double(answer{2});
    poi=['Latitude= ',num2str(plat),', Longitude= ',num2str(plong)];
    if plong < 0
        plong=plong+360;
    end    
    %
    [~,mm]=size(vstore);
    [~,jmax]=min(vstore,[],2);
    jmax=jmax-1;
    ys=repmat(yearstore(nt),[1,mm]);
    datestore=datenum(double(ys),monthstore(nt,:),daystore(nt,:), hourstore(nt,:),0,0); %#ok<DATNM>
    [vs1, dir2, dayk1] = pointseriesn(latstore(nt,:),longstore(nt,:), ...
        vstore(nt,:),rmstore(nt,:),vsestore(nt,:),rmsestore(nt,:), ...
        datestore,uinc(nt,:),vinc(nt,:),plat,plong,timeres);
end    
%
if exist('vs1','var') == 0 || max(vs1) == 0
    msgbox('This storm does not affect point in question')
else
    %
    [n,m]=size(ut);
    utmod=ut(nt,:);
    vtmod=vt(nt,:);
    psize=max(size(plat));
    for i=1:psize
        if plong(i) < 0
            plong(i)=plong(i)+360;
        end
    end    
    %
    ush=zeros(1,m);
    vsh=ush;
    vdrift=1.5*3600/1852;
    vdrift=vdrift*latstore(nt,1)/(abs(latstore(nt,1))+1e-8);
    if exist('u850store','var') == 1
        ush=5.*knotfac.*(utmod-u850store(nt,:));
        vsh=5.*knotfac.*(vtmod-vdrift.*cos(pifac.*latstore(nt,:))-v850store(nt,:));
    end    
    lat=latstore(nt,:);
    long=longstore(nt,:);
    v=vstore(nt,:);
    vse=vsestore(nt,:);
    mon=monthstore(nt,:);
    day=daystore(nt,:);
    hour=hourstore(nt,:);
    year=yearstore(nt);
    %  
    % Scale and randomize radii of maximum wind
    %
    rfac=1+zeros(nt,m);
    if strcmp(randfac,'y')
        [nrm,mrm]=size(rmstore);
        rfac=outer(nrm,mrm);
    end  
    %
    rm=rmstore(nt,:).*rfac(nt,:).*magfac;
    rmse=rmsestore(nt,:).*rfac(nt,:).*magfac;
    %
    h=zeros(psize,1);
    hx=zeros(psize,1);
    hy=zeros(psize,1);
    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:psize
        ib=1+floor(toporesi.*plong(i));
        ibp=ib+1;
        if ibp > ntopo
            ibp=ibp-ntopo;
        end    
        ibs=1+floor(toporesi.*plong(i)-0.5);
        ibsp=ibs+1;
        plongs=plong(i);
        if ibs < 0
            ibs=ntopo;
            plongs=plong(i)+360;
        end    
        if ibsp > ntopo
            ibsp=ibsp-ntopo;
        end
        jb=1+floor(toporesi.*(plat(i)+90));
        jbs=1+floor(toporesi.*(plat(i)+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(i)+90)-jb;
        delx=1+toporesi.*plong(i)-ib;
        d1=(1.-delx).*(1.-dely);
        d2=dely.*(1.-delx);
        d3=delx.*(1.-dely);
        d4=delx.*dely;
        h(i)=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(i)+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)=(b1.*d1+b2.*d2+b3.*d3+b4.*d4)./cos(pifac.*plat(i));
        %
        b1=dhdy(ibs,jbs);
        b2=dhdy(ibs,jbs+1);
        b3=dhdy(ibsp,jbs);
        b4=dhdy(ibsp,jbs+1);
        hy(i)=(b1.*d1+b2.*d2+b3.*d3+b4.*d4);    
        %
    end  
    clear bathy
    if exist('T600store','var')
        clear q900 q600
        [ q900, q600 ] = qs900b( T600store(nt,:), vstore(nt,:) );
        ys=repmat(yearstore(nt),[1,mm]);
        datestore=datenum(double(ys),mon,day,hour,0,0); %#ok<DATNM>
        [wq,dayk1] = pointwnq(lat,long,q900,v,rm, ...
            vse,rmse,datestore,utmod,vtmod,ush,vsh,plat, ...
             plong,h,hx,hy,timeres,wrad);
    else
        [w,dayk1] = pointwn(lat,long,v,rm,vse,rmse,year,mon,day,hour,utmod,vtmod,ush, ...
            vsh,plat,plong,h,hx,hy,timeres);
        wq=max(w-wrad,0).*q900;
        clear w
    end   
    %
    rainrate=eprecip.*m_to_mm.*3600.*rowa_over_rowl.*wq;
    qr=max(size(rainrate));
    rain=zeros(1,qr);
    for i=1:qr
        rain(i)=timeres.*sum(rainrate(1:i));
    end    
    %
    datea(1)=upper(datea(1));
    for i=1:psize
        nplotr=2;
        while nplotr > 1
            nplotr=menu('Rain Probability Plots','Stop','Rainfall Rate', ...
                'Accumulated Rainfall');
            if nplotr == 1    
                nplotr=0;
                break
            elseif nplotr == 2    
                x=dayk1;
                [~,imax]=max(x);
                y=runitsfac*rainrate;
                h=plot(x(1:imax),y(1:imax));
                yvmax=10.*ceil(0.1.*max(y))+10;
                set(h,'linewidth',2)
                set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
                set(gca,'XLim',[x(1) max(x)])
                datetick('x','dd')
                ylabl=char(strcat('Rain rate (',runitslabel,'/hr)'));
                ylabel(ylabl)
                %
                month1=datestr(x(1),'mmmm');
                month2=datestr(x(imax),'mmmm');
                day1=datestr(x(1),'dd');
                day2=datestr(x(imax),'dd');
                line2=[char(poi) char(', Track Number ') char(ntstring)];
                titstring2={datea;line2};
                title(line2,'interpreter','none')
                if exist('yearstore','var')
                    if strcmp(month1,month2)
                        titstring3=[month1 ' ' num2str(yearstore(nt))];
                    else
                        titstring3=[month1 ' - ' month2 ' ' num2str(yearstore(nt))];
                    end
                else
                    if strcmp(month1,month2)
                        titstring3=month1;
                    else
                        titstring3=[month1 ' - ' month2];
                    end
                end    
                xlabel(titstring3)
                %
            elseif nplotr == 3    
                %
                x2=dayk1;
                [~,imax]=max(x2);
                y2=runitsfac*rain;
                h=plot(x2(1:imax),y2(1:imax));
                yvmax=10.*ceil(0.1.*max(y2))+10;
                set(h,'linewidth',2)
                set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
                %set(gca,'XLim',[x2(1) x2(end)]) 
                set(gca,'XLim',[x2(1) max(x2)])
                datetick('x','dd')
                ylabl=char(strcat('Accumulated Rainfall (',runitslabel,')'));
                ylabel(ylabl)
                %
                month1=datestr(x2(1),'mmmm');
                month2=datestr(x2(imax),'mmmm');
                day1=datestr(x2(1),'dd');
                day2=datestr(x2(imax),'dd');
                line2=[char(poi) char(', Track Number ') char(ntstring)];
                %title(titstring2,'interpreter','none')
                title(line2,'interpreter','none')
                if exist('yearstore','var')
                    if strcmp(month1,month2)
                        titstring3=[month1 ' ' num2str(yearstore(nt))];
                    else
                        titstring3=[month1 ' - ' month2 ' ' num2str(yearstore(nt))];
                    end
                else
                    if strcmp(month1,month2)
                        titstring3=month1;
                    else
                        titstring3=[month1 ' - ' month2];
                    end
                end    
                xlabel(titstring3)
            end
            %
            % Add button to access documentation of script
            %
            if strcmp(Docdisp,'y')
                try
                   delete(hButton)
                catch
                end    
                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            
        end     
    end    
    if exist('tmr7','var')
        stop(tmr7)
        delete(tmr7)
        clear tmr7
    end
end
%
% Clean up
%
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 