%
%  This script compares best track to synthetic cumulative frequency
%  distributions 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 July, 2020
%--------------------------------------------------------------------------
%
vhisto  % This script is run to be sure synthetic data have been generated
if strcmp(Docdisp,'y')
    delete(hButton)
end    
clear alpha
%
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;
%u=transpose([y2;y]);
%
if strcmp(hbins,'ss')
    u=[y2;y]';
    h1=bar(xplot,u,1);
else    
    h1=plot(wunitsfac*xplot,y2,'.b','markersize',24);
    hold on
    h2=plot(wunitsfac*xplot,y,'r','linewidth',1.5);
    hold off
end
%
set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
%
if strcmp(hbins,'ss') 
    set(gca,'xticklabel',{'TS' '1' '2' '3' '4' '5'})
    xlabel('Saffir-Simpson Category','fontweight','bold')
else    
    clear u us us2
    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);
if strcmp(histf,'cum')
    ylabel('Annual Exceedence Frequency','fontweight','bold')
else
    ylabel('Annual Frequency Density','fontweight','bold')
end     
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','Northeast')
else    
    legend([h1,h2],tit1,tit2,'location','Northeast')
end    
%
if strcmp(ebar,'y')
    if strcmp(hbins,'ss')
        hold on
        y3=y-y3;
        y4=y4-y;
        j=errorbar(xplot,y,y3,y4);
        set(j,'linewidth',2,'color','g')
        yt=x*NaN;  %  This is to fool script into providing error line in legend
        jkk=plot(x,yt,'g','linewidth',2);
        hv=[h1,jkk];
        legend(hv,tit1,tit2,'90% of BT-size subsamples','location','NorthEast')
        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)];
            j=patch(wunitsfac*xtot,ytot,[0 0 1],'edgecolor','none');
            alpha(halpha)
        end    
        hv=[h1,h2,j];
        legend(hv,tit1,tit2,'90% of BT-size subsamples','location','NorthEast')
        hold off
        set(gca,'xlim',[min(wunitsfac*x)-wunitsfac*delx max(wunitsfac*x)+wunitsfac*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 