%
%  Plots time series of various quantities.
%  Copyright WindRiskTech, L.L.C., December, 2011
%  Modified March, 2016
%  Modified January, 2020 to include outer radius calculation
%  Modified May, 2020 to speed up considerably
%  Modified September, 2021 to delete option of using previously saved file
%  Modified 4 February 2024 to substantially speed up code. Note that the
%  variables nland, nhur and nmajor, which were and are calculated but not
%  plotted, may not have been accurate before this. 
%  Modified May, 2024 to plot sampling confidence intervals and to smooth
%  series in time. Also included some new variables. Dropped additional
%  smoothing options. 
%  Modified March, 2026 to alow for a lower threshold in duration
%--------------------------------------------------------------------------
%
params    %  Load parameters
jcrit=round(durcrit/2);      % Critical number of 2-hour pts
%
clf('reset')
days=[31 28 31 30 31 30 31 31 30 31 30 31];
daysum=zeros(1,12);
p=0.05:0.05:0.95;
ps=max(size(p));
for i=1:12
    daysum(i)=sum(days(1:i));
end    
pi=acos(-1);
pi180=pi/180;
pi180i=1/pi180;
dfac=60*1852/(4*3600);
pdifacl=(1852/3600).^3;
pdifac=2.*3600.*pdifacl;
%
clear yearstore yearset passus udtemp vs pass T600store datea
warning ('off','all');
load temp yearstore passi
warning ('on','all');
if exist('yearstore','var') == 0
    errordlg('No time series available for this file','File Error');
    return
end    
load temp freqyear
%
% Detrend freqyear
%
%[b,bint,stats]=regressm(1:length(freqyear),freqyear);
%f=b(1)+b(2)*(1:length(freqyear));
%freqyear=freqyear+mean(f)-f;
%
tmr9=timer;
set(tmr9,'startdelay',3);
set(tmr9,'timerfcn','delete(hms)')
hms=msgbox('Calculating time series; please wait for menu to appear');
start(tmr9)
load sorted
%
%  Determine filename to store or retrieve time series data to accelerate routine (added
%  February, 2014)
%
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 
%
if exist('shearstore','var') == 0 
    load temp
end  
%
basin=lower(bas);
[nn,m]=size(vstore);
nni=1/nn;
%
%  Find the set of years in this event set
%
yearmax=max(yearstore);
yearmin=min(yearstore);
yeartemp=yearmin:yearmax;
aa=ismember(yeartemp,yearstore);
yearset=nonzeros(aa.*single(yeartemp));
yearset=transpose(yearset);
n=max(size(yearset));
x2=yearset;
%
% Initialize time series quantitites
%
ntc=zeros(1,n);
ncat1=zeros(1,n);
ncat2=zeros(1,n);
ncat3=zeros(1,n);
ncat4=zeros(1,n);
ncat5=zeros(1,n);
pdi=zeros(1,n);
%potdam=zeros(1,n);
pdimax=zeros(1,n);
pdiland=zeros(1,n);
vmaxsyn=zeros(1,n);
nshortdur=zeros(1,n);
duration=zeros(1,n);
durhur=zeros(1,n);
pressland=zeros(1,n);
nland=zeros(1,n);
nlandhur=zeros(1,n);
nlandmajor=zeros(1,n);
routyear=zeros(1,n);
rmyear=zeros(1,n);
if exist('vs','var')
    %
    udtemp=zeros(1,nn);
    vdtemp=zeros(1,nn);
    [vsmax,jsmax]=max(vs,[],2);
    for i=1:nn
        udtemp(i)=-vsmax(i).*sin(pi180*dir1(i,jsmax(i)));
        vdtemp(i)=-vsmax(i).*cos(pi180*dir1(i,jsmax(i)));
    end
    dirseries=zeros(1,n);
    %
    vsseries=zeros(1,n);
