function [rain,rainrate,dayk] = raingen(plat,plong)
%
%  This script calculates accumulated rain and rainrates at specified
%  latitudes (plat) and longitudes (plong) for the active event set. The
%  rainrates are a function of time with the time given by the matlab
%  datenumber array dayk. plat and plong may be scalars or 1-D vectors
%  containing a series of lat-long points. 
%       October 2018
%       Copyright 2018 WindRiskTech L.L.C.
%       Revised December 2020 to remove cos(latitude) factor from hy and to
%       allow for higher resolution topography
%
%-------------------------------------------------------------------------
%
params    %  Load parameters
load temp %  Load active event set
%
if strcmp(randfac,'y')
    %disp('Caution: Random draw for wind radii turned on; results will differ with each call')
    if exist('tmr7','var')
        stop(tmr7) %#ok<NODEF>
        delete(tmr7)
        clear tmr7
    end   
    tmr7=timer;
    set(tmr7,'startdelay',6);
    set(tmr7,'timerfcn','try; delete(hms); catch err; end')
    [~]=msgbox('Caution: Random draw for wind radii turned on; results will differ with each call');
    start(tmr7)
end 
%
sx=max(size(plat));
fs3='sorted';
load(fs3) %#ok<*LOAD>
%
if strcmpi(topores,'high') %#ok<NODEF>
    load bathymetry_high
else
    load bathymetry
end    
[ntopo,~]=size(bathy); %#ok<NODEF>
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);
if min(plong) < 0
    plong=plong+360;
end  
%
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    
lat=latstore;
long=longstore;
v=vstore;
vse=vsestore;
%  
% 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;
%
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
years=repmat(yearstore',[1,m]);
datestore=datenum(double(years),monthstore,daystore,hourstore,0,0);
clear q900 q600
[ q900, ~ ] = qs900b( T600store, vstore );
[wq,dayk] = pointwshortnqdx(lat,long,q900,v,rm,datestore,vse,rmse,ut,vt,ush,vsh,plat,plong,h,hx,hy,timeres,wrad);
%
rainrate=eprecip.*m_to_mm.*3600.*rowa_over_rowl.*wq;
clear rain
rain(:,:)=timeres.*sum(rainrate,2);