%
% This script compares the annual cycle of synthetic storm tracks to best
% track data. The synthetic track cycle can be normalized so that the 
% synthetic annual total is the same as for the best track data. This 
% version includes 90% confidence limits based on the assumption that the
% best track data for each month have a Poisson distribution.
%   Copyright 2009 WindRiskTech L.L.C.
%   Last modified April, 2020 to include storms just at the critical wind
%   speeds
%-------------------------------------------------------------------------
%
annual % This script is run to be sure synthetic data have been generated
%
bestproc   %  Process best tracks
%
clear alpha
[m,n]=size(vbest);
[~,jmax]=min(vbest,[],2);
jmax=jmax-1;
%
[~,jmin]=min(max((startv-vbest),0),[],2);
%
yearbeg=min(bestyears);
yearend=max(bestyears);
%
y2=zeros(1,12);
x=1:12;
qp=max(size(bestyears));
[~,monthbest,~,~,~,~]=datevec(datebest);
monthtempb=zeros(1,m);
for i=1:m
    monthtempb(i)=monthbest(i,jmin(i));
end
vmask=ceil(min(max(vmaxb+1-vcritann,0),1)); % Added 1, April, 2020, to include events at the threshold intensity
for k=1:12
    y2(k)=sum(vmask.*ismember(monthtempb,k));
end
%
dlow=zeros(1,12);
dhigh=zeros(1,12);
%
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:12
        ab=y(i).*qp;
        dlow(i)=poissinv(0.05,ab);
        dhigh(i)=poissinv(0.95,ab);
    end
else
    for i=1:12
        ab=y(i).*qp;
        ab=max(ab,0.00001);
        dum=zeros(1,1000);
        for j=1:1000
            dum(j)=Poissrnd(double(ab));
        end
        dum2=sort(dum);
        dlow(i)=dum2(50);
        dhigh(i)=dum2(950);
    end
end    
denom=1./qp;
%
y2=y2.*denom;
%
y3=dlow.*denom;
y4=dhigh.*denom;
%
% Normalize, if desired
%
acal=sum(y2)./sum(y);
ynorm=y.*acal;
if strcmp(normalize,'y')
    y3=y3.*sum(y2)./sum(y);
    y4=y4.*sum(y2)./sum(y);
    y=ynorm;
end
%
error=sum(abs(y2-ynorm));
errornorm=100*error./sum(y2);
aberr=['mean absolute value of difference= ' num2str(error) ' storms per month'];
aberrnorm=['Mean difference= ' num2str(errornorm) ' %'];
%
h2=plot(x,y,'r','linewidth',2);
hold on
h1=plot(x,y2,'.b','markersize',24);
%h3=plot(x,y2,'b','linewidth',1.5);
hold off
set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
set(gca,'xlim',[0.5 12.5])
xlabel('Month','fontweight','bold')
ylabel('Number of storms per month','fontweight','bold')
tit1=[num2str(sum(vmask)) ' best tracks, '  num2str(yearbeg) ' to ' ...
    num2str(yearend)];
tit2=[num2str(msyn) ' synthetic tracks'];
legend([h1,h2],tit1,tit2,'location','best')
if strcmp(ebar,'y')
    hold on
    bmax=max(size(x));
    for kd=1:bmax-1
        delx=0.5*(x(2)-x(1));
        xtot=[x(kd) x(kd+1) x(kd+1) x(kd)];
        ytot=[y3(kd) y3(kd+1) y4(kd+1) y4(kd)];
        j=patch(xtot,ytot,[0 0 1],'edgecolor','none');
        alpha(halpha)
    end
    legend([h1,h2,j],tit1,tit2,'90% of BT-size subsamples','location','best')
    hold off
end
if isempty(storm)
    datea(1)=upper(datea(1));
    tit1=datea;
else
    tit1=[storm(1:2) ' ' storm(3:4) ' ' datea(1:4) ' ' datea(5:6) ' ' datea(7:8)...
    ' ' datea(9:10) ' GMT'];
end    
title(tit1,'fontsize',12,'fontweight','bold','interpreter','none')
%
% Add button to access documentation of script
%
if strcmp(Docdisp,'y')
    delete(hButton)
    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 
%hms=msgbox(aberrnorm);