%
%  This script calculates key statistics of secondary (SEs) and primary
%  (PEs) eyewall replcements in the synthetic events. The first plot shows
%  the probability of encountering a particular number of SEs. The second
%  plot shows the means and standard deviations of the wind speeds
%  associated with particular numbers of SEs, and the third plot shows the
%  means and standard deviations of the power dissipation over each event,
%  associated with particular numbers of SEs.  The next three plots show
%  the same quantities for primary eyewall replacement cycles.
%  Copyright WindRiskTech, L.L.C.
%  Last modified June, 2026
%--------------------------------------------------------------------------
secrit=5;  % Critical sample size, in each bin of secondary eyewall counts, 
%            below which data are not plotted
%--------------------------------------------------------------------------
warning('off','MATLAB:load:variableNotFound')
load temp rmsestore vsestore vnetse vnet
if exist('vsestore','var') == 0 || sum(nonzeros(vsestore)) == 0
    errordlg('No secondary eyewalls in this event set');
    warning('on','MATLAB:load:variableNotFound')
    return
end    
params
clf('reset')
[n,m]=size(vnet);
%
f=5e-5;
[~,jmax]=min(vstore,[],2);
jmax=jmax-1;
jmax=max(jmax,2);
vmax=max(vnet,[],2);
pdi=sum(vnet.^3,2);
vms=vstore.*1850/3600;
M=sqrt(1000*rmstore.*vms+0.5*f*1e6*rmstore.^2);
senum=zeros(1,n);
penum=zeros(1,n);
for i=1:n
    delM=M(i,2:jmax(i))-M(i,1:jmax(i)-1);
    senum(i)=sum(ismember(round(delM),50:1000));
    penum(i)=sum(ismember(round(delM),-1000:-50));
end    
semax=max(senum);
pemax=max(penum);
%
x=0:semax;
xe=x-0.5;
xe(end+1)=xe(end)+0.5;
qm=max(size(x));
%
xp=0:pemax;
xep=xp-0.5;
xep(end+1)=xep(end)+0.5;
qmp=max(size(xp));
%
[us,~] = histcounts(senum,xe);
x(us<secrit)=NaN;  % Do not graph stats for sample size less than secrit
us(us<secrit)=NaN;
[up,~] = histcounts(penum,xep);
xp(up<secrit)=NaN;  % Do not graph stats for sample size less than secrit
up(up<secrit)=NaN;
%
nentry=2;
while nentry > 1
    nentry=menu('Secondary and Primary Eyewall Replacement Stats','Stop', ...
        'Secondary Eyewall Replacement Probability', ...
        'Secondary Eyewall Replacement vs Wind Speed', ...
        'Secondary Eyewall Replacement vs Power Dissipation', ...
        'Primary Eyewall Replacement Probability', ...
        'Primary Eyewall Replacement vs Wind Speed', ...
        'Primary Eyewall Replacement vs Power Dissipation');
    if nentry == 1
        break
    elseif nentry == 2
        y=us./n;
        bar(x,y,0.8)
        set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        xlabel('Number of Secondary Replacement Cycles')
        ylabel('Probability')
        title('Probability of Secondary Eyewall Replacements')
        %
    elseif nentry == 3
        vse=zeros(1,qm);
        vsestd=zeros(1,qm);
        for i=1:qm
            vse(i)=mean(vmax(senum'== i-1));
            vsestd(i)=std(vmax(senum'== i-1),1);
        end    
        %
        bar(x,vse,0.8)
        hold on
        errorbar(x,vse,min(vsestd,vse,"includemissing"),vsestd, ...
            'linestyle','none','linewidth',2,'color','g')
        hold off
        set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        xlabel('Number of Secondary Replacement Cycles')
        ylabel('Mean Lifetime Maximum Wind Speed (knots)')
        title('# Secondary Cycles vs. Mean Lifetime Maximum Wind Speed')
        y=vse;
        %
    elseif nentry == 4     
        %
        pdise=zeros(1,qm);
        pdisestd=zeros(1,qm);
        for i=1:qm
            pdise(i)=mean(pdi(senum'== i-1));
            pdisestd(i)=std(pdi(senum'== i-1),1);
        end    
        bar(x,pdise,0.8)
        hold on
        errorbar(x,pdise,min(pdisestd,pdise,"includemissing"),pdisestd, ...
            'linestyle','none','linewidth',2,'color','g')
        hold off
        set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        xlabel('Number of Secondary Replacement Cycles')
        ylabel('Mean Power Dissipation Index')
        title('# Secondary Cycles vs. Power Dissipation')
        y=pdise;
        %
    elseif nentry == 5
        y=up./n;
        bar(xp,y,0.8)
        set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        xlabel('Number of Secondary Replacement Cycles')
        ylabel('Probability')
        title('Probability of Primary Eyewall Replacements')
        %
    elseif nentry == 6
        vpe=zeros(1,qmp);
        vpestd=zeros(1,qmp);
        for i=1:qmp
            vpe(i)=mean(vmax(penum'== i-1));
            vpestd(i)=std(vmax(penum'== i-1),1);
        end    
        bar(xp,vpe,0.8)
        hold on
        errorbar(xp,vpe,min(vpestd,vpe,"includemissing"),vpestd, ...
            'linestyle','none','linewidth',2,'color','g')
        hold off
        set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        xlabel('Number of Primary Replacement Cycles')
        ylabel('Mean Lifetime Maximum Wind Speed (knots)')
        title('# Primary Cycles vs. Mean Lifetime Maximum Wind Speed')
        y=vpe;
        %
    elseif nentry == 7     
        %
        pdipe=zeros(1,qmp);
        pdipestd=zeros(1,qmp);
        for i=1:qmp
            pdipe(i)=mean(pdi(penum'== i-1));
            pdipestd(i)=std(pdi(penum'== i-1),1);
        end 
        bar(xp,pdipe,0.8)
        hold on
        errorbar(xp,pdipe,min(pdipestd,pdipe,"includemissing"),pdipestd, ...
            'linestyle','none','linewidth',2,'color','g')
        hold off
        set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        xlabel('Number of Primary Replacement Cycles')
        ylabel('Mean Power Dissipation Index')
        title('# Primary Cycles vs. Power Dissipation')
        y=pdipe;
    end    
    %
    % Add button to access documentation of script
    %
    if strcmp(Docdisp,'y')
        try
            delete(hButton)
        catch
        end
        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    
end    
%
% Clean up
%
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 