%
%  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. 
%-------------------------------------------------------------------------
%
params    %  Load parameters
%
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
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
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);
pdihigh=zeros(1,n);
pdilow=zeros(1,n);
pdimax=zeros(1,n);
vmaxsyn=zeros(1,n);
nshortdur=zeros(1,n);
duration=zeros(1,n);
durhur=zeros(1,n);
pdiland=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);
%
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(passi.*vnet.^3,2));
vtemp3=(1852.*vtemp./3600).^3;
vuu=ismember(round(vnet),64:220);
vuus=sum(vuu,2);
vuus=vuus';
%vuus=max(vuus,0.0001);
for k=1:n    
    pass=ismember(single(yearstore),yearset(k)); % Identifies members of the kth year
    fac=1./sum(pass);
    ntc(k)=freqyear(k);
    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;
    pdi(k)=freqyear(k).*fac.*pdifac.*sum(ss);
    pdihigh(k)=quantile(freqyear(k)*nonzeros(pdifac*ss),0.85);
    pdilow(k)=quantile(freqyear(k)*nonzeros(pdifac*ss),0.15);
    pdimax(k)=freqyear(k).*fac.*sum(pass.*vtemp3);
    pdiland(k)=freqyear(k).*fac.*pdifac.*sum(pass.*pdil);
    %nland(k)=freqyear(k).*fac.*sum(pass.*sum(min(jland',1),1)); % slow
    nland(k)=freqyear(k).*fac.*sum(pass.*min(pressl,1));
    %nlandhur(k)=freqyear(k).*fac.*sum(pass.*sum(min(jlandhur',1),1)); % slow
    nlandhur(k)=freqyear(k).*fac.*sum(pass.*uhur);
    %nlandmajor(k)=freqyear(k).*fac.*sum(pass.*sum(min(jlandmajor',1),1)); % slow
    nlandmajor(k)=freqyear(k).*fac.*sum(pass.*umajor);
    vmaxsyn(k)=max(pass.*vtemp);
    %vmaxsyn(k)=mean(nonzeros((pass.*vtemp)));
    pressland(k)=freqyear(k)*fac*sum((nonzeros(pass.*(pressl+0.001))));
    nshortdur(k)=freqyear(k).*fac.*sum(min(max(pass.*(24-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;
%
%  Save variables for next time routine called (added February, 2014)
%
%{
filename=char(strcat('time_series_files/',lower(tit1),'_time_series'));
filename=char(filename);
save(filename, 'nhur', 'ncat1', 'ncat2', 'ncat3', 'ncat4', 'ncat5', 'yearset', ...
    'x2', 'ntc', 'pdi', 'pdimax', 'pdiland', 'nland', 'vmaxsyn', 'nshortdur', ...
    'duration', 'yearmin', 'yearmax', 'nn', 'm', 'freqyear', 'nmajor', 'bas', ...
    'shape','rmyear','eset','pressland','vcrit') 
if exist('vs','var')
    save(filename,'vsseries','-append')
end    
if exist('router','var')
    save(filename,'routyear','-append')
end    
if exist('clat','var')
    save(filename,'city_radius','clat','clong','-append')
end    
%}
%
%
%-------------------------------------------------------------------------
bestproc   %  Process best tracks
%
%  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')
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);
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 Julian day of first track points
    %
    [yearbest,~,~,~,~,~]=datevec(datebest);
    dateyear=datenum(yearbest,1,1,0,0,0);
    julianbest=datebest-dateyear;
    clear dateyear 
    %
    [~,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);
    %
    yearbest1=transpose(yearbest(:,1));
    n=max(size(bestyears));
    if vswitch == 2
        vtemp=vmaxb;
    else
        vtemp=vsmmax;
    end    
%    
    pdfacb=transpose(sum(passib.*vbest.^3,2));
    vtemp3=(1852.*vtemp./3600).^3;
    for k=1:n    
        pass=ismember(yearbest1,bestyears(k)); % Identifies members of the kth year
        %ntcbest(k)=sum(pass);
        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.*sum(min(jlandhurbest',1),1));
        nlandmajorbest(k)=sum(pass.*sum(min(jlandmajorbest',1),1));
        vmaxbest(k)=max(pass.*vtemp);
        a=mean((nonzeros(pass.*(pressbest+0.001))));
        a(isempty(a))=1013;
        presslandbest(k)=a;
        nshortdurbest(k)=sum(min(max(pass.*(24-jmax),0),1));
        if max(nonzeros(pass.*jmax)) > 0
            durationbest(k)=2.*mean(nonzeros(pass.*jmax));
        end
        %
    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
    nplot=menu('Plot Options','Stop','Frequency','# Hurricanes','# Cat1s','# Cat2s','# Cat3s','# Cat4s', ...
        '# Cat5s','# Majors','Power dissipation','Storm max power dissipation','Duration', ...
        'All landfall power dissipation','Frequency of Short Duration (< 2 days) events', ...
        'Frequency of long duration (> 2 days) events', ...
        'Mean radius of maximum winds at time of peak intensity', ...
        'Mean outer radius at time of peak intensity', ...
        'Apply 1-3-4-3-1 smoother to last series plotted', ...
        'Apply 10-year running average to last series plotted', ...
        'Apply n-year running average to last series plotted (Enter n)');
    if nplot == 1
        stop(tmr9)
        delete(tmr9)
        clear tmr9
        break       
    elseif nplot == 2
        y=ntcbest;
        y2=ntc;
        ylab='Annual Frequency';
    elseif nplot == 3
        y=nhurbest;
        y2=nhur;
        ylab='Frequency of hurricanes';    
    elseif nplot == 4
        y=ncat1best;
        y2=ncat1;
        ylab='Frequency of Cat 1 storms';   
    elseif nplot == 5
        y=ncat2best;
        y2=ncat2;
        ylab='Frequency of Cat 2 storms';   
    elseif nplot == 6
        y=ncat3best;
        y2=ncat3;
        ylab='Frequency of Cat 3 storms';   
    elseif nplot == 7
        y=ncat4best;
        y2=ncat4;
        ylab='Frequency of Cat 4 storms';               
    elseif nplot == 8
        y=ncat5best;
        y2=ncat5;
        ylab='Frequency of Cat 5 storms';    
    elseif nplot == 9
        y=nmajorbest;
        y2=nmajor;
        ylab='Frequency of Major Hurricanes';    
    elseif nplot == 10
        y=pdibest;
        y2=pdi;
        ylab='Power dissipation index';    
    elseif nplot == 11
        y=pdimaxbest;
        y2=pdimax;
        ylab='Storm max power dissipation';  
    elseif nplot == 12
        y=durationbest/24;
        y2=duration/24;
        ylab='Mean duration (days)';    
    elseif nplot == 13
        y=pdilandbest;
        y2=pdiland;
        ylab='All landfall power dissipation';     
    elseif nplot == 14
        y=nshortdurbest;
        y2=nshortdur;
        ylab='Frequency of short duration events';      
    elseif nplot == 15
        y=ntcbest-nshortdurbest;
        y2=ntc-nshortdur;
        ylab='Frequency of long duration events'; 
    elseif nplot == 16
        y=rmyear;
        y2=rmyear;
        ylab='Mean radius of maximum winds (km) at time of peak intensity';         
    elseif nplot == 17
        y=routyear;
        y2=routyear;
        ylab='Mean outer radius (km) at time of peak intensity'; 
    end    
    %
    if strcmp(besttrackplot,'y') && nplot ~= 16 && nplot ~= 17
        [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));
            tlabel=['r^2 = ' num2str(r2)];
        end
        %r2=std(y2(i2)-y(i1));
        if (nplot >= 2 && nplot <= 9) || nplot == 14 ...
                || nplot == 15 
            h1=bar(x,y,0.5);
            hold on
            g=plot(x2,y2,'linewidth',2,'color','r');
            hold off
        elseif nplot < 16    
             h=plot(x,y,'b',x2,y2,'r');
             set(h,'linewidth',2)
        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 < 16
            legend('Observed','MIT Synthetic','location','best')
        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])
    %
    if nplot == 18
        hold on
        xs=smooth134(x2);
        y3=smooth134(y2);
        j6=plot(xs,y3);
        set(j6,'linewidth',2,'color','g');
        hold off
        if mean(y) == mean(routyear) || mean(y) == mean(rmyear)
            legend('Original','1-3-4-3-1 smoothed','location','best')
        else
            legend('Observed','MIT Synthetic','1-3-4-3-1 smoothed','location','best')
        end    
    end    
    if nplot == 19
        hold on
        xs=runav(x2,10);
        y3=runav(y2,10);
        j=plot(xs,y3,'g');
        set(j,'linewidth',2);
        hold off
        if mean(y) == mean(routyear) || mean(y) == mean(rmyear)
            legend('Original','10-year smoothed','location','best')
        else
            legend('Observed','MIT Synthetic','10-year smoothed','location','best')
        end            
    end    
    if nplot == 20
        hold on
        prompt = {'Enter number of years to average over'};
        dlg_title = 'Years to Average Over';
        num_lines = 1;
        answer = inputdlg(prompt,dlg_title,num_lines);
        if isempty(answer) || isempty(answer{1})
            return
        end    
        nsmooth=str2double(answer{1});
        xs=runav(x2,nsmooth);
        y3=runav(y2,nsmooth);
        j=plot(xs,y3,'g');
        set(j,'linewidth',2);
        hold off
        ltit=char(strcat(answer(1),'-year smoothed'));
        if mean(y) == mean(routyear) || mean(y) == mean(rmyear)
            legend('Original',ltit,'location','best')
        elseif strcmp(besttrackplot,'y')
            legend('Observed','MIT Synthetic',ltit,'location','best')
        end            
    end  
end 
%