%% Magnitude - Frequency Analysis of Sediment Transport - Grain Size and
% Critical Shear Analysis
% Sediment Tansport Calculations Using MPM Bedload Relation 
% with Log-Normal Daily Flow PDF
%
% Note that the cvX vector must be changed to differeing lengths for
% particular figures.
%
% Joel Sholtes, PhD Student
% Dept. of Civil and Environmental Engineering
% July 25, 2013
% Updated October 31, 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)
% Med. sand, coarse sand, med. gravel, v. coarse gravel
d50 = [0.001 0.008 0.016 0.032 0.064];             
% Critical Shield's Parameter at bankfull Q for gravels - User Defined
tau_star_c = 0.03; 
% 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.5:.5:5; % To plot the Effectiveness Curves
cvX = 0.5:.1:6; % To plot the RI and Qsp Figures
% Std Dev of daily flows (cms)
stdX = meanX*cvX;       

% 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),5000);
for i = 1:length(cvX)                        
        Q(i,:) = logspace(-8,log10(Qmax(i)),5000);
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)
    % set max index to index of max flow value in each row
    [val index] = max(Q(i,:));                          
    for j = 1:index;
        % Frequency of daily discharge value
        freq_y(i,j) = lognpdf(Q(i,j),meanY(i),stdY(i)); 
        % LN2 CDF
        CDF(i,j) = logncdf(Q(i,j),meanY(i),stdY(i));
        % Prob of Exceedence
        PofE(i,j) = 1 - CDF(i,j);                       
    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 = 0.2;   
%  Exponent for in-channel at-a-station depth HG relation
b1 = 0.4;
%  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)
        % Depth (m) in-channel flow at-a-station HG
        h(i,j) = a1*Q(i,j)^b1;   
    end
end

%% Meyer-Peter Mueller Sediment Transport Relation
for i = 1:length(d50)
    % Dimensionless Shear Stress
    tau_star(:,:,i) = h.*S./(1.65.*d50(i));
end

% Qs in m^3/s
for i = 1:length(cvX)
   for k = 1:length(d50)
       QsvMPM(i,:,k) = W.*sqrt(1.65.*9.81.*d50(k).^3)*8.*...
                    (tau_star(i,:,k) - tau_star_c).^1.5;
    end
end

%% Parker (1979) Sediment Transport Relation
for i = 1:length(d50)
    % Dimensionless Shear Stress
    tau_star(:,:,i) = h.*S./(1.65.*d50(i));
end

% Qs in m^3/s
for i = 1:length(cvX)
    for k = 1:length(d50)
            QsvP(i,:,k) = W*sqrt( 1.65*9.81*d50(k) ) * d50(k) * 11.2 * ...
            (tau_star(i,:,k) - tau_star_c).^4.5 ./ tau_star(i,:,k).^3;
    end
end

% Qs in tonnes/day
Qs = Qsv.*2.65.*86400;

%% Calculate Qeff 
% Effectiveness Curves 
Eff = zeros(length(cvX),size(Q,2),length(d50));
for i = 1:length(cvX)
    for k = 1:length(d50)
        % Broken power law Qeff curve
        Eff(i,:,k) = freq_y(i,:).*Qs(i,:,k);
    end
end

% Amount of sediment transported under Eff curves.
for i = 1:length(cvX)
    for k = 1:length(d50)
        [Eff_val(i,k) Eff_index(i,k)] = max(Eff(i,:,k));
        Qeff(i,k) = Q(i,Eff_index(i,k));
        Qs_tot(i,k) = trapz(Q(i,:),Eff(i,:,k));
        ckQ = (Q(i,:) > Qeff(i,k));
        if Eff_val(i,k) < 1E-10
            Qs_tot_abv(i,k) = 0;
            Qsp(i,k) = 0;
        else
            Qs_tot_abv(i,k) = trapz(Q(i,ckQ),Eff(i,ckQ,k));
            Qsp(i,k) = Qs_tot_abv(i,k)/Qs_tot(i,k);
        end
    end
end

% Calculate the RI of Qeff (years)
for i = 1:length(cvX)
    for k = 1:length(d50)
        RI_Qeff(i,k) = 1/365/(1-normcdf((log(Qeff(i,k))-meanY(i))/stdY(i)));
    end
end

