%
%  This script calculates key statistics of secondary eyewalls (SEs) 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, nad the third plot shows the means and standard deviations of the
%  integrals of wind speed over the life of each event, associated with
%  particular numbers of SEs.  Copyright WindRiskTech, L.L.C.
%  Last modified June, 2017
%--------------------------------------------------------------------------
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;
fac2=1/3.8; % In fortran code, rmstore and rmsestore are increased by 3.8.
M=sqrt(fac2*1000*1852*rmstore.*vstore/3600+0.5*f*1e6*fac2.^2*rmstore.^2); 
delM=M(:,2:end)-M(:,1:end-1);
delMp=max(min(round(1e-1*delM),1),0);
delMm=-min(max(round(1e-1*delM),-1),0);
senum=sum(delMp,2);
penum=sum(delMm,2)-1;
semax=max(senum);
pemax=max(penum);
%
x=0:12;
xe=x-0.5;
xe(end+1)=xe(end)+0.5;
qm=max(size(x));
%
[us,~] = histcounts(senum,xe);
facs=zeros(1,qm)+1;
facs(us<secrit)=0;  % Do not graph stats for sample size less than secrit
[up,~] = histcounts(penum,xe);
facp=zeros(1,qm)+1;
facp(up<secrit)=0;  % Do not graph stats for sample size less than secrit
%
nentry=2;
while nentry > 1
    nentry=menu('Secondary Eyewall Stats','Stop','SER Probability','SER vs Wind Speed', ...
        'SER vs Lifetime-Integrated Wind Speed','PER Probability','PER vs Wind Speed', ...
        'PER vs Lifetime-Integrated Wind Speed');
    if nentry == 1
        break
    elseif nentry == 2
        y=us./n;
        bar(x,y,0.8)
        set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        set(gca,'xlim',[x(1)-1 semax+1])
        xlabel('Number of Secondary Replacement Cycles')
        ylabel('Probability')
        title('Probability of Secondary Eyewall Replacements')
        %
    elseif nentry == 3 || nentry == 4   
        %
        vse=zeros(1,qm);
        vsestd=zeros(1,qm);
        vsum=zeros(1,qm);
        vsumstd=zeros(1,qm);
        vmaxa=max(vnet,[],2);
        sumvmax=sum(vnet,2);
        dur=2.*sum(min(vnet,1),2)/24;
        for i=1:qm
            ss=ismember(senum,x(i));
            ss=ss*facs(i);
            vsum(1,i)=mean(nonzeros(sumvmax.*ss));
            vsumstd(1,i)=std(nonzeros(sumvmax.*ss));
            vse(1,i)=mean(nonzeros(vmaxa.*ss));
            vsestd(1,i)=std(nonzeros(vmaxa.*ss));
        end
        vsum(vsum==1)=0;
        vsumstd(vsumstd==1)=0;
        vsum=vsum*2/24;
        vsumstd=vsumstd.*2./24;
        %
        if nentry == 3
            bar(x,vse,0.8)
            hold on
            errorbar(x,vse,min(vsestd,vse),vsestd,'linestyle','none','linewidth',2,'color','g')
            hold off
            set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
            set(gca,'xlim',[x(1)-1 semax+1])
            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     
            %
            bar(x,vsum,0.8)
            hold on
            errorbar(x,vsum,min(vsumstd,vsum),vsumstd,'linestyle','none','linewidth',2,'color','g')
            hold off
            set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
            set(gca,'xlim',[x(1)-1 semax+1])
            xlabel('Number of Secondary Replacement Cycles')
            ylabel('Mean Integral Lifetime Wind Speed (knot-days)')
            title('# Secondary Cycles vs. Mean Integral Lifetime Wind Speed')
            y=vsum;
        end    
        elseif nentry == 5
            y=up./n;
            bar(x,y,0.8)
            set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
            set(gca,'xlim',[x(1)-1 pemax+1])
            xlabel('Number of Primary Replacement Cycles')
            ylabel('Probability')
            title('Probability of Primary Eyewall Replacements')
            %
        elseif nentry == 6 || nentry == 7   
            %
            qm=max(size(x));
            vse=zeros(1,qm);
            vsestd=zeros(1,qm);
            vsum=zeros(1,qm);
            vsumstd=zeros(1,qm);
            vmaxa=max(vnet,[],2);
            sumvmax=sum(vnet,2);
            dur=2.*sum(min(vnet,1),2)/24;
            for i=1:qm
                ss=facp(i).*ismember(penum,x(i));
                vsum(1,i)=mean(nonzeros(sumvmax.*ss));
                vsumstd(1,i)=std(nonzeros(sumvmax.*ss));
                vse(1,i)=mean(nonzeros(vmaxa.*ss));
                vsestd(1,i)=std(nonzeros(vmaxa.*ss));
            end
            vsum=vsum*2/24;
            vsumstd=vsumstd.*2./24;
        %
        if nentry == 6
            bar(x,vse,0.8)
            hold on
            errorbar(x,vse,min(vsestd,vse),vsestd,'linestyle','none','linewidth',2,'color','g')
            hold off
            set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
            set(gca,'xlim',[x(1)-1 pemax+1])
            xlabel('Number of Primary Replacement Cycles')
            ylabel('Mean Lifetime Maximum Wind Speed (knots)')
            title('# Primary Cycles vs. Mean Lifetime Maximum Wind Speed')
            y=vse;
        %
        elseif nentry == 7     
            %
            bar(x,vsum,0.8)
            hold on
            errorbar(x,vsum,min(vsumstd,vsum),vsumstd,'linestyle','none','linewidth',2,'color','g')
            hold off
            set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
            set(gca,'xlim',[x(1)-1 pemax+1])
            xlabel('Number of Primary Replacement Cycles')
            ylabel('Mean Integral Lifetime Wind Speed (knot-days)')
            title('# Primary Cycles vs. Mean Integral Lifetime Wind Speed')
            y=vsum;
        end    
    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 