function [rainrate,rain,dayk,rx,f,rp] = rainpdfx( POI,cal,rmin,rmax )
%
%  This function calculates pdfs of tropical cyclone rain rate and
%  accumulated rainfall at specified points of interest, using calibration
%  factor cal. If POI = 'default' then clat and clong are used. Otherwise,
%  POI(:,1) is a vector of latitudes, and POI(:,2) is a vector of
%  longitudes of points of interest at which times series of wind speed and
%  direction are calculated.
%
%  This function operates on the active event set, produced by running
%  prep.m.
%
%  cal is a calibration factor. cal=1 leaves the dataset unaltered. 
%
%  rmin and rmax are the minimum and maximun storm total rainfall (mm) to
%  be considered in constructing the storm total rainfall return periods.
%  Note that rmin itself is forced to be larger than a critical value set
%  below.
%
%  rainrate is a matrix containing the rain rates in mm/hr at the POI. It
%  is dimensioned n x k x s, where n is the number of storms, k is the
%  length of the time series, which is set internally, and s is the number
%  of points of interest. (Note that if s=1, n x k matrices are returned.)
%  The matlab datenumbers corresponding to these series are contained in
%  dayk.
%
%  rain is a vector of dimension n x s containing the storm total rainfall
%  in mm at the POI.
%
%  rx is the range of storm total rainfall (mm) considered for histograms
%  of rainfall at the POI, determined by the paramters rmin and rmax
%  described above.
%
% f is the annual exceedence frequency, and rp is the return period
% (inverse of the annual probability) of storm total rainfall at the POI.
% These are s x m matrices, where s is the number of points of interest and
% m is the size of the range of rain rates in the histogram.
%
%  Note that this function can be slow if many points of interest are used.
%      
%  Copyright WindRiskTech LLC, 2016
%       Revised December 2020 to remove cos(latitude) factor from hy and to
%       allow for higher resolution topography
%
%-------------------------------------------------------------------------
%  It is very important to put lower bounds on rain for data sets in which
%  the number of storms causing rain at the POI is a fraction of the total
%  number of storms. No harm is done by having such bounds when that
%  fraction is one or close to one. 
%
%-------------------------------------------------------------------------
rminlower=20;     % Minimum storm rainfall (mm) to consider
%------------------------------------------------------------------------
%
rmin=max(rmin,rminlower);
params    %  Load parameters
load temp %#ok<LOAD>
load sorted %#ok<LOAD>
%
if ischar(POI)
    if isnan(clat)
        disp('Must specify points of interest for this event set')
        return
    end    
    plat=clat;
    plong=clong;
else
    plat=POI(:,1)';
    plong=POI(:,2)';
end
%
sx=max(size(plong));
for i=1:sx
    if plong(i) < 0
        plong(i)=plong(i)+360;
    end
end   
%
%
if strcmpi(topores,'high') %#ok<NODEF>
    load bathymetry_high bathy
else
    load bathymetry bathy
end    
[ntopo,~]=size(bathy);
topores=360/ntopo;
toporesi=1/topores;
%
clear rain rtemp rainrate rainratec
sfac=1./(topores.*60.*1852);
pifac=acos(-1)./180;
%
knotfac=1852./3600;
%
[n,m]=size(ut);
%
ush=zeros(n,m);
vsh=ush;
vdrift=1.5.*3600./1852;
vdrift=vdrift.*latstore(1,1)./(abs(latstore(1,1))+1e-8);
if exist('u850store','var') == 1
    ush=5.*knotfac.*(ut-u850store);
    vsh=5.*knotfac.*(vt-vdrift.*cos(pifac.*latstore)-v850store);
end    
%  
% Scale and randomize radii of maximum wind
%
[nrm,mrm]=size(rmstore);
rfac=magfac.*(1+zeros(nrm,mrm));
if strcmp(randfac,'y')
    rfac=magfac.*outer(nrm,mrm);
end  
%
rm=rmstore.*rfac;
rmse=rmsestore.*rfac;
%
% Define parameters for an array vnew that contains all the peak winds in a
% given year
%
if max(yearstore) == min(yearstore)
    freqyear(1)=freq;
end   
%
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));
%
h=zeros(1,sx);
hx=zeros(1,sx);
hy=zeros(1,sx);
%
for i=1:sx
    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
ys=repmat(yearstore',[1,mrm]);
if exist('T600store','var')
    clear q900 q600
    [ q900, q600 ] = qs900b( T600store,vstore ); %#ok<ASGLU>   
    datestore=datenum(double(ys),monthstore,daystore,hourstore,0,0); %#ok<DATNM>
    [wq,dayk] = pointwnq(latstore,longstore,q900,vstore,rm, ...
       vsestore,rmse,datestore,ut,vt,ush,vsh,plat, ...
       plong,h,hx,hy,timeres,wrad);
else
    [w,dayk] = pointwn(latstore,longstore,vstore,rmstore, ...
       vsestore,rmsestore,ys,monthstore,daystore,hourstore,ut,vt,ush,vsh,plat, ...
       plong,h,hx,hy,timeres);
    wq=max(w-wrad,0).*q900; %#ok<NODEF>
    clear w
end   
%
wpfac=1;
%{
if wprofile == 2
    wpfac=2;
end
%}    
rainrate=wpfac.*eprecip.*m_to_mm.*3600.*rowa_over_rowl.*wq;
rainrate=real(rainrate);
rain(:,:)=timeres.*sum(rainrate,2);
%
%xmax=10*ceil(0.1*max(rain));
xmax=rmax;
xmax=max(xmax,10);
%
%maxrain=10*ceil(0.1*max(rain));
%
xint=0.02*xmax;
xi=(rmin:xint:xmax);         % Changed 12/11/13 and 10/14/2018
qx=max(size(xi));
xil=xi;
xil(end+1)=100000;
xil(1)=max(xil(1),0.01); % This is essential to exclude masked events
%
f=zeros(sx,qx);
year=min(yearstore):max(yearstore);
qy=max(size(year));
for j=1:sx
    rat=sum(histcounts(rain(:,j),[rmin, 100000]))./n;
    for i=1:qy
        % 
        % Define rainnew as the array containing the point accumulated rainfall in each
        %    year i. This was added in October, 2016
        %
        mask=ismember(yearstore,year(i));
        rainnew=mask'.*rain(:,j);
        u=histcounts(rainnew,xil);
        fd=freqyear(i)*rat*cal*u./(0.001+sum(u));
        fdc=zeros(1,qx);
        for k=1:qx
            fdc(k)=sum(fd(k:end));
        end
        f(j,:)=f(j,:)+fdc/qy; 
    end    
end    
%f=max(f,1e-6);
%  Smooth pdf, if desired
%{
fs=f;
fs(3:end-2)=smooth134(f);
fs(3:end-2)=smooth134(fs);
fs(3:end-2)=smooth134(fs);
f=fs;         
%}
rp=1./(1-exp(-f));    % Changed 12/10/2013
rx=xi;
%