%
%  This script compares best track to synthetic return periods by intensity.
%  This version plots bounds within which lie 90% of subsamples of the 
%  synthetic tracks, randomly drawn from a Poisson distribution based on the 
%  total filtered best track count.
%  January 2009  
%  Copyright 2009 WindRiskTech L.L.C.
%  Last modified December, 2020
%--------------------------------------------------------------------------
%
return_period  % This script is run to be sure synthetic data have been generated
if strcmp(Docdisp,'y')
    delete(hButton)
end    
%
if exist('shearstore','var') == 0
    shape='circ';  %#ok<NASGU> %Default shape, used for older event sets
    fs='temp.mat';
    load(fs, 'bas','shape','clat','clong','city_radius','vbeg')
end
%
bestproc   %  Process best tracks
%
if (strcmp(shape,'poly') || city_radius >= 2000) && htype == 3
    htype=2;
end    
%
if htype == 1
    vtemp=vsmmax;
elseif htype == 2
    vtemp=vmaxb;
else
    vtemp=vsmaxb;
end  
u=ismember(round(vtemp),[round(vbeg):300]); %#ok<NBRAK>
vtemp=vtemp.*u;
%
[vhistbest,~]=histcounts(vtemp,xb);  % Replaced histc with histcounts, May, 2017
y2=zeros(1,histsize);
sum5=zeros(1,histsize);
sum95=zeros(1,histsize);
qb=max(size(bestyears));
pswitch=1;
try
    s=poissrnd(2);
catch err
    pswitch=0;
    msgbox('poissrnd.m missing from MATLAB path; using free (but slower) function Poissrnd instead')
end    
if pswitch == 1   
    for i=1:histsize
        if strcmp(histf,'cum')
            y2(i)=sum(vhistbest(1,i:end));
        else     
            y2(i)=vhistbest(1,i);
        end    
        ab=y(i).*qb;
        sum5(i)=poissinv(0.05,ab);
        sum95(i)=poissinv(0.95,ab);        
    end    
else
    for i=1:histsize
        if strcmp(histf,'cum')
            y2(i)=sum(vhistbest(1,i:end));
        else     
            y2(i)=vhistbest(1,i);
        end  
        ab=y(i).*qb;
        ab=max(ab,0.00001);
        dum=zeros(1,1000);
        for j=1:1000
            dum(j)=Poissrnd(double(ab));
        end
        dum2=sort(dum);
        sum5(i)=dum2(50);
        sum95(i)=dum2(950);
    end     
end    
%
nbc=sum(vhistbest);
denom=1/qb;
y2=y2.*denom;
y3=sum5.*denom;
y4=sum95.*denom;
%
if strcmp(hbins,'ss')
    u=[1./(1-exp(-y2));1./(1-exp(-y))]';
    h1=bar(xplot,u);
    set(gca,'yscale','log')
    ylimits=get(gca,'ylim');
    ylimits(1)=1;
    ylimits(2)=10^(0.1*ceil(10.*log10(ylimits(2))));
    set(gca,'ylim',ylimits)
    set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
    set(gca,'xticklabel',{'TS','1' '2' '3' '4' '5'})
    xlabel('Saffir-Simpson Category','fontweight','bold')
else    
    h1=semilogy(wunitsfac*xplot,1./(1-exp(-y2)),'.b','markersize',24);
    hold on
    h2=semilogy(wunitsfac*xplot,1./(1-exp(-y)),'r','linewidth',1.5);
    hold off
    set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
    set(gca,'xlim',wunitsfac*[min(x)-10 max(x)+10])
    xlabl=char(strcat('Maximum Wind Speed (',wunitslabel,')'));
    xlabel(xlabl,'fontweight','bold')
end
%
yearbeg=min(bestyears);
yearend=max(bestyears);
ylabel('Return Period (years)','fontweight','bold')
tit1=[num2str(nbc) ' best tracks, '  num2str(yearbeg) ' to ' ...
    num2str(yearend)];
tit2=[num2str(floor(sumtotal)) ' synthetic tracks'];
if strcmp(hbins,'ss')
    legend(h1,tit1,tit2,'location','best','AutoUpdate','off')
else    
    legend([h1,h2],tit1,tit2,'location','best','AutoUpdate','off')
end    
%
ymax=get(gca,'ylim');
if strcmp(ebar,'y')
    if strcmp(hbins,'ss')
        hold on
        ycentral=1./(1-exp(-y));
        ylower=ycentral-max(1./(1-exp(-y4)),1);
        yupper=min(1./(1-exp(-y3)),ymax(2))-ycentral;
        j=errorbar(xplot,ycentral,ylower,yupper);
        set(j,'linewidth',2,'color','g')
        yt=x*NaN;  %  This is to fool script into providing error line in legend
        jkk=semilogy(x,1./(1-exp(-yt)),'g','linewidth',2);
        hv=[h1,jkk];
        legend(hv,tit1,tit2,'90% of BT-size subsamples','location','best','AutoUpdate','off')
        set(j,'linestyle','none')
        hold off
    else
        hold on
        bmax=max(size(x))-1;
        for kd=1:bmax
            delx=0.5*(xplot(2)-xplot(1));
            xtot=[xplot(kd) xplot(kd+1) xplot(kd+1) xplot(kd)];
            ytot=[y3(kd) y3(kd+1) y4(kd+1) y4(kd)];
            ytot=(1./(1-exp(-ytot)));
            ytot=min(ytot,ymax(2));
            j=patch(wunitsfac*xtot,ytot,[0 0 1],'edgecolor','none');
            alpha(halpha)
            if strcmp(rgridx,'y')
                set(gca,'xgrid','on')
            end
            if strcmp(rgridy,'y')
                 set(gca,'ygrid','on')
            end  
        end    
        hv=[h1,h2,j];
        legend(hv,tit1,tit2,'90% of BT-size subsamples','location','best','AutoUpdate','off')
        hold off
        set(gca,'xlim',wunitsfac*[min(x)-delx max(x)+delx])
    end
end 
%
if isempty(storm)
    tit1=datea;
else
    tit1=[storm(1:2) ' ' storm(3:4) ' ' datea(1:4) ' ' datea(5:6) ' ' datea(7:8)...
    ' ' datea(9:10) ' GMT'];
end    
if strcmp(shape,'circ') 
    if htype == 2
        tit2=['All events within ' num2str(city_radius) ' km of '];
        fsize=11;
    elseif htype == 3
        tit2='All events at POI, ';
        fsize=11;
    else
        tit2='Storm Lifetime Maximum Wind, All Events in ';
        fsize=11;
    end    
else
    if strcmp(polytype,'closed')
        tit2='All events within closed polygon, event set ';
        fsize=11;
    else    
        tit2='All events at points of intersection with line segments in poly.in, event set ';
        fsize=10;   
    end
end
titnet=[tit2 tit1];
title(titnet,'fontsize',fsize,'fontweight','bold','interpreter','none')
%
% 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 