%% Calculate Q_IQR
% 1: CDF of Effectiveness Curve
Cum_Eff = zeros(length(cvX),size(Q,2),length(d50));
Frct_Cum_Eff = zeros(length(cvX),size(Q,2),length(d50));
for i = 1:length(cvX)
    for k = 1:length(d50)
        Cum_Eff(i,:,k) = cumsum(Eff(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(d50)
            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(d50)
            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,j);
            spread = Q_75(i,j)-Q_25(i,j);
            Qeff_Q50 = Qeff(i,j)/Q_50(i,j);
            Qbf_Qeff = Q_bf/Qeff(i,j);
        end
end

%% Plot Results in Color
color = flipud(summer(length(d50)+3));
lines = cellstr(char('--',':','-.','--','-',':','-.','--','-',':','--',':','-.','--','-.','-','-.','--','-',':'));
markers = cellstr(char('o','*','s','^','d','h','>','.','*','+','x','*'));

% Plot Effectiveness Curves
figure
c = [210/255 240/255 90/255;210/255 220/255 90/255;160/255 200/255 90/255;110/255 170/255 80/255;65/255 140/255 80/255];
for k = [1 3 5];
    j = 1;
for i = [1 6 16 26 36] % the step is calibrated for step in cvX of 1
    subplot(1,5,k)
semilogy(Q(i,:),Eff(i,:,k),[lines{j}],'color',color(j+2,:),'linewidth',2)
hold on
axis([0 50 .01 1000])
title(['D_{50} = ' num2str(d50(k)*1000) ' mm'])
j = j + 1;
end
end
subplot(1,5,1)
    ylabel('Sediment Transport [tonnes/day]','fontsize',12)
subplot(1,5,3)
    xlabel('Discharge [m^3/s]','fontsize',12)
subplot(1,5,5)
    legend('\eta = 0.5','\eta = 1','\eta = 2','\eta = 3','\eta = 4',...
        '\eta = 5','location','eastoutside')


% Plot Rating Curves
figure
% Semilog Axes
% for i = 1:length(d50)
%     semilogx(Q(length(cvX),:),Qs(length(cvX),:,i),[lines{i}],'color',color(i+5,:),...
%          'LineWidth',1.5)
%     hold on
% end
% axis([.01 1000 0 3E4])
% text(.011,2.1E4,'\bf Bed Material D_{50}')
for i = 1:length(d50)
    loglog(Q(length(cvX),:),Qs(length(cvX),:,i),[lines{i}],'color',color(i+2,:),...
         'LineWidth',2)
    hold on
end
text(1.5E-4,4E7,'Bed Material D_{50}')
%eqtext = '$$Q_s\propto\sqrt{(G-1)gd_{50}^3}(\tau^*-\tau^*_c)^{3/2}$$';
%text(10^-2, 10^-3, eqtext, 'Interpreter', 'Latex', 'FontSize', 12, 'Color', 'k');

axis([.0001 10^5 0.0001 10^8])
h_leg = legend('    1 mm','    8 mm','  16 mm','  32 mm','  64 mm',...
               'Location','northwest');
set(h_leg,'box','off');
xlabel('Discharge [m^3/s]','FontSize',12)
ylabel('Sediment Transport [tonnes/day]','FontSize',12)


%% Plot the Value and RI of Qeff vs. cvX
figure
c1 = [153/255 223/255 232/255;0 204/255 204/255;0 153/255 204/255;0 153/255 153/255;0 102/255 120/255];
% Qeff vs. cvX
subplot(1,2,1)
for k = 1:length(d50)
    semilogy(cvX,Qeff(:,k),[lines{k}],...
                          'Color',c1(k,:),...
                          'LineWidth',2)
    hold on
end
xlabel('Daily Flow Coefficient of Variation, \eta','FontSize',12)
ylabel('Q_{EFF} [m^3/s]','FontSize',12)
axis square
axis([0.5 5 10^-2 10^2])
text(4.5,10^10,'Bed Material D_{50}')

% Plot the RI of Qeff vs. cvX
subplot(1,2,2)
for k = 1:length(d50)
    semilogy(cvX,RI_Qeff(:,k),[lines{k}],...
                          'Color',c1(k,:),...
                          'LineWidth',2)
    hold on
end
xlabel('Daily Flow Coefficient of Variation, \eta','FontSize',12)
ylabel('Return Interval of Q_{EFF} [yr]','FontSize',12)
axis square
axis([0.5 5 10^-4 10^10])
text(4.5,10^10,'Bed Material D_{50}')
h_leg = legend('    1 mm','   8 mm','  16 mm','  32 mm','  64 mm',...
               'Location','northeast');
set(h_leg,'box','off');


%% Plot Qs_+ vs cvX
figure
cvX_plot = .5:0.01:5;
% Vogel et al 2003's fraction of sediment transported above Qeff function
f_LN = 1 - normcdf(-1.*sqrt(log(1 + cvX_plot.^2)));
p(1) = plot(cvX_plot,f_LN,'k-','LineWidth',2);
hold on
k = 2;
for i = [1 3 4 5]
    for j = 1:3:size(Qsp,1)
    p(k) = plot(cvX(j),Qsp(j,i),[markers{i}],...
                          'markersize',7,...
                          'Color',c1(i,:),...
                          'LineWidth',1.5);
    k = k + 1;
    end
end
axis square
xlabel('Daily Flow Coefficient of Variation, \eta','FontSize',12)
ylabel('Fraction of Sediment Load Transported Above Q_{EFF}','FontSize',12)
text(4,.73,'Bed Material D_{50}')
h_leg = legend([p(2) p(22) p(40) p(75)],...
               '    1 mm','  16 mm','  32 mm','  64 mm',...
               'Location','southeast');
set(h_leg,'box','off');
str = {'Percent of Sediment Transported above Q_{EFF} vs.',...
    'Daily Flow Coefficient of Variation'};
title(str,'FontSize',12)
axis([0.5 5 0.6 1])

%% Plot Q_IQR vs. CvX and d50
figure
k = 2;
hold on
for i = 1:length(d50)
    for j = 1:3:size(Q_IQR,1)
    p(k) = plot(cvX(j),Q_IQR(j,i),[markers{i}],...
                          'markersize',7,...
                          'Color',c1(i,:),...
                          'LineWidth',1.5);
    k = k + 1;
    end
end
box on
axis square
xlabel('Daily Flow Coefficient of Variation, \itC_v','FontSize',12)
ylabel('Interquartile Effectiveness Range','FontSize',12)
text(4,.32,'Bed Material D_{50}')
h_leg = legend([p(2) p(22) p(40) p(75) p(90)],...
               '  1 mm','  8 mm','16 mm','32 mm','64 mm',...
               'Location','northwest');
set(h_leg,'box','off');
axis([0.5 5 0 3])