%% License

% Matlab code that accompanies the paper: "Towards a conceptualization of the hydrological processes behind changes
% of young water fraction with elevation: a focus on mountainous alpine catchments"
% submitted to Hydroloy and Earth System Science (HESS) 
% Written by Alessio Gentile, Politecnico and Università degli Studi di
% Torino, alessio.gentile@polito.it
% 
% Creative Commons Attribution Share Alike 3.0 license
% Original author: bettina.schaefli@giub.unibe.ch
% 
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
% 
% * Redistributions of source code must retain the above copyright
%   notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
%   notice, this list of conditions and the following disclaimer in
%   the documentation and/or other materials provided with the distribution
% 
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
%%
clc
clear 
clearvars
close all
%% Set-up

ID = 'VdN'; %Catchment ID

numiter = 1000; %Set of random parameters A_Q,A_P,Phi_Q,Phi_P for evaluating the uncertainty of alpha
numiter_str = num2str(numiter); % for the output name

wrt_rslt = 1; %1: write results 0: do not write results

%% Import data 
%Streamwater and precipitation isotopes are stored in two different files.
%Both files must be called as "ID".txt. For example: ERL.txt
%These files must be formatted as those reported in the folders:
%isotopes_streamflow.zip and isotopes_precipitation.zip downloadable from https://zenodo.org/record/4057967#.ZBlwVHbMKUk

%Reference: Staudinger, M., Seeger, S., Herbstritt, B., Stoelzle, M., Seibert, J., Stahl, K., and Weiler, M.:
%The CH-IRP data set: a decade of fortnightly data on ?2H and ?18O in streamflow and precipitation in Switzerland,
%Earth System Science Data, 12, 3057–3066, https://doi.org/10.5194/essd-12-3057-2020, 2020.