end    
%
% Calculate outer radius at time of maximum wind speed
%
[~,jmax]=max(vnet,[],2);
routmax=zeros(1,nn);
rmmax=zeros(1,nn);
for i=1:nn
    rmmax(i)=rmstore(i,jmax(i));
end    
if exist('router','var')
    for i=1:nn
        routmax(i)=router(i,jmax(i));
    end
end
%
% Calculate Julian day of first track points
%
julian=zeros(1,nn);
for i=1:nn
   if monthstore(i,1) == 1
       julian(i)=daystore(i,1);
   else
       julian(i)=daysum(monthstore(i,1)-1)+daystore(i,1);
   end
end    
[~,jmax]=min(vstore,[],2);
jmax=jmax-1;
jmax=transpose(jmax);
%  
[pdil, jland, jlandhur, jlandmajor] = landfallpdi(latstore,longstore,vnet,jmax);
[pressl,~] = landfallpress(latstore, longstore, pstore, jmax);
uhur=max(jlandhur,[],2);
uhur=uhur';uhur=min(uhur,1);
umajor=max(jlandmajor,[],2);
umajor=umajor';umajor=min(umajor,1);
%{
pdam=zeros(1,nn);
for i=1:nn
    for k=1:20
        if jland(i,k) > 1
            pdam(i)=pdam(i)+damage(vnet(i,jland(i,k)-1));
        end
    end
end
%}
fac=n.*nni;
clear vtemp
if vswitch == 2
    vtemp(1,:)=vmax;
else
    vtemp(1,:)=vnetmax;
end    
%delv=max(vnet(:,13:end)-vnet(:,1:end-12),[],2);
%delv=delv';
delv=vnet(:,13:end)-vnet(:,1:end-12);
pfac=transpose(sum(vnet.^3,2));
vtemp3=(1852.*vtemp./3600).^3;
vuu=ismember(round(vnet),64:220);
vuus=sum(vuu,2);
vuus=vuus';
%
% For non-integer quantities, estimate sampling error by sub-sampling
%
if strcmpi(samperr,'y')

    [pdilow,pdihigh] = subsample(pfac*pdifac,freqyear,double(yearstore),yearset,rsmooth);
    [pdimaxlow,pdimaxhigh] = subsample(vtemp3,freqyear,double(yearstore),yearset,rsmooth);
    [pdilandlow,pdilandhigh] = subsample(pdifac*pdil,freqyear,double(yearstore),yearset,rsmooth);
