%% Magnitude - Frequency Analysis of Sediment Transport - Broken Power Law
% Sediment Tansport Calculations Using Einstein-Brown Bedload Relation 
% and Channel-Floodplain Break along with Log-Normal Daily Flow PDF
%
% Joel Sholtes, PhD Student
% Dept. of Civil and Environmental Engineering
% March 18, 2013
% Updated June 17, 2013
% Updated November 1, 2013

clear all; close all; 
format long;

%% Channel Geometry and Flow Depth Parameters

% User defined bankfull discharge (cms)
Q_bf = 10;              
% User defined median bed material size (meters)
d50 = 0.01;             
% d50 = [.001 .005 .01 .05 .1]; 
% Critical Shield's Parameter at bankfull Q for gravels - User Defined
tau_star_c = 0.047; 
% Slope from correlation relationship from Torrizo and Pitlick 2004
S = 0.022*Q_bf^-0.49; 

%% Define Q Vector & Flow Regime Parameters

% Daily mean discharge: correlation between Qb and Qmean from gages used in Torizzo and Pitlick (2004)
meanX = 0.1512*Q_bf^0.8732; 
% User defined CV of flows
cvX = 0.4:0.3:5;           
% Std Dev of daily flows (cms)
stdX = meanX*cvX;       

% Vogel et al 2003's fraction of sediment transported above Qeff function
f_LN = 1 - normcdf(-1.*sqrt(log(1 + cvX.^2)));

% Log normal parameters per Salas (yellow binder).. Method of Moments
meanY = (1/2)*log((meanX^2)./(1+cvX.^2));
varY = log(1 + cvX.^2);
stdY = varY.^0.5;

% Qmax is set to 1 - 1/10^-12 PofE of lognormal distribution  
Qmax = logninv(1-10^-16,meanY,stdY); 

% Create Q vector (abscissa of Daily Flow PDF for each cvX)
Q = zeros(length(cvX),1000);
for i = 1:length(cvX)                        % Creates vector with 1 CMS increment from Qbf to Qmax for each cvX
        Q(i,:) = logspace(-8,log10(Qmax(i)),1000);
end

% Lognormal PDF and CDF - generates flow statistics
freq_y = zeros(size(Q));
CDF = zeros(size(Q));
PofE = zeros(size(Q));
for i = 1:length(cvX)
    [val index] = max(Q(i,:));                          % set max index to index of max flow value in each row
    for j = 1:index;
        freq_y(i,j) = lognpdf(Q(i,j),meanY(i),stdY(i)); % Frequency of daily discharge value
        CDF(i,j) = logncdf(Q(i,j),meanY(i),stdY(i));    % LN2 cum distrib function
        PofE(i,j) = 1 - CDF(i,j);                       % Prob of Exceedence
    end
%   area(i) = trapz(Q(i,1:index),freq_y(i,:));
end


%% Define Sediment Transport Inputs
% Hydraulic geometry constants
%  Coefficient for at-a-station depth HG relation 
a1_1 = 0.2;   
%  Exponent for in-channel at-a-station depth HG relation
b1_1 = [0.4 0.4 0.4];  
%  Exponents for overbank at-a-station depth HG relation
b1_2 = [0.1 0.2 0.3];  
%  Coeff for at-a-station depth HG relation for overbank floods set to match h at Q_bf for in-channel flow
for i = 1:length(b1_1)
    a1_2(i) = a1_1*Q_bf^(b1_1(i)-b1_2(i));
end
%  Coeff for downstream width HG relation set to match h at Qbf for in-channel flow (Torrizo & Pitlick 2004)
a2 = 3.25;          
%  Exponent for downstream width HG relation (Torrizo and Pitlick 2004)
b2 = 0.56;          
% Calculate rectangular channel width (meters)
W = a2*Q_bf*b2;     

