%
%  This script calculates and graphs pdfs of tropical cyclone
%  wind speed at specified points of interest.
%       December 2011
%       Copyright 2011 WindRiskTech L.L.C.
%       Revised October, 2018
%       Last revised September 2020
%-------------------------------------------------------------------------
%
histtype='hc';    % Histogram type (kernal smoothed ('ks') or histcounts ('hc')
smoothing='y';    % Simple smoothing of return periods
%
clear poi plat plong vtemp fa
clf('reset')
%
params    %  Load parameters
pfield='wind';
%
if strcmp(randfac,'y')
    %disp('Caution: Random draw for wind radii turned on; results will differ with each run')
    if exist('tmr5','var')
        stop(tmr5)
        delete(tmr5)
        clear tmr5
    end   
    tmr5=timer;
    set(tmr5,'startdelay',2);
    set(tmr5,'timerfcn','try; delete(hms); catch err; end')
    hms=msgbox('Caution: Random draw for wind radii turned on; results will differ with each call');
    start(tmr5)
end  
%
if exist('shearstore','var') == 0
    shape='circ';
    load('temp.mat')
end
fs3='sorted';
load(fs3)
load temp filtyears
sfac=1./(0.25.*60.*1852);
pifac=acos(-1)./180;
%
if exist('plat','var') == 0
    if strcmp(shape,'poly') == 1 || city_radius >= 2000
        %plat=input('Enter POI latitude');
        %plong=input('Enter POI longitude');
        %tmp=input('Enter POI name','s');
        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
        prompt = {'Use default POI (y or n; default is y)?'};
        dlg_title = 'Use default POI?';
        num_lines = 1;
        answer = inputdlg(prompt,dlg_title,num_lines);
        pswitch=answer{1};
        %pswitch=input('Use default POI (y or n; default is y)?','s');
        if isempty(pswitch)
            pswitch='y';
        end    
        if pswitch == 'n'
            %plat=input('Enter POI latitude');
            %plong=input('Enter POI longitude');
            %tmp=input('Enter POI name','s');
            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);
            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
    end
end    
%
knotfac=1852./3600;
%
if plong < 0
    plong=plong+360;
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;
%
[~,jmaxd]=min(vstore,[],2);
jmaxd=jmaxd-1;
[n,m]=size(ut);
%
[vs2, ~] = pointshortn(latstore,longstore,vstore,rm,vsestore,rmse, ...
    uinc,vinc,plat,plong,timeres);
vmax2=max(vs2,[],2);
%
% 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)
    qyb=max(yearstore)-min(yearstore)+1;
    yearset=filtyears;
    [~,my]=ismember(yearset,yearstore);
    my(end+1)=n;
    vnew=zeros(1,my(2)-1); 
    delmy=my(2:end)-my(1:end-1);
    qy=max(size(yearset));
else
    qy=1;
    my=1;
    my(2)=n+1;
    freqyear(1)=freq;
    delmy=n;
end    
%
tit1='Peak Wind Speed at ';
%
vtemp(:,1)=vmax2(:,1,1);
xmax=10*ceil(0.1*max(vtemp));
xmax=max(xmax,vbeg+40);
xint=0.01*xmax;
xi=(xint:xint:xmax);
xib=xi;
xib(end+1)=xmax+xint; 
%
% Statement below added October 2018 to account for storm wind peak winds <
% xi
%
rat=sum(histcounts(vtemp,[xint, 300]))./nrm;
%
qx=max(size(xi));
f=zeros(1,qx);
%
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
    %
        vnew=vtemp(my(i):my(i+1)-1);
    %
    if exist('ksdensity','file') && 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=f+f2*freqyear(i)*rat.*delmy(i)/(n*f2(1));
    else
        ftt=histcounts(vnew,xib);
        fa=zeros(1,qx);
        for kj=1:qx
            %fa(kj)=sum(ftt(1:kj));
            fa(kj)=sum(ftt(kj:qx));
        end
        %fa=fa./fa(qx);
        %f=f-delmy(i)*freqyear(i).*rat.*log(fa)/n;
        f=f+fa*freqyear(i)*rat.*delmy(i)/(n*fa(1));
    end    
end   
[~,nbeg]=min(max(vbeg-xi,0));  % Find index value of xi corresponding to wind speed vbeg
rp=1./(1-exp(-f(nbeg:end)));
if strcmpi(smoothing,'y')
    rp2=rp;
    rp2(3:end-2)=smooth134(rp);
    rp2(3:end-2)=smooth134(rp2);
    rp=rp2;
end
h=semilogy(wunitsfac*xi(nbeg:end),rp,'.b');     
y=rp;
x=wunitsfac*xi(nbeg:end);
set(h,'markersize',12)
set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
ylabel('Return Period (Years)','fontweight','bold')
xlabl=char(strcat('Peak storm wind speed (',wunitslabel,')'));
xlabel(xlabl,'fontweight','bold')
if strcmp(rgridx,'y')
     set(gca,'xgrid','on')
end
if strcmp(rgridy,'y')
     set(gca,'ygrid','on')
end  
clear titnet
titnet=[tit1 char(poi(1))];
title(titnet,'fontweight','bold','fontsize',11,'interpreter','none');
xadj=x;
yadj=y;
%
% 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 