%% Qs+ discharge analysis for Werbylo and Sholtes paper
% This script samples from distributions of lognormal daily flow pdfs to
% determine f+ for different flow record lengths and CVs
% Updated May 23, 2013

clc
clear all

cd('U:\Research\Qeff_Paper\MatLab')

save_flag = 1; % 1 - save, 0 - don't

% # trials of experiment
N = 100;

% Flow regime assumptions
meanX = 2; % Mean flow in m3/s
cvX = 1:0.2:4; % Coefficient of variations of flows... Olden and Poff 
FRL = 1000000; % Assumed flow record length (FRL) in yrs

% Sediment transport assumptions
alpha = 0.02; %arbitrary value, does not affect this analysis
beta = 0.5:0.5:5; %range of betas to consider (from Nash 1994)

% Start analysis
for x = 1:size(cvX,2)
    
    % Calculate standard deviation of flows from assumptions
    stdX = cvX(x)*meanX;
    
    % Log normal parameters
    meanY = (1/2)*log((meanX^2)/(1+cvX(x)^2));
    varY = log(1 + cvX(x)^2);
    stdY = varY^0.5;
    
    % Calculate fraction above for log-normal (from Vogel eq 15)
    f_LN(1,1,x) = 1 - normcdf(-1*sqrt(log(1 + cvX(x)^2)));
    
    % Create synthetic flow record and calculate Qs+ from data
    for n = 1:N
        % Probability of daily flows from standard uniform distribution
        ProbDaily = rand(FRL,365);
        
        % Flow values by finding inverse of log normal distribution 
        FlowDaily = logninv(ProbDaily,meanY,stdY);
        
        % Determine histogram of FlowDaily
        [Freq Mid] = hist(FlowDaily(:),50);
        Freq(1,51) = 0;
        Mid(1,51) = Mid(1,50) + (Mid(1,50)-Mid(1,49));
        Freq = Freq/sum(Freq(:));
        
        % Calculate Qe and Qs+ for the range of beta values
        for i = 1:size(beta,2)
            % Calculate sediment transport using mid points of bins
            Qs(i,:) = alpha.*Mid.^(beta(i));
            
            % Calculate effective curve using mid points and fQ
            Eff(i,:) = Qs(i,:).*Freq;
                 
            % Calculate effective discharge as bin with max effective curve
            [~,c(n,i,x)] = find(Eff(i,:) == max(Eff(i,:)));
            Qe(n,i,x) = Mid(c(n,i,x));

            % Calculate daily sediment transport using flow data
            Qs_new = alpha.*FlowDaily.^(beta(i));
            
            % Check for flows above effective discharge
            ckQ = (FlowDaily > Qe(n,i,x));
            
            % Calculate Qs+ using daily flow data
            Qs_plus(n,i,x) = sum(Qs_new(ckQ))/sum(Qs_new(:));
        
        end
        
        Mid_trials(n,:,x) = Mid;
        Freq_trials(n,:,x) = Freq;
   
    end
   
end

fname = [num2str(FRL),'_yr_RL'];
AnalysisResults.Qs_plus = Qs_plus;
AnalysisResults.Qe = Qe;
AnalysisResults.Mid_trials = Mid_trials;
AnalysisResults.Freq_trials = Freq_trials;
AnalysisResults.meanX = meanX;
AnalysisResults.cvX = cvX;
AnalysisResults.FRL = FRL;
AnalysisResults.alpha = alpha;
AnalysisResults.beta = beta;

cd('U:\Research\Qeff_Paper\MatLab\Record Length Results\Run1_50bins')
if save_flag == 1
    save(fname,'AnalysisResults')
end

cd('U:\Research\Qeff_Paper\MatLab')