% Depth-Discharge Relation
for i = 1:length(cvX)
    for j = 1:size(Q,2)
        for k = 1:length(b1_1)
            % Depth (m) in-channel flow at-a-station HG
            h_ch(i,j,k) = a1_1*Q(i,j)^b1_1(k);   
            % Depth (m) for Q > Qbf (overbank flow) at-a-station HG
            h_fp(i,j,k) = a1_2(k)*Q(i,j)^b1_2(k);   
        end
    end
end
% Creates new depth vector, h (m) with overbank HG at first
h = h_fp;
% Now replaces h with h_ch values for all in-channel flows (Q < Q_bf)
for i = 1:length(cvX)
    for k = 1:length(b1_1)
        j = 1;
        while Q(i,j) < 10
            h(i,j,k) = h_ch(i,j,k);
            j = j + 1;
        end
    end
end

%% Einstein-Brown Sediment Transport Relation
% Approximation for  0.18 < tau_star < 0.52 

% Broken Power Law Sediment Transport
tau_star = h.*S/1.65/d50;
% E-B Q_sv approximation
Qsv = W*40*sqrt(1.65*9.81*d50^3).*tau_star.^3;   
% Q_s in tonnes/day
Qs = Qsv.*2650*86400;     

% Single Power Law Sediment Transport
tau_star_ch = h_ch.*S/1.65/d50;
% E-B Q_sv approximation
Qsv_ch = W*40*sqrt(1.65*9.81*d50^3).*tau_star_ch.^3;
% Q_s in tonnes/day
Qs_ch = Qsv_ch.*2650*86400;                             

%% Create Single Power Law Function for Qs based on In-channel Transport
% Case 1: based on within-channel flow using linear regression
Q_ch = Q;
logQ = log(Q_ch);               % Take logarithm of data for log-linear regression
logQs = log(Qs_ch);             % Take logarithm of data for log-linear regression
Qs_model = zeros(length(cvX),size(Q,2));

for i = 1:length(cvX)
    for j = 1:size(Q,2);
        logQfit(j) = logQ(i,j);
        logQsfit(j) = logQs(i,j);
    end
    a(i,:) = polyfit(logQfit,logQsfit, 1);
    Qs_model(i,:) = exp(a(i,2)).*Q_ch(i,:).^a(i,1);
end
% Qs_model(isnan(Qs_model))= 0;
Qs_upper = max(Qs_model(length(cvX),:));
Q_upper = max(Qmax);

%% Calculate Qeff with Broken Power Law Qs Function (BPL) & Single Power Law
% Effectiveness Curves with Broken Power Law & Single Power Law
for i = 1:length(cvX)
    for k = 1:length(b1_1)
        % Broken power law Qeff curve
        Eff1(i,:,k) = freq_y(i,:).*Qs(i,:,k); 
        % Single power law Qeff curve
        Eff2(i,:,k) = freq_y(i,:).*Qs_model(i,:); 
    end
end

% Amount of sediment transported under Eff curves.
% -> Note that while the Effectiveness curves are different from SPL to the 
%   three BPL case, the values of the effective discharge does not change
%   since it is below the bankfull discharge
for i = 1:length(cvX)
    for k = 1:length(b1_1)
        [Eff_val(i,k) Eff_index(i,k)] = max(Eff1(i,:,k));
        Qeff(i) = Q(i,Eff_index(i,k));
        Qs_tot1(i,k) = trapz(Q(i,:),Eff1(i,:,k));
        Qs_tot2(i,k) = trapz(Q(i,:),Eff2(i,:,k));
    end
end