end
%
for k=1:n    
    pass=ismember(single(yearstore),yearset(k)); % Identifies members of the kth year
    pass=pass.*ismember(jmax,jcrit:400);
    jn=floor(0.5*(rsmooth-1));
    kmin=max(k-jn,1);
    kmax=min(k+jn,n);
    passci=ismember(single(yearstore),yearset(kmin:kmax));
    fac=1./sum(pass);
    ntc(k)=freqyear(k).*fac.*sum(min(max(pass.*(vtemp+1-str2double(vcrit)),0),1));
    ncat1(k)=freqyear(k).*fac.*sum(min(max(pass.*(vtemp-63),0),1));
    ncat2(k)=freqyear(k).*fac.*sum(min(max(pass.*(vtemp-82),0),1));
    ncat3(k)=freqyear(k).*fac.*sum(min(max(pass.*(vtemp-95),0),1));
    ncat4(k)=freqyear(k).*fac.*sum(min(max(pass.*(vtemp-112),0),1));
    ncat5(k)=freqyear(k).*fac.*sum(min(max(pass.*(vtemp-136),0),1));
    ss=pass.*pfac;
    ssci=passci.*pfac;
    pdi(k)=freqyear(k)*mean(nonzeros(pdifac.*ss));
    pdimax(k)=freqyear(k)*mean(nonzeros(pass.*vtemp3));
    pdiland(k)=freqyear(k)*pdifac.*mean(nonzeros(pass.*pdil));
    %potdam(k)=freqyear(k).*mean(nonzeros(pass.*pdam));
    nland(k)=freqyear(k).*fac.*sum(pass.*min(pressl,1));
    nlandhur(k)=freqyear(k).*fac.*sum(pass.*uhur);
    nlandmajor(k)=freqyear(k).*fac.*sum(pass.*umajor);
    vmaxsyn(k)=max(pass.*vtemp);
    aa=pass.*pressl;
    if max(aa) > 0
        pressland(k)=min(nonzeros(aa));
    end    
    nshortdur(k)=freqyear(k).*fac.*sum(min(max(pass.*(36-jmax),0),1));
    %
    durhur(k)=2.*mean(nonzeros(pass.*vuus));
    %
    if max(nonzeros(pass.*jmax)) > 0
        duration(k)=2.*mean(nonzeros(pass.*jmax));
    end
    %
    if exist('vs','var')
        vsseries(k)=mean(nonzeros(pass.*vsmax'));
        udd=mean(nonzeros(pass.*udtemp));
        vdd=mean(nonzeros(pass.*vdtemp));
        dirseries(k)=180+atan2(udd,vdd)*pi180i;
    end  
    rmyear(k)=mean(nonzeros(pass.*rmmax));
    if exist('router','var')
        routyear(k)=mean(nonzeros(pass.*routmax));
    end    
    %
end
nmajor=ncat3;
clear udtemp vdtemp udd vdd
%
% Transform from cumulative numbers to binned numbers
%
nts=ntc-ncat1;
nhur=ncat1;
ncat1=ncat1-ncat2;
ncat2=ncat2-ncat3;
ncat3=ncat3-ncat4;
ncat4=ncat4-ncat5;
%
%-------------------------------------------------------------------------
bestproc   %  Process best tracks
vcrit={num2str(str2double(vcrit)-5)};  % Subtract 5 knots from syntehtic critical value
%
%  Now check whether event years and best track years overlap
%
bestover=sum(ismember(yearset,bestyears));
if bestover >= 3
    besttrackplot='y';
else
    besttrackplot='n';
    msgbox('Not enough overlappiing best track years, proceeding without best tracks')
    samperr='n';
end    
if strcmp(besttrplot,'n')  %  If desired, omit best tracks in any case
    besttrackplot='n';
end    
%
yearset=bestyears;
nb=max(size(yearset));
x=yearset;
%
% Initialize best track time series quantitites
%
ntcbest=zeros(1,nb);
nhurbest=zeros(1,nb);
ncat1best=zeros(1,nb);
ncat2best=zeros(1,nb);
ncat3best=zeros(1,nb);
ncat4best=zeros(1,nb);
ncat5best=zeros(1,nb);
pdibest=zeros(1,nb);
pdimaxbest=zeros(1,nb);
vmaxbest=zeros(1,nb);
nshortdurbest=zeros(1,nb);
durationbest=zeros(1,nb);
durhurbest=zeros(1,nb);
pdilandbest=zeros(1,nb);
presslandbest=zeros(1,nb);
nlandbest=zeros(1,nb);
nlandhurbest=zeros(1,nb);
nlandmajorbest=zeros(1,nb);
nmajorbest=zeros(1,nb);
%
% If there are at least 3 overlapping years, load best tracks and find
% time series
%
if strcmpi(besttrackplot,'y')
    %
    % Calculate Year of first track points
    %
    [yearbest,~,~,~,~,~]=datevec(datebest);
    %
    [~,jmax]=min(vbest,[],2);
    jmax=jmax-1;
    jmax=max(jmax,1);
    jmax=transpose(jmax);
    %  
    [pdilbest, jlandbest, jlandhurbest, jlandmajorbest] = landfallpdi(latbest,longbest,vbest,jmax);
    [pressbest,~] = landfallpress(latbest, longbest, pbest, jmax);
    uhurbest=max(jlandhurbest,[],2);
    uhurbest=uhurbest';uhurbest=min(uhurbest,1);
    umajorbest=max(jlandmajorbest,[],2);
    umajorbest=umajorbest';umajorbest=min(umajorbest,1);
    %
    yearbest1=transpose(yearbest(:,1));
    n=max(size(bestyears));
    if vswitch == 2
        vtemp=vmaxb;
    else
        vtemp=vsmmax;
    end    
%    
    pdfacb=transpose(sum(vbest.^3,2));
    vtemp3=(1852.*vtemp./3600).^3;
    vuu=ismember(round(vbest),64:220);
    vuus=sum(vuu,2);
    vuus=vuus';
    for k=1:n    
        pass=ismember(yearbest1,bestyears(k)); % Identifies members of the kth year
        pass=pass.*ismember(jmax,jcrit:400);
        ntcbest(k)=sum(min(max(pass.*(vtemp-str2double(vcrit)),0),1));
        ncat1best(k)=sum(min(max(pass.*(vtemp-63),0),1));
        ncat2best(k)=sum(min(max(pass.*(vtemp-82),0),1));
        ncat3best(k)=sum(min(max(pass.*(vtemp-95),0),1));
        ncat4best(k)=sum(min(max(pass.*(vtemp-112),0),1));
        ncat5best(k)=sum(min(max(pass.*(vtemp-136),0),1));
        pdibest(k)=pdifac.*sum(pass.*pdfacb);
        pdimaxbest(k)=sum(pass.*vtemp3);
        pdilandbest(k)=pdifac.*sum(pass.*pdilbest);
        nlandbest(k)=sum(pass.*sum(min(jlandbest',1),1));
        nlandhurbest(k)=sum(pass.*uhurbest);
        nlandmajorbest(k)=sum(pass.*umajorbest);
        vmaxbest(k)=max(pass.*vtemp);
        a=min((nonzeros(pass.*pressbest)));
        a(isempty(a))=1013;
        presslandbest(k)=a;
        nshortdurbest(k)=sum(min(max(pass.*(36-jmax),0),1));
        if max(nonzeros(pass.*jmax)) > 0
            durationbest(k)=2.*mean(nonzeros(pass.*jmax));
        end
        durhurbest(k)=2.*mean(nonzeros(pass.*vuus));
        %
    end
    presslandbest(isnan(presslandbest))=1013;
    nmajorbest=ncat3best;
    %clear vtemp
    %
    % Transform from cumulative numbers to binned numbers
    %
    nhurbest=ncat1best;
    ncat1best=ncat1best-ncat2best;
    ncat2best=ncat2best-ncat3best;
    ncat3best=ncat3best-ncat4best;
    ncat4best=ncat4best-ncat5best;
%
end
%
nplot=2;
while nplot > 1
    clear Errbar
    int=0;
    r2p=1;
    nplot=menu('Plot Options','Stop','Frequency','# Hurricanes','# Cat1s','# Cat2s','# Cat3s','# Cat4s', ...
        '# Cat5s','# Majors','# Landfalling TCs','# Landfalling hurricanes','# Landfalling major hurricanes', ...
        'Power dissipation','Storm max power dissipation', ...
        'All landfall power dissipation','Annual maximum wind speed', 'Annual minimum landfall pressure', ...
        'Duration','Frequency of Short Duration (< 2 days) events', ...
        'Frequency of long duration (> 2 days) events', ...
        'Frequency of duration at hurricane strength', ...
        'Mean radius of maximum winds at time of peak intensity', ...
        'Mean outer radius at time of peak intensity');
    if nplot == 1
        stop(tmr9)
        delete(tmr9)
        clear tmr9
        break       
    elseif nplot == 2
        y=ntcbest;
        y2=ntc;
        ylab='Annual Frequency';
        int=1;
    elseif nplot == 3
        y=nhurbest;
        y2=nhur;
        ylab='Frequency of hurricanes';    
        int=1;
    elseif nplot == 4
        y=ncat1best;
        y2=ncat1;
        ylab='Frequency of Cat 1 storms';   
        int=1;
    elseif nplot == 5
        y=ncat2best;
        y2=ncat2;
        ylab='Frequency of Cat 2 storms';   
        int=1;
    elseif nplot == 6
        y=ncat3best;
        y2=ncat3;
        ylab='Frequency of Cat 3 storms';   
        int=1;
    elseif nplot == 7
        y=ncat4best;
        y2=ncat4;
        ylab='Frequency of Cat 4 storms';     
        int=1;
    elseif nplot == 8
        y=ncat5best;
        y2=ncat5;
        ylab='Frequency of Cat 5 storms';    
        int=1;
    elseif nplot == 9
        y=nmajorbest;
        y2=nmajor;
        ylab='Frequency of Major Hurricanes';    
        int=1;
    elseif nplot == 10
        y=nlandbest;
        y2=nland;
        ylab='Frequency all landfalling TCs';    
        int=1;
     elseif nplot == 11
        y=nlandhurbest;
        y2=nlandhur;
        ylab='Frequency all landfalling hurricanes';    
        int=1;
     elseif nplot == 12
        y=nlandmajorbest;
        y2=nlandmajor;
        ylab='Frequency all landfalling major hurricanes';    
        int=1;
    elseif nplot == 13
        y=pdibest;
        y2=pdi;
        ylab='Power dissipation index';    
        int=2;
    elseif nplot == 14
        y=pdimaxbest;
        y2=pdimax;
        ylab='Storm max power dissipation';  
        int=3;
    elseif nplot == 15
        y=pdilandbest;
        y2=pdiland;
        ylab='All landfall power dissipation'; 
        int=4;
    elseif nplot == 16
        y=vmaxbest;
        y2=vmaxsyn;
        ylab='Annual maximum wind speed (knots)';
    elseif nplot == 17
        y=presslandbest;
        y2=pressland;
        ylab='Mimimum landfall pressure (hPa)';     
    elseif nplot == 18
        y=durationbest/24;
        y2=duration/24;
        ylab='Mean duration (days)';    
        int=1;
    elseif nplot == 19
        y=nshortdurbest;
        y2=nshortdur;
        ylab='Frequency of short duration (< 36 hr) events';      
        int=1;
    elseif nplot == 20
        y=ntcbest-nshortdurbest;
        y2=ntc-nshortdur;
        ylab='Frequency of long duration (> 36 hr) events'; 
        int=1;
    elseif nplot == 21
        y=durhurbest/24;
        y2=durhur/24;
        ylab='Mean duration at hurricane intensity (days)';    
        int=1;
    elseif nplot == 22
        y=rmyear;
        y2=rmyear;
        ylab='Mean radius of maximum winds (km) at time of peak intensity';         
    elseif nplot == 23
        y=routyear;
        y2=routyear;
        ylab='Mean outer radius (km) at time of peak intensity'; 
    end    
    %
    ym=mean(y); y2m=mean(y2);
    y=movmean(y,rsmooth); y2=movmean(y2,rsmooth);  % Apply smoothing
    y=y*ym/mean(y);   y2=y2*y2m/mean(y2);          % Unfortunately, movmean changes sample mean. This restores it. 
    %y=y*mean(y2)/mean(y);                         % If you want time series normalized so that observed and modeled 
                                                   % means are equal
    if strcmp(besttrackplot,'y') && nplot < 22 
        [x3,i1,i2]=intersect(x,x2);
        if max(size(i1)) > 3 && max(size(nonzeros(y(i1)))) > 1
            [c,f,r2]=regressnpl(y(i1),y2(i2));
            r2=0.01*round(100*r2);
            tlabel=['r^2 = ' num2str(r2)];
        end
        if strcmpi(samperr,'n')
            int=0;
        end    
        if int == 1
            Errbar(1,:)=poissinv(0.95,rsmooth*y2)/rsmooth-y2;
            Errbar(2,:)=y2-poissinv(0.05,rsmooth*y2)/rsmooth;
            k2=shadedErrorBar(x2,y2,Errbar,{'r','linewidth',2},1);
            hold on
            g2=plot(x,y,'b','linewidth',2);
            lvec=[g2 k2.mainLine k2.patch];
            hold off
        elseif int == 2
            Errbar(1,:)=pdihigh-y2;
            Errbar(2,:)=y2-pdilow;
            k2=shadedErrorBar(x2,y2,Errbar,{'r','linewidth',2},1);
            hold on
            g2=plot(x,y,'b','linewidth',2);
            lvec=[g2 k2.mainLine k2.patch];
            hold off
        elseif int == 3
            Errbar(1,:)=pdimaxhigh-y2;
            Errbar(2,:)=y2-pdimaxlow;
            k2=shadedErrorBar(x2,y2,Errbar,{'r','linewidth',2},1);
            hold on
            g2=plot(x,y,'b','linewidth',2);
            lvec=[g2 k2.mainLine k2.patch];
            hold off
        elseif int == 4
            Errbar(1,:)=pdilandhigh-y2;
            Errbar(2,:)=y2-pdilandlow;
            k2=shadedErrorBar(x2,y2,Errbar,{'r','linewidth',2},1);
            hold on
            g2=plot(x,y,'b','linewidth',2);
            lvec=[g2 k2.mainLine k2.patch];
            hold off
        else    
             h=plot(x,y,'b',x2,y2,'r');
             set(h,'linewidth',2)
             lvec=h;
        end    
        maxvar=max(max(y),max(y2));
        minvar=min(min(y),min(y2));
        ypos=double(minvar+0.8*(maxvar-minvar));
        xpos=min(x)+0.2*(max(max(x),max(x2))-min(min(x),min(x2)));
        xpos=double(xpos);
        if max(size(i1)) > 3 && max(size(nonzeros(y(i1)))) > 1 ...
                && mean(y) ~= mean(routyear) && mean(y) ~= mean(rmyear)
            text(xpos,ypos,tlabel,'fontweight','bold')
        end
        if nplot < 22
            tit1=[num2str(100*plow),'^{th} to ',num2str(100*phigh),'^{th} percentile'];
            if (nplot < 16 || nplot > 17) && strcmpi(samperr,'y')
                legend(lvec,'Observed','MIT Synthetic',tit1,'location','best')
            else
                legend(lvec,'Observed','MIT Synthetic','location','best')
            end    
        end
    else
        h=plot(x2,y2,'r');
        set(h,'linewidth',2)
    end
    set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
    xlabel('Year','fontweight','bold')
    ylabel(ylab,'fontweight','bold')
    title(eset,'interpreter','none')
    if strcmp(besttrackplot,'y')     
        yearbeg=min(min(x2),min(x));
        yearend=max(max(x2),max(x));
    else
        yearbeg=min(x2);
        yearend=max(x2);
    end    
    set(gca,'xlim',[yearbeg-2 yearend+2])
    % 
end 
%--------------------------------------------------------------------------
function [y05,y95] = subsample(y,freqyear,yearstore,yearset,rsmooth)

%  Written by K. Emanuel and ChatGPT, February, 2026

ny = numel(yearset);

% Collect indices for each year
idx = cell(ny,1);
L   = zeros(ny,1);
for i = 1:ny
    idx{i} = find(yearstore == yearset(i));
    L(i)   = numel(idx{i});
end
Lmax = max(L);

% Build matrix: rows = within-year index, cols = year (padded with NaN)
Y = NaN(Lmax, ny);
for i = 1:ny
    Y(1:L(i), i) = freqyear(i) * y(idx{i});
end

% Smooth across years (dim 2), then quantiles over within-year samples (dim 1)
Ysm = movmean(Y, rsmooth, 2, 'omitnan');
y05 = quantile(Ysm, 0.05, 1);
y95 = quantile(Ysm, 0.95, 1);

end