clear all
close all
clc

load nino_ind3 % load ENSO index 'tim' vs 'ind'
po=ind; time=tim; %redefine variables (not really necessary, but to make it
% consistent with the rest of the script

[p1,coef]=lanczos(po,1*365.25*24,30.4375*1440); % 1-yr filtered
p1=detrend(p1); % detrend 1-yr filtered
[p3,coef]=lanczos(po,3*365.25*24,30.4375*1440); % 3-yr filtered
p3=detrend(p3); % detrend 3-yr filtered
[p5,coef]=lanczos(po,5*365.25*24,30.4375*1440);
p5=detrend(p5); % detrend 5-yr filtered

figure % generate Figure 1, upper panel, with different filtered series
set(gcf, 'Position',  [200, 50, 1100, 850])
ax1=axes('Position',[0.15 0.775 0.8 0.22]);
axes(ax1);
plot(time,po,'k','LineWidth',1) % ENSO
hold on
plot(time,p1,'m','LineWidth',3) % 1-yr filtered
plot(time,p3,'r','LineWidth',4) % 3-yr filtered
plot(time,p5,'b','LineWidth',5) % 5-yr filtered
xt=[1880:10:2020]; xticks(xt); % ticks to draw on the 'x' axis
grid on;  grid minor;  set(gca,'FontSize',18,'XTickLabel',[])
axis([1870 2024 -3 3])
text(1910,2.6,'Monthly','color','k','FontSize',18)
text(1945,2.6,'1-yr filtered','color','m','FontSize',18)
text(1980,2.6,'3-yr filtered','color','r','FontSize',18)
text(1920,-2.5,'5-yr filtered','color','b','FontSize',18)
text(1875,-2.25,'a)','FontSize',18)

% still on Figure 1
% Periods for the Fits
TN=18.61; % Nodal Precession
TA=8.85; % Apsidal Precession
AP=6; % Argument of Perigee
S1=10.5; %Solar Activity 1 
S2=11.5; % Solar Activity 2
TN2=TN/2; % Half Nodal Precession
TNN=TN*2; TN3=TN*3; TN6=TN*6;
TA2=TA/2; % Half Apsidal Precession
AP2=AP/2; %Half Argument of Perigee
TP2=(TN+TA)/2; %Average Precession
TNA=1/(1/TA-1/TN); % ~16.87 % Apsidal - Nodal Modulation
TNS1=1/(1/S1-1/TN); %Solar 1 - Nodal Modulation
TNS2=1/(1/S2-1/TN); %Solar 2 - Nodal Modulation
TAS1=1/(1/TA-1/S1); %Solar 1 - Apsidal Modulation
TAS2=1/(1/TA-1/S2); %Solar 2 - Apsidal Modulation
TPS1=1/(1/AP-1/S1); %Solar 1 - Argument of Perigee
TPS2=1/(1/AP-1/S2); %Solar 2 - Argument of Perigee

per=1./[TN,TPS2,TNS2,TNA,S1*2,TN6,S1,S2,TN2,TN3,TP2,TA,AP]; m=length(per); % multiples better
omega=2*pi*per';

%add best fits, derived from 'halessqfit', on top of filtered series for
% the other panels of Figure 1
cc=[.1 .1 .5]; % color of fitted line
ax1=axes('Position',[0.15 0.1 0.8 0.22]);
axes(ax1);
plot(time,p5,'b','LineWidth',3) % 5-yr filtered 
[meanS,amp,phase,r5,s5]=halessqfit(omega,time,p5,1,'b',cc); % best fit
axis([1870 2024 -.7 .7])
  xlabel('year')
 text(1910,-.5,'ENSO 5-yr filtered','color','b','FontSize',17)
 text(1945,.5,strcat('Reconstruction R^2 = ' ,num2str(r5,2)),'color',cc,'FontSize',17)
 set(gca,'FontSize',18)
 text(1875,-0.55,'d)','FontSize',18)
 
%panel for 3-yr filtered
cc=[.6 .1 .1]; % color of fitted line
ax1=axes('Position',[0.15 0.325 0.8 0.22]);
axes(ax1);
[mean3,am3,phas3,r3,s3]=halessqfit(omega,time,p3,1,'r',cc);
axis([1870 2024 -1 1])
 text(1910,-.8,'ENSO 3-yr filtered','color','r','FontSize',17)
 text(1945,.75,strcat('Reconstruction R^2 = ' ,num2str(r3,2)),'color',cc,'FontSize',17)
 text(1875,-0.8,'c)','FontSize',18)
set(gca,'FontSize',18,'XTickLabel',[])

%panel for 1-yr filtered
cc=[.1 .7 .1]; % color for filtered line
ax1=axes('Position',[0.15 0.55 0.8 0.22]);
axes(ax1);
[mean1,am1,phas1,r1,s1]=halessqfit(omega,time,p1,1,'m',cc);
axis([1870 2024 -2.5 2.5])
 text(1910,-1.75,'ENSO 1-yr filtered','color','m','FontSize',17)
 text(1920,1.8,strcat('Reconstruction R^2 = ' ,num2str(r1,2)),'color',cc,'FontSize',17)
set(gca,'FontSize',18,'XTickLabel',[])
text(-.07,-0.4,'ENSO Index','Rotation',90,'Units','normalized','FontSize',18)
text(1875,-2.1,'b)','FontSize',18)

% Figure 4 projection until year 2050
nx=length(time);
timep=[1870:.1:2100]';
gomexp=0*timep+meanS; % to reconstruct time series 5-yr until 2100
gomex3=0*timep+mean3; % to reconstruct time series 3-yr until 2100
  m=length(omega);
  for i=1:m
    gomexp=gomexp+amp(i)*sin(omega(i)*timep+phase(i)*pi/180);
    gomex3=gomex3+am3(i)*sin(omega(i)*timep+phas3(i)*pi/180);
  end
  
figure % draw the figure
set(gcf, 'Position',  [200, 100, 1100, 650])
  plot(timep,gomexp,'w','LineWidth',5);
  ylabel('ENSO Index')
  xlabel('year')
  hold on
  plot(time,p5,'w','LineWidth',3)
  plot(timep,gomex3,'color',[.6 .1 .1],'LineWidth',3); % reconstructed 3-yr
  fill([time(1),time',time(nx)],[-0.8,p5',-0.8],'b','LineStyle','None')
  alpha(0.55)
  plot(timep,gomexp,'color',[.1 .1 .6],'LineWidth',3); % reconstructed 5-yr
    axis([1870 2050 -.8 1])
  text(1940,.65,'3-yr filtered projection','color',[.6 .1 .1],'FontSize',18)
  text(1940,.8,'5-yr filtered projection','color',[.1 .1 .6],'FontSize',18)
  text(1900,-.625,'Period of observations -- 5-yr filtered','color','w','FontSize',17)
  quiver(1870,-0.7,166,0,'color','w','LineWidth',2,'MaxHeadSize',.003);
  quiver(2024,-0.7,-166,0,'color','w','LineWidth',2,'MaxHeadSize',.003);
  grid on
  grid minor
  set(gca,'FontSize',18)

  % Figure 2a AMPLITUDES OF HARMONICS
figure
cc=[.6 .1 .1]; % color for 3-yr bars
cb=[.1 .1 .6]; % color for 5-yr bars
set(gcf, 'Position',  [200, 100, 1300, 650])
la=['  N';'AS2';'NS2';' NP';'2S1';' 6N';' S1';' S2';' N2';' 3N';'TP2';'  P';'  A']; % labels
P=1./per; % periods of each harmonic to draw (defined on line 54)
  semilogx(P,amp,'*','color',cb,'LineWidth',10); % plot stars on top of each
  % 5-yr amplitude bar
  axis([5 120 0 .2])
  xticks([5 6 7 8 9 10 20 30 40 50 60 70 80 90 100 120])
  xticklabels({'5','6','7','8','9','10','20','30','40','50','60','70','80','90','100','120'})
  ylabel('Amplitude of Contribution to ENSO Index')
  xlabel('Period of Contributing Harmonic (y)')
  hold on
  for k=1:m % draw each bar of 5-yr amplitude
      semilogx([P(k),P(k)],[0,amp(k)],'color',cb,'LineWidth',12)
  end
  semilogx(P,am3,'o','color',cc,'LineWidth',5) % draw 'o' on top of each
  % 3-yr amplitude bars
  for k=1:m; % draw each bar of 3-yr amplitude
      semilogx([P(k),P(k)],[0,am3(k)],'color',cc,'LineWidth',6)
  end
  grid on
  grid minor
  text(40,.17,'ENSO 5-yr filtered','color',cb,'FontSize',18)
  text(40,.15,'ENSO 3-yr filtered','color',cc,'FontSize',18)
  text(5.5,.18,'a)','FontSize',18)
  text(P,am3*1.1,la,'verticalAlignment','bottom','horizontalAlignment','center','FontSize',18)
  set(gca,'FontSize',18)
 % end of Figure 2a 
  
 % reconstruction of 5-yr ENSO signal with contours -- for Figure 2b
cc=[.1 .1 .5];
 
for k=1:m % fit with each harmonic; from 1 to m harmonics
  ome=omega(1:k);
    [meanx,ampx,phasx,rx,skx]=halessqfit2(ome,time,p5,1,'b',cc); % best fit
     for i=1:length(ome); % reconstruct signal with i harmonics
      meanx=meanx+ampx(i)*sin(ome(i)*time+phasx(i)*pi/180); % reconstruct signal with i harmonics
     end
   sl(:,k)=meanx; % build matrix of reconstructed signal
 r(k)=rx; % build vectors of R^2
end

sl(:,m+1)=p5; % last row of matrix is the 5-yr ENSO
y=[1:1:m+1]; r(m+1)=1; % complete 'y' axis

% plot contours;  figure 2b
 figure % plot contours of increasingly better reconstructions
set(gcf, 'Position',  [100, 100, 1450, 550])
ax1=axes('Position',[0.05 .15 0.79 0.75]);
axes(ax1);
colormap ('redblue')

pcolor(tim,y,sl')
shading interp
alpha(0.7)
 grid on;  grid minor; 
set(gca,'FontSize',18)
hold on
contour(time,y,sl',[.0001,0],'w')
yticks(y)
yticklabels({'N','AS2','NS2','NP','2S1','6N','S1','S2','N2','3N','TP2','P','A','O'})
title('5-yr-ENSO Index)') 
colorbar('westoutside') ; alpha(0.7)
text(1875,2,'b)','FontSize',18)
xlabel('Year')

% plot R^2 as a function of increased number of harmonics; still Fig 2b
ax1=axes('Position',[0.85 0.15 0.1 0.75]);
axes(ax1);
plot(r,y,'+','LineWidth',5)
text(.7,2,'R^2','FontSize',18,'FontWeight','bold')
axis([0 1 1 m+1]);  grid on;  grid minor; 
set(gca,'FontSize',18,'YTickLabel',[])
ylabel('Name of Harmonic on Fit')
set(gca,'YAxisLocation','right')