% Calculate Fraction of Sediment Transported Above Qeff
for i = 1:length(cvX);
    for k = 1:length(b1_1)
        j = 1;
        % Calc amount of sed transported above Qeff for each curve to curtailment
        while Q(i,Eff_index(i,k)-1+j) < max(Q(i,:))
            Q_abv(i,j,k) = Q(i,Eff_index(i,k)-1+j);
            Eff_abv1(i,j,k) = Eff1(i,Eff_index(i)-1+j,k);
            Eff_abv2(i,j,k) = Eff2(i,Eff_index(i)-1+j,k);
            j = j+1;
        end
        Qs_tot_abv1(i,k) = trapz(Q_abv(i,:,k),Eff_abv1(i,:,k));
        Qs_tot_abv2(i,k) = trapz(Q_abv(i,:,k),Eff_abv2(i,:,k));

        % Percent of total sed transported above Qeff for each curve relative to single power law function
        pct_broke(i,k) = Qs_tot_abv1(i,k)/Qs_tot2(i,k);
        pct_power(i,k) = Qs_tot_abv2(i,k)/Qs_tot2(i,k);
     end
end

%% Calculate Q_IQR
% Append all effectiveness curves into one matrix
Eff_Curves(:,:,1) = Eff2(:,:,1);
for i = 1:size(Eff2,3)
    Eff_Curves(:,:,i+1) = Eff1(:,:,i);
end

% 1: CDF of Effectiveness Curve
Cum_Eff = zeros(length(cvX),size(Q,2),length(b1_1)+1);
Frct_Cum_Eff = zeros(length(cvX),size(Q,2),length(b1_1)+1);
for i = 1:length(cvX)
    for k = 1:length(b1_1)+1
        Cum_Eff(i,:,k) = cumsum(Eff_Curves(i,:,k));
        Frct_Cum_Eff(i,:,k) = Cum_Eff(i,:,k)./max(Cum_Eff(i,:,k));
    end
end

%2: ID indices and upper and lower bound values or Frct_Cum_Eff qtiles
qtile = [0.25 0.50 0.75];
% Index and actual values for upper and lower bounds of quartiles of 
% Frct_Cum_Eff for an awkward linear interpolation
for i = 1:length(cvX)
        for j = 1:length(b1_1)+1
            for k = 1:length(qtile)
            i_lower(i,j,k)  = find(Frct_Cum_Eff(i,:,j) <= qtile(k),1,'last');
            i_upper(i,j,k) = find(Frct_Cum_Eff(i,:,j) >= qtile(k),1,'first');
            lower_Q(i,j,k) = Q(i,i_lower(i,j,k));
            upper_Q(i,j,k) = Q(i,i_upper(i,j,k));
            lower_Frct(i,j,k) = Frct_Cum_Eff(i,i_lower(i,j,k),j);
            upper_Frct(i,j,k) = Frct_Cum_Eff(i,i_upper(i,j,k),j);
            end
        end
end

%3: Calc discharges assoc. with qtiles of cum sed transport and Q_IQR
% using linear interpolation
for i = 1:length(cvX)
        for j = 1:length(b1_1)+1
            Q_25(i,j) = lower_Q(i,j,1)+(upper_Q(i,j,1)-lower_Q(i,j,1))...
                 *(0.25-lower_Frct(i,j,1))/(upper_Frct(i,j,1)-...
                 lower_Frct(i,j,1));
            Q_50(i,j) = lower_Q(i,j,2)+(upper_Q(i,j,2)-lower_Q(i,j,2))...
                 *(0.50-lower_Frct(i,j,2))/(upper_Frct(i,j,2)-...
                 lower_Frct(i,j,2));
            Q_75(i,j) = lower_Q(i,j,3)+(upper_Q(i,j,3)-lower_Q(i,j,3))...
                 *(0.75-lower_Frct(i,j,3))/(upper_Frct(i,j,3)-...
                 lower_Frct(i,j,3));
            Q_IQR(i,j) = (Q_75(i,j)-Q_25(i,j))/Qeff(i);
            spread(i,j) = Q_75(i,j)-Q_25(i,j);
            Qeff_Q50(i,j) = Qeff(i)/Q_50(i,j);
            Qbf_Qeff(i,j) = Q_bf/Qeff(i);
        end
end

%% Plot Results in Gray Scale
colors = gray(length(cvX)+8);

