% uses lanczos.m & frlsqfit.m
clear all
close all

load sydney % data in Matlab format as provided by Chari

%organize variables 
mo=(sydney(:,1)-1)/12; % month vector
yr=sydney(:,2); % year vector
tim=yr+mo; % time vector (year and fraction thereof)

fmin=sydney(:,5); % monthly minima water level
fmax=sydney(:,6); % monthly maxima water level
fmea=sydney(:,7); % monthly mean sea level

% convert nans to zeros to produce detrended and filtered series
fmin(isnan(fmin)==1)=nanmean(fmin);
fmax(isnan(fmax)==1)=nanmean(fmax);
fmea(isnan(fmea)==1)=nanmean(fmea);

% detrend the 3 series
fmn=detrend(fmin);
fmx=detrend(fmax);
fmm=detrend(fmea);

% filter detrended series at sf years
sf=5; dt=365.25/12; % filter period and sampling interval in days
[fm5,coef]=lanczos(fmm,sf*365.25*24,dt*1440); % mean 5-yr
[fm3,coef]=lanczos(fmm,3*365.25*24,dt*1440); % mean 3-yr
[fn5,coef]=lanczos(fmn,sf*365.25*24,dt*1440); % min 5-yr
[fn3,coef]=lanczos(fmn,3*365.25*24,dt*1440); % min 3-yr
[fx5,coef]=lanczos(fmx,sf*365.25*24,dt*1440); % max 5-yr
[fx3,coef]=lanczos(fmx,3*365.25*24,dt*1440); % max 3-yr

% fit known periodicities to the 3 filtered time series; plot original and
% astronomic periods
TN=18.61; % Nodal Precession
TA=8.85; % Apsidal Precession
AP=6; % Argument of Perigee
S1=10.5; %Solar Activity 1 10.5
S2=11.5; % Solar Activity 2 11.5
TN2=TN/2; % Half Nodal Precession
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

% maximum
 per=[TN,AP,S1,S2,TA,3*TN,2*TN,TNS1,TNS2,TPS1,TPS2,TNA,6*TN]; omega=2*pi./per'; % actual periods to fit
 [meanx,ampx,phasx,rx,skx,ffx]=frlsqfit(omega,tim,fx5,1,'r','TN,TA,AP,S1,S2,TA2,AP2,TP2,TNS1,TNS2,TAS1,TAS2,TPS2,S1S2','Year');
lae=strcat('Sydney, Mean Monthly Maxima ', ' R^2 = ' ,num2str(rx,2) ,' ;  Skill = ' ,num2str(skx,2)); text(1920,.07,lae,'FontSize',18)%

% min
[meann,ampn,phasn,rn,skn,ffn]=frlsqfit(omega,tim,fn5,1,'b','TN,TA,AP,S1,S2,TA2,AP2,TP2,TNS1,TNS2,TAS1,TAS2,TPS2,S1S2','Year');
lae=strcat('Sydney, Mean Monthly Minima ',' R^2 = ' ,num2str(rn,2) ,' ;  Skill = ' ,num2str(skn,2)); text(1920,.07,lae,'FontSize',18)%

% mean
[meanm,ampm,phasm,rm,skm,ffm]=frlsqfit(omega,tim,fm5,1,'m','TN,TA,AP,S1,S2,TA2,AP2,TP2,TNS1,TNS2,TAS1,TAS2,TPS2,S1S2','Year');
lae=strcat('Sydney, Monthly Mean ',' R^2 = ' ,num2str(rm,2) ,' ;  Skill = ' ,num2str(skm,2)); text(1920,.07,lae,'FontSize',18)%
amp5x=ampx; amp5n=ampn; amp5m=ampm;% for max, min & mean

% figure 4c, d
% ==========================================================
% reconstruct 5-yr signals with parameters from fit
    m=length(per);
    for i=1:m
    meanx=meanx+ampx(i)*sin(omega(i)*tim+phasx(i)*pi/180);
    meann=meann+ampn(i)*sin(omega(i)*tim+phasn(i)*pi/180);
    meanm=meanm+ampm(i)*sin(omega(i)*tim+phasm(i)*pi/180);
    end
      
figure
set(gcf, 'Position',  [20, 250, 1250, 420])
ax1=axes('Position',[0.56 0.18 0.37 0.8]);
axes(ax1);
% plot means 5-yr
plot(tim,meanm,'m','LineWidth',3.) % means fit 5-yr
hold on
 plot(tim,fm5,'m','LineWidth',1) % original mean 5-yr
 s=0; pp=find(fm5 > 0); nn=find(fm5 < 0); % shade positives and negatives
