function [vs2,dir,dayk,v,f,rp] = windpdfx( POI,cal )

% This function calculates time series and probability distributions of
% tropical cyclone wind speed at specified points of interest, using
% calibration factor cal. If POI = 'default' then clat and clong are used.
% If POI='poly' then the maximum wind along line segments is 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. 
%
% vs2 contains time series of wind speed (knots) at the the POI while dir
% is the direction from which the wind is blowing (degrees). These are n x
% k x s matrices 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.
%
% v is the range of wind speeds (knots) considered for histograms of
% maximum wind speeds at the POI, determined by the paramters vbeg,
% vend, and vres in the params.m file.
%
% f is the annual exceedence frequency, and rp is the return period
% (inverse of the annual probability) of storm peak winds 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 windspeeds in the histogram.
%
% Note that this function can be slow if many points of interest are used.
%
% Copyright WindRiskTech LLC, 2016
% September, 2016. Revised May, 2019.
%-------------------------------------------------------------------------
histtype='ks';    % Histogram type (kernal smoothed ('ks') or histcounts ('hc')
%
params    %  Load parameters
%vres=4;   %  Override value from params
%
load temp %#ok<LOAD>
fs3='sorted';
load(fs3) %#ok<LOAD>
%
if ischar(POI)
    if isnan(clat)
        disp('Point of interest must be specified 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    
%  
% 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;  
%
clear vtemp
if ischar(POI) && strcmp(POI,'poly')
    vtemp(:,1)=vmax';
    vs2=[];dir=[];dayk=[];
else    
    %
    yearbig=repmat(yearstore',[1,mrm]);
    datestore=datenum(double(yearbig),monthstore,daystore,hourstore,0,0); %#ok<DATNM>
    %
    [vs2, dir, dayk] = pointseriesn(latstore,longstore, ...
        vstore,rm,vsestore,rmse, ...
        datestore,uinc,vinc,plat,plong,timeres);
    % 
    %vs2=vs2.*cos(pi*(dir-130)/180);
    vmax2=max(vs2,[],2);
    vtemp(:,1,:)=vmax2(:,1,:);
end  
%
% Define parameters for an array vnew that contains all the peak winds in a
% given year
%
if max(size(yearstore)) > 1 && max(yearstore) > min(yearstore)
    qy=max(yearstore)-min(yearstore)+1;
    atemp=ismember((min(yearstore):max(yearstore)),yearstore);
    atemp=cast(atemp,'like',yearstore);
    yearset=nonzeros(atemp.*(min(yearstore):max(yearstore)));
    [~,my]=ismember(yearset,yearstore);
    delmy=my(2)-my(1);
else
    qy=1;
    my(1)=1;
    freqyear(1)=freq;
    delmy=nrm;
end
%
xi=(vbeg:vres:vend);
xib=xi;
xib(end+1)=vend+vres;
%
qx=max(size(xi));
f=zeros(sx,qx);
%
for k=1:sx      % Cycle through points of interest
    %
    % Statement below added October 2018 to account for storm wind peak
    % winds < xi
    %
    rat=sum(histcounts(vtemp(:,1,k),[vbeg, 300]))./nrm;
    %
    for i=1:qy  % Cycle through years in data set
        % 
        % Define vnew as the array containing the point peak wind speeds in each
        %    year i. This was added in October, 2016
        %
        if i < qy
            vnew=vtemp(my(i):my(i+1)-1,1,k);
        else
            vnew=vtemp(my(i):end,1,k);
        end    
        %
        if exist('ksdensity','file') > 0 && strcmpi(histtype,'ks')
            % Line below replaced with code below that, January 19 2018. Does
            % not make much difference. 
            %f=f-freqyear(i)*delmy*log(ksdensity(vnew,xi,'function','cdf'))/n;
            fm=-log(1-ksdensity(vnew,xi));
            f2=zeros(1,qx);
            for kk=1:qx
                f2(kk)=sum(fm(kk:qx));
            end    
            f(k,:)=f(k,:)+f2*freqyear(i)*cal*rat.*delmy/(nrm*f2(1));
        else
            ftt=histcounts(vnew,xib);
            fa=zeros(1,qx);
            for kj=1:qx
                fa(kj)=sum(ftt(1:kj));
            end
            fa=fa./fa(qx);
            f(k,:)=f(k,:)-delmy*freqyear(i).*cal*rat.*log(fa)/nrm;
        end       
    end
end
rp=1./(1-exp(-f));
v=xi;
%
end