% Plot Rating Curves
figure
subplot(1,2,1)
% BPL Curves with different overbank at-a-station HG exponents
loglog(Q(length(cvX),:),Qs(length(cvX),:,3),'linewidth',5,'color',[.8 .8 .8])
hold on
loglog(Q(length(cvX),:),Qs(length(cvX),:,2),'linewidth',5,'color',[.5 .5 .5])
loglog(Q(length(cvX),:),Qs(length(cvX),:,1),'linewidth',5,'color',[0 0 0])
% Single Power Law Sediment Rating Curve
loglog(Q(length(cvX),:),Qs_model(length(cvX),:),'color',[.6 .6 .6],...
    'linewidth',1.5,'linestyle','--')
vertical_bar = logspace(-10,15,1000);
loglog(10,vertical_bar,'color',[0 0 0],'linewidth',1,'linestyle',':')
title('Broken & Single Power Law Discharge-Sediment Transport Relations')
axis square
axis([.001 Q_upper 1 Qs_upper])
xlabel('Discharge [m^3/s]')
ylabel('Sediment Discharge [tonnes/day]')
h_leg = legend('BPL, b = 0.3','BPL, b = 0.2','BPL, b = 0.1','SPL',...
               'location','se');
set(h_leg,'box','off');
eqtext = '$$Q_s=\psi\sqrt{(G-1)gd_{50}^3}\tau_{*}^3$$';
% Add the string containing the Latex expression to the plot
text(10^2, 10^6.5, eqtext, 'Interpreter', 'Latex', 'FontSize', 12, 'Color', 'k');
% Text string for Q_BF designation
q_text = '$$Q_{BF}=10 m^3/s$$';
text(10^1.25, 10^4, q_text, 'Interpreter', 'Latex', 'FontSize', 12, 'Color', 'k');


% Plot Qs_+ vs cvX
subplot(1,2,2)
plot(cvX,f_LN,'k-','LineWidth',1.5)
hold on
% plot(cvX,pct_broke(:,3),'k--','LineWidth',1.5)
% plot(cvX,pct_broke(:,2),'k-.','LineWidth',1.5)
% plot(cvX,pct_broke(:,1),'k:','LineWidth',1.5)
plot(cvX,pct_broke(:,3),'o','MarkerEdgeColor',[.8 .8 .8],...
                        'MarkerFaceColor',[.8 .8 .8],...
                        'MarkerSize',5)
plot(cvX,pct_broke(:,2),'+','color',[.5 .5 .5],'LineWidth',1.5,'MarkerSize',7)
plot(cvX,pct_broke(:,1),'x','color',[0 0 0],'LineWidth',1.5,'MarkerSize',7)
plot(cvX,pct_power(:,1),'o','color',[.6 .6 .6],'MarkerSize',6,'LineWidth',1.5)
axis square
xlabel('Daily Flow Coefficient of Variation, \eta')
ylabel('Fraction of Sediment Load Transported Above \it{Q_{eff}}')
h_leg = legend('Analytical Value','BPL, b = 0.3','BPL, b = 0.2',...
               'BPL, b = 0.1','SPL','Analytical Value',...
               'location','southwest');
set(h_leg,'box','off');
str = {'Percent of Sediment Transported above {\itQ}_{EFF} vs.',...
    'Daily Flow Coefficient of Variation'};
title(str)
axis([0.5 5 0.6 1])

%% Plot Results in Color
% Plot Rating Curves
col = summer(5);
figure
subplot(1,2,1)
% BPL Curves with different overbank at-a-station HG exponents
loglog(Q(length(cvX),:),Qs(length(cvX),:,3),'linewidth',5,'color',[0.75,0.875,0.4])
hold on
loglog(Q(length(cvX),:),Qs(length(cvX),:,2),'linewidth',5,'color',[.5 .75 .4])
loglog(Q(length(cvX),:),Qs(length(cvX),:,1),'linewidth',5,'color',[0 .5 .4])
% Single Power Law Sediment Rating Curve
loglog(Q(length(cvX),:),Qs_model(length(cvX),:),'color',[.6 .6 .6],...
    'linewidth',1.5,'linestyle','--')