IsoQ = importIsoQ(strcat('C:\Users\aless\Desktop\Anonymous_referee1\isotopes_streamflow_filtered\',ID,'.txt')); %Table with Streamwater isotopes
IsoP = importIsoP(strcat('C:\Users\aless\Desktop\Anonymous_referee1\Precipitation_filtered\',ID,'.txt')); %Table with Precipitation isotopes

%"ID"_wQ.txt: Discharge (mm/d) at the sampling date. For example: ERL_wQ.txt
%This file must be formatted with a single column containing discharge (mm/d) at the sampling date.
%Thus, length(IsoQ.date) = length(ww_Q). It will be used for the flow-weighted fit of streamwater isotopes.

ww_Q = importQ(strcat('C:\Users\aless\Desktop\Anonymous_referee1\Weight_Q\',ID,'_wQ.txt')); 

%"ID"_wP.txt: Precipitation Volume (mm) between two sampling dates. For example: ERL_wP.txt
%This file must be formatted with a single column containing precipitation volume (mm) between two sampling dates.
%If no information about the precipitation volume of the first sampling date, i.e., ww_P(1), are available,
%set-up ww_P(1) as the precipitation volume fallen during the previous 30 days.  
%Thus, length(IsoP.date) = length(ww_P). It will be used for the volume-weighted fit of precipitation isotopes.

ww_P = importP(strcat('C:\Users\aless\Desktop\Anonymous_referee1\Weight_P_InterVolume\',ID,'_wP.txt')); 

%% Preliminary operations 

%%%%%%%%%%%%%%%%%%%%% Calculate decimal years %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%leap years in the studied period
leap_years = [2008;2012;2016;2020];

%leap years positions
idx_leapyears_Q = year(IsoQ.date) == leap_years(1) |...
                  year(IsoQ.date) == leap_years(2) |...
                  year(IsoQ.date) == leap_years(3) |...
                  year(IsoQ.date) == leap_years(4);

idx_leapyears_P = year(IsoP.date) == leap_years(1) |...
                  year(IsoP.date) == leap_years(2) |...
                  year(IsoP.date) == leap_years(3) |...
                  year(IsoP.date) == leap_years(4);
              
%Day of Year of stre   amwater isotopes sampling date           
DOY_Q = day(IsoQ.date,'dayofyear'); 

%Day of Year of precipitation isotopes sampling date 
DOY_P = day(IsoP.date,'dayofyear'); 

%initialize the number of total days 
Tot_days_Q = zeros(length(IsoQ.date),1); 
Tot_days_P = zeros(length(IsoP.date),1); 

% Define the number of total days. Leap year: 366 Not-Leap year: 365
Tot_days_Q(idx_leapyears_Q == 1) = 366;
Tot_days_Q(idx_leapyears_Q == 0) = 365;

Tot_days_P(idx_leapyears_P == 1) = 366;
Tot_days_P(idx_leapyears_P == 0) = 365;

%Calculate decimal years
dy_Q = (DOY_Q./Tot_days_Q) ; 
dy_P = (DOY_P./Tot_days_P) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%Put streamwater and precipitation isotopic composition in two different
%variables

%Isotopic composition d18O(‰) of streamwater samples
C_Q = IsoQ.d18O_mean; 

%Isotopic composition d18O(‰) of Precipitation samples
C_P = IsoP.d18O_prcp;
C_P(C_P == -9999) = NaN;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% remove NaN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

idxNaN_Q= isnan(C_Q);

C_Q(idxNaN_Q) = [];
dy_Q(idxNaN_Q) = [];
ww_Q(idxNaN_Q) = [];

idxNaN_P= isnan(C_P);

C_P(idxNaN_P) = [];
dy_P(idxNaN_P) = [];
ww_P(idxNaN_P) = [];

%% Fitting data using IRLS technique

%IRLS.m script is available from Kirchner, J. W. and Knapp, J. L. A.:
%Technical note: Calculation scripts for ensemble hydrograph separation, 
%Hydrology and Earth System Sciences, 24, 5539–5558, https://doi.org/10.5194/hess-24-5539-2020, 2020.

%IRLS_fit.m script is based on IRLS.m script. Some changes have been applied to adapt
%the code at this specific application. Where changes have been applied,
%the original code as been commented (please search in IRLS_fit.m with keywords "original
%code", where modifications were made). Changes have been made following
%the R code provided by the anonymous referee 1.

[Sin_Q,gof_Q] = IRLS_fit(C_Q,dy_Q,ww_Q,'Cauchy');
[Sin_P,gof_P] = IRLS_fit(C_P,dy_P,ww_P,'Cauchy');


%% Calculate standard error of A,Phi,k parameters

CB = 0.95; % 95% confidence

ci_Q = confint(Sin_Q, CB);              %Confidence Interval  
t_Q = tinv((1+CB)/2, gof_Q.dfe);        %Student's inverse cumulative distribution function 
SE_Q= (ci_Q(2,:)-ci_Q(1,:))./ (2*t_Q);  %Standard Error

ci_P = confint(Sin_P, CB);              %Confidence Interval  
t_P = tinv((1+CB)/2, gof_P.dfe);        %Student's inverse cumulative distribution function 
SE_P= (ci_P(2,:)-ci_P(1,:)) ./ (2*t_P); %Standard Error

%% Calculate Fyw and SE

% Following Kirchner (2016a), calculate flow-weighted young water fraction 
% as the amplitude ratio As/Ap

Fyw = abs(Sin_Q.A)/abs(Sin_P.A); 

% Standard error of F*yw
SE_Fyw= (1/abs(Sin_P.A))*sqrt((SE_Q(1)^2)+((Fyw^2)*SE_P(1)^2)); 

%% Calculate alpha and threshold age

rng(1) %for reproducibility

% Define a range of shape factors
alpha = 0.01:0.01:20; %this range is the same suggested in the R code provided by the anonymous referee 1

%Gaussian error propagation for alpha and tauyw
AmplitudeQ = Sin_Q.A;
AmplitudeP = Sin_P.A;

Phi_Q_ransampl = normrnd(Sin_Q.Fi,SE_Q(1,2),numiter,1); %generate 'numiter' Phi_Q parameters according to a Gaussian distribution
Phi_P_ransampl = normrnd(Sin_P.Fi,SE_P(1,2),numiter,1); %generate 'numiter' Phi_P parameters according to a Gaussian distribution
A_Q_ransampl = normrnd(Sin_Q.A,SE_Q(1,1),numiter,1);    %generate 'numiter' A_Q parameters according to a Gaussian distribution
A_P_ransampl = normrnd(Sin_P.A,SE_P(1,1),numiter,1);    %generate 'numiter' A_P parameters according to a Gaussian distribution


for i = 1:numiter % i: index of rows, it identifies a parameters set
    
        for j = 1:length(alpha) % j: index of columns, it identifies the alpha value
            
            % Define the error function according to Eq. (11) of Kirchner (2016a)
            % This is the error obtained for a certain set of parameters
            % (rows)with varying alpha (columns)
            
            if AmplitudeQ < 0 && AmplitudeP > 0     
            error(i,j) = abs((Phi_Q_ransampl(i)+ pi)  - (Phi_P_ransampl(i)) - (alpha(j).*atan(sqrt((abs(A_Q_ransampl(i)./A_P_ransampl(i)).^(-2./alpha(j)))-1))));
            
            elseif AmplitudeQ > 0 && AmplitudeP < 0  
            error(i,j) = abs((Phi_Q_ransampl(i))  - (Phi_P_ransampl(i)+ pi) - (alpha(j).*atan(sqrt((abs(A_Q_ransampl(i)./A_P_ransampl(i)).^(-2./alpha(j)))-1))));
            
            elseif AmplitudeQ < 0 && AmplitudeP < 0 
            error(i,j) = abs((Phi_Q_ransampl(i)+pi)  - (Phi_P_ransampl(i)+ pi) - (alpha(j).*atan(sqrt((abs(A_Q_ransampl(i)./A_P_ransampl(i)).^(-2./alpha(j)))-1))));
            
            elseif AmplitudeQ > 0 && AmplitudeP > 0             
            error(i,j) = abs((Phi_Q_ransampl(i)) - (Phi_P_ransampl(i)) - (alpha(j).*atan(sqrt((abs(A_Q_ransampl(i)./A_P_ransampl(i)).^(-2./alpha(j)))-1))));
            end
            
        end
        
  % Given a set of parameters, find the shape factor alpha for which 'error' is minimized
    alpha_opt(i,1) = alpha(error(i,:) == min(error(i,:)));
       
end
   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %Calculate alpha best using the best-fitting parameters
    if AmplitudeQ < 0 && AmplitudeP > 0     
    error_best = abs((Sin_Q.Fi + pi)  - (Sin_P.Fi) - (alpha.*atan(sqrt((abs(Sin_Q.A/Sin_P.A).^(-2./alpha))-1))));
    
    elseif AmplitudeQ > 0 && AmplitudeP < 0
    error_best = abs((Sin_Q.Fi)  - (Sin_P.Fi + pi) - (alpha.*atan(sqrt((abs(Sin_Q.A/Sin_P.A).^(-2./alpha))-1))));
    
    elseif AmplitudeQ < 0 && AmplitudeP < 0
    error_best = abs((Sin_Q.Fi + pi)  - (Sin_P.Fi + pi) - (alpha.*atan(sqrt((abs(Sin_Q.A/Sin_P.A).^(-2./alpha))-1))));
    
    elseif AmplitudeQ > 0 && AmplitudeP > 0    
    error_best = abs((Sin_Q.Fi) - (Sin_P.Fi) - (alpha.*atan(sqrt((abs(Sin_Q.A/Sin_P.A).^(-2./alpha))-1))));
    end
    
    alpha_best = alpha(error_best == min(error_best));
        
    %Calculate alpha standard deviation
    
    SD_alpha_best = std(alpha_opt);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %Calculate Thr. age (years) using alpha_best according to Eq. (14) of Kirchner (2016a)
    tauyw_Kirch_best = 0.0949 + (0.1065.*alpha_best) - (0.0126.*(alpha_best).^2);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Calculate Thr. age (years) using alpha_opt according to Eq. (14) of Kirchner (2016a)
    tauyw_Kirch_opt = 0.0949 + (0.1065.*alpha_opt) - (0.0126.*(alpha_opt).^2); 
   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %Calculate Thr. age (years) standard deviation
    SD_tauyw_Kirch_best = std(tauyw_Kirch_opt);
   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
%% Write results

Catch = cellstr(ID);

A_Q = abs(Sin_Q.A);
SE_A_Q = SE_Q(1,1);
A_P = abs(Sin_P.A);
SE_A_P = SE_P(1,1);

if AmplitudeQ < 0 && AmplitudeP > 0 
   Phi_Q = Sin_Q.Fi + pi;
   Phi_P = Sin_P.Fi;
   
elseif AmplitudeQ > 0 && AmplitudeP < 0
   Phi_Q = Sin_Q.Fi;
   Phi_P = Sin_P.Fi + pi;
   
elseif AmplitudeQ < 0 && AmplitudeP < 0
   Phi_Q = Sin_Q.Fi + pi;
   Phi_P = Sin_P.Fi + pi;
   
elseif AmplitudeQ > 0 && AmplitudeP > 0
   Phi_Q = Sin_Q.Fi;
   Phi_P = Sin_P.Fi;   
end

   SE_Phi_Q = SE_Q(1,2);
   SE_Phi_P = SE_P(1,2);


%Create a table with the output variables
outtable = table(Catch,A_Q,SE_A_Q,A_P,SE_A_P,Phi_Q,SE_Phi_Q,Phi_P,SE_Phi_P,Fyw,SE_Fyw,alpha_best,SD_alpha_best,tauyw_Kirch_best,SD_tauyw_Kirch_best)

%% Figure (1)

figure(1)

X = (0:0.001:1);
Y = [-10 10];
Xmax_P = X(Sin_P(X) == max(Sin_P(X)));

f = fittype('A*sin(2*pi*x - Fi)+ K');

Sin_Q_plot = cfit(f,A_Q,Phi_Q,Sin_Q.K);
Sin_P_plot = cfit(f,A_P,Phi_P,Sin_P.K);

Xmax_Q= X(Sin_Q_plot(X) == max(Sin_Q_plot(X)));

p(1) = scatter(dy_P,C_P-Sin_P_plot.K,25,[0 0 1],'filled','MarkerEdgeColor',[0 0 0],'MarkerFacealpha',0.75);
hold on
p(2) = scatter(dy_Q,C_Q-Sin_Q_plot.K,25,[0.3010 0.7450 0.9330],'filled','MarkerEdgeColor',[0 0 0],'MarkerFacealpha',0.75);
hold on

p(3) = plot(X,(Sin_P_plot(X)-Sin_P_plot.K),'Color',[0 0 1],'Linewidth',2);
hold on
p(4) = plot(X,(Sin_Q_plot(X)-Sin_Q_plot.K),'Color',[0.3010 0.7450 0.9330],'Linewidth',2);
hold on
plot(X,repelem(0,length(X)),'-k','Linewidth',0.5)
hold on
plot([(Sin_P_plot.Fi/(2*pi)) (Sin_P_plot.Fi/(2*pi))],[Y(1) 0],'--','Color',[1 0 0],'Linewidth',1)
hold on
plot([(Sin_Q_plot.Fi/(2*pi)) (Sin_Q_plot.Fi/(2*pi))],[Y(1) 0],'--','Color',[0.9290 0.6940 0.1250],'Linewidth',1)
hold on

plot([Xmax_P  Xmax_P], [0 Sin_P_plot.A],'-','Color',[1 0 0],'Linewidth',2.5)
hold on
plot([Xmax_Q  Xmax_Q], [0 Sin_Q_plot.A],'-','Color',[0.9290 0.6940 0.1250],'Linewidth',2.5)

hold on
% text(Xmax_P-0.1,- 1, strcat('A_{P} =',num2str(round(Sin_P_plot.A,2)),'‰'),'Fontsize',10,'fontweight','bold','Color',[1 0 0]);
text(Xmax_P-0.1,- 1, 'A_{P}','Fontsize',10,'fontweight','bold','Color',[1 0 0]);

hold on
% text(Xmax_Q-0.05,-1, strcat('A_{Q} =',num2str(round(Sin_Q_plot.A,2)),'‰'),'Fontsize',10,'fontweight','bold','Color',[0.9290 0.6940 0.1250]);
text(Xmax_Q-0.05,-1, 'A_{Q}','Fontsize',10,'fontweight','bold','Color',[0.9290 0.6940 0.1250]);

hold on
% text(Sin_P_plot.Fi/(2*pi)-0.15,Y(1)+1, strcat('\phi_{P} =',num2str(round(Sin_P_plot.Fi/(2*pi),3)),'y'),'Fontsize',10,'fontweight','bold','Color',[1 0 0]);
text(Sin_P_plot.Fi/(2*pi)-0.15,Y(1)+1, '\phi_{P}','Fontsize',10,'fontweight','bold','Color',[1 0 0]);

hold on
% text(Sin_Q_plot.Fi/(2*pi)+0.025,Y(1)+1, strcat('\phi_{Q} =',num2str(round(Sin_Q_plot.Fi/(2*pi),3)),'y'),'Fontsize',10,'fontweight','bold','Color',[0.9290 0.6940 0.1250]);
text(Sin_Q_plot.Fi/(2*pi)+0.025,Y(1)+1,'\phi_{Q}','Fontsize',10,'fontweight','bold','Color',[0.9290 0.6940 0.1250]);

xlim([0 1])
ylim(Y)


xlabel('year','FontSize',11,'fontweight','bold')
ylabel('\delta^{18}O_{P/S} + k_{P/S} (‰)','FontSize',11,'fontweight','bold')

box on

legend([p(1), p(2),p(3),p(4)],'\delta^{18}O_{P}','\delta^{18}O_{Q}','A_{P}sin(2\pift - \phi_{P})','A*_{Q}sin(2\pift - \phi*_{Q})')


title(ID)


%% Write results
if wrt_rslt == 1
    %Write outtable to a .xlsx file
    writetable(outtable,strcat('C:\Users\aless\Desktop\Anonymous_referee1\#OUT_Fyw_alpha_tau\',ID,'_Fyw_alpha_tau_',numiter_str,'.xlsx'));

    %Save figure(1)
    saveas(figure(1),strcat('C:\Users\aless\Desktop\Anonymous_referee1\#OUT_Fyw_alpha_tau\',ID,'_Fyw_alpha_tau_',numiter_str,'.fig'));
else
    disp('The results have not been written')
end