fill([tim(1),tim(pp)',2024],[s,fm5(pp)'+s,s],'m','LineStyle','None')
fill([tim(1),tim(nn)',2024],[s,fm5(nn)'+s,s],'c','LineStyle','None')
alpha(.1)
% plot minima 5-yr
s=-.05;
plot(tim,meann+s,'b','LineWidth',3.) % fit 5-yr
plot(tim,fn5+s,'b','LineWidth',1) % min 5-yr
pp=find(fn5 > 0); nn=find(fn5 < 0); 
fill([tim(1),tim(pp)',2024],[s,fn5(pp)'+s,s],'m','LineStyle','None')
fill([tim(1),tim(nn)',2024],[s,fn5(nn)'+s,s],'c','LineStyle','None')
alpha(.1)
% plot maxima 5-yr
s=.05;
plot(tim,meanx+s,'r','LineWidth',3.); % fit 5-yr
plot(tim,fx5+s,'r','LineWidth',1) % max 5-yr
pp=find(fx5 > 0); nn=find(fx5 < 0); 
fill([tim(1),tim(pp)',2024],[s,fx5(pp)'+s,s],'m','LineStyle','None')
fill([tim(1),tim(nn)',2024],[s,fx5(nn)'+s,s],'c','LineStyle','None')
alpha(.1)
axis([1900 2025 -.1 .1]);  grid on;  grid minor;  
xlabel('Year'); ylabel('5-yr filtered (m)')
text(1917,0.085,strcat('Max; ',' R^2 = ' ,num2str(rx,2),'; S = ' ,num2str(skx,2)),'color','red','FontSize',18 )
text(1910,-0.09,strcat('Min; ',' R^2 = ' ,num2str(rn,2),'; S = ' ,num2str(skn,2)),'color','blue','FontSize',18 )
text(2024.75,0.04,'Mean;','color','magenta','FontSize',18 )
text(2023.25,0.025,strcat(' R^2 = ' ,num2str(rm,2)),'color','magenta','FontSize',18 )
text(2023.25,0.01,strcat(' S = ' ,num2str(skm,2)),'color','magenta','FontSize',18 )
set(gca,'FontSize',18)
text(2015,-.085,'d)','FontSize',18,'FontWeight','bold')
%===================================================================
% Do the same for 3-yr fits
%================================================================
[meanx,ampx,phasx,rx,skx,ffx]=frlsqfit(omega,tim,fx3,1,'r','TN,AP,S1,S2,TN2,TAS1,TAS2,TNS1,TNS2,TPS1,TPS2,TNA,S1S2','Year');
[meann,ampn,phasn,rn,skn,ffn]=frlsqfit(omega,tim,fn3,1,'b','TN,AP,S1,S2,TN2,TAS1,TAS2,TNS1,TNS2,TPS1,TPS2,TNA,S1S2','Year');
[meanm,ampm,phasm,rm,skm,ffm]=frlsqfit(omega,tim,fm3,1,'m','TN,AP,S1,S2,TN2,TAS1,TAS2,TNS1,TNS2,TPS1,TPS2,TNA,S1S2','Year');
amp3x=ampx; amp3n=ampn; amp3m=ampm;% for max, min & mean

% reconstruct 3-yr signals with parameters from fit
    m=length(per);
    for i=1:m
    meanx=meanx+ampx(i)*sin(omega(i)*tim+phasx(i)*pi/180);
    meann=meann+ampn(i)*sin(omega(i)*tim+phasn(i)*pi/180);
    meanm=meanm+ampm(i)*sin(omega(i)*tim+phasm(i)*pi/180);
    end
    
ax1=axes('Position',[0.06 0.18 0.37 0.8]);
axes(ax1);
s=0; % plot means 3-yr
plot(tim,meanm,'m','LineWidth',3) % fit 3-yr
hold on
 plot(tim,fm3,'m','LineWidth',1) % mean 3-yr
 pp=find(fm3 > 0); nn=find(fm3 < 0); 
fill([tim(1),tim(pp)',2024],[s,fm3(pp)'+s,s],'m','LineStyle','None')
fill([tim(1),tim(nn)',2024],[s,fm3(nn)'+s,s],'c','LineStyle','None')
alpha(.1)
%
s=-.08; % plot minima 3-yr
plot(tim,meann+s,'b','LineWidth',3.) % fit 3-yr
plot(tim,fn3+s,'b','LineWidth',1) % min 3-yr
pp=find(fn3 > 0); nn=find(fn3 < 0); 
fill([tim(1),tim(pp)',2024],[s,fn3(pp)'+s,s],'m','LineStyle','None')
fill([tim(1),tim(nn)',2024],[s,fn3(nn)'+s,s],'c','LineStyle','None')
alpha(.1)
%
s=.08; % plot maxima 3-yr
plot(tim,meanx+s,'r','LineWidth',3.); % fit 3-yr
plot(tim,fx3+s,'r','LineWidth',1) % max 3-yr
pp=find(fx3 > 0); nn=find(fx3 < 0); 
fill([tim(1),tim(pp)',2024],[s,fx3(pp)'+s,s],'m','LineStyle','None')
fill([tim(1),tim(nn)',2024],[s,fx3(nn)'+s,s],'c','LineStyle','None')
alpha(.1)
axis([1900 2025 -.15 .15]);  grid on;  grid minor;  
xlabel('Year'); ylabel('3-yr filtered (m)')
text(1917,0.13,strcat('Max;',' R^2 = ' ,num2str(rx,2),'; S = ' ,num2str(skx,2)),'color','red','FontSize',18 )
text(1920,-0.135,strcat('Min;',' R^2 = ' ,num2str(rn,2),'; S = ' ,num2str(skn,2)),'color','blue','FontSize',18 )
text(2024.75,0.05,'Mean;','color','magenta','FontSize',18 )
text(2023.25,0.015,strcat(' R^2 = ' ,num2str(rm,2)),'color','magenta','FontSize',18 )
text(2023.25,-0.01,strcat(' S = ' ,num2str(skm,2)),'color','magenta','FontSize',18 )
set(gca,'FontSize',18)
text(1980.25,0.135,'Sydney','color','black','FontSize',22,'FontWeight','bold' )
text(2015,-.13,'c)','FontSize',18,'FontWeight','bold')

P=per;   save amplitudes_sydney P amp5m amp5n amp5x amp3m amp3n amp3x