vertical_bar = logspace(-10,15,1000);
loglog(10,vertical_bar,'color',[0 0 0],'linewidth',1,'linestyle',':')
title('Broken & Single Power Law Discharge-Sediment Transport Relations')
axis square
axis([.001 Q_upper 1 Qs_upper])
xlabel('Discharge [m^3/s]')
ylabel('Sediment Discharge [tonnes/day]')
h_leg = legend('BPL, b = 0.3','BPL, b = 0.2','BPL, b = 0.1','SPL',...
               'location','se');
set(h_leg,'box','off');
% Add the string containing the Latex expression to the plot
eqtext = '$$Q_s=\psi\sqrt{(G-1)gd_{50}^3}\tau_{*}^3$$';
text(10^2, 10^6.5, eqtext, 'Interpreter', 'Latex', 'FontSize', 12, 'Color', 'k');
% Text string for Q_BF designation
q_text = '$$Q_{BF}=10 m^3/s$$';
text(10^1.25, 10^4, q_text, 'Interpreter', 'Latex', 'FontSize', 12, 'Color', 'k');

% Plot Qs_+ vs cvX
subplot(1,2,2)
plot(cvX,f_LN,'k-','LineWidth',1.5)
hold on
% plot(cvX,pct_broke(:,3),'k--','LineWidth',1.5)
% plot(cvX,pct_broke(:,2),'k-.','LineWidth',1.5)
% plot(cvX,pct_broke(:,1),'k:','LineWidth',1.5)
plot(cvX,pct_broke(:,3),'o','MarkerEdgeColor',[0.75,0.875,0.4],...
                        'MarkerFaceColor',[0.75,0.875,0.4],...
                        'MarkerSize',5)
plot(cvX,pct_broke(:,2),'+','color',[.5 .75 .4],'LineWidth',1.5,'MarkerSize',7)
plot(cvX,pct_broke(:,1),'x','color',[0 .5 .4],'LineWidth',1.5,'MarkerSize',7)
plot(cvX,pct_power(:,1),'o','color',[.6 .6 .6],'MarkerSize',6,'LineWidth',1.5)
axis square
xlabel('Daily Flow Coefficient of Variation')
ylabel('Fraction of Sediment Load Transported Above \it{Q_{eff}}')
h_leg = legend('Analytical Value','BPL, b = 0.3','BPL, b = 0.2',...
               'BPL, b = 0.1','SPL','Analytical Value',...
               'location','southwest');
set(h_leg,'box','off');
str = {'Percent of Sediment Transported above {\itQ}_{EFF} vs.',...
    'Daily Flow Coefficient of Variation'};
title(str)
axis([0.5 5 0.6 1])

%% Plot Q_IQR vs. CvX and BPL
figure
hold on
plot(cvX,Q_IQR(:,1),'o','color',[.6 .6 .6],'MarkerSize',10,'LineWidth',2)
plot(cvX,Q_IQR(:,4),'o','MarkerEdgeColor',[0.75,0.875,0.4],...
                        'MarkerFaceColor',[0.75,0.875,0.4],...
                        'MarkerSize',8)
plot(cvX,Q_IQR(:,3),'+','color',[.5 .75 .4],'LineWidth',1.5,'MarkerSize',8)
plot(cvX,Q_IQR(:,2),'x','color',[0 .5 .4],'LineWidth',1.5,'MarkerSize',8)
axis square
xlabel('Daily Flow Coefficient of Variation','FontSize',12)
ylabel('Interquartile Sediment Yield Range','FontSize',12)
h_leg = legend('SPL, b_{fp} = 0.4','BPL, b_{fp} = 0.3','BPL, b_{fp} = 0.2',...
               'BPL, b_{fp} = 0.1','location','southeast');
set(h_leg,'box','off');
box on