%% Creates histogram sample figures (random)

cd('U:\Research\Qeff_Paper\MatLab\Flow Regime Analysis')

clear all; close all; 

% Flow regime assumptions
meanX = 2; % Mean flow in m3/s
cvX1 = 1:3;   % Coefficient of variations of flows... Olden and Poff 

% Sediment transport assumptions
alpha = 0.02; %arbitrary value, does not affect this analysis
beta = 1:4; %range of betas to consider (from Nash 1994)
    
% Std Dev of daily flows (cms)
stdX1 = meanX.*cvX1;

% Log normal parameters per Salas (yellow binder).. Method of Moments
meanY1 = (1/2)*log((meanX.^2)./(1+cvX1.^2));
varY1 = log(1 + cvX1.^2);
stdY1 = varY1.^0.5;

% create flow histograms for different Cv vals for 100 yr FRL
FRL = 100;
bin = 50;
densCV(length(cvX1),bin+1) = 0;
midsCV(length(cvX1),bin+1) = 0;
for i = 1:length(cvX1)
    % Probability of daily flows from standard uniform distribution
    ProbDaily = rand(FRL,365);

    % Flow values by finding inverse of log normal distribution 
    FlowDaily = logninv(ProbDaily,meanY1(i),stdY1(i));

    % Determine histogram of FlowDaily
    [Freq Mid] = hist(FlowDaily(:),bin);
    bin_width = Mid(1,2) - Mid(1,1);
    Freq(1,bin+1) = 0;
    Mid(1,bin+1) = Mid(1,bin) + bin_width;
    Density = Freq ./ sum(Freq(:)) ./ bin_width;
       
    % save midpoints and density vals in matrix
    densCV(i,:) = Density;
    midsCV(i,:) = Mid;
    
    areaCV(i) = trapz(midsCV(i,:), densCV(i,:));
end

% creates flow histograms for different FRL for Cv = 1
FRL = [10 100 100000]; % Assumed flow record length (FRL) in yrs
densFRL(length(FRL),bin+1) = 0;
midsFRL(length(FRL),bin+1) = 0;

cvX2 = 1;
% Std Dev of daily flows (cms)
stdX2 = meanX*cvX2;
% Log normal parameters per Salas (yellow binder).. Method of Moments
meanY2 = (1/2)*log((meanX^2)/(1+cvX2^2));
varY2 = log(1 + cvX2^2);
stdY2 = varY2^0.5;

for i = 1:length(FRL)
    % Probability of daily flows from standard uniform distribution
    ProbDaily = rand(FRL(i),365);

    % Flow values by finding inverse of log normal distribution 
    FlowDaily = logninv(ProbDaily,meanY2,stdY2);

    % Determine histogram of FlowDaily
    [Freq Mid] = hist(FlowDaily(:),bin);
    bin_width = Mid(1,2) - Mid(1,1);
    Freq(1,bin+1) = 0;
    Mid(1,bin+1) = Mid(1,bin) + bin_width;
    Density = Freq ./ sum(Freq(:)) ./ bin_width;
       
    % save midpoints and density vals in matrix
    densFRL(i,:) = Density;
    midsFRL(i,:) = Mid;
    
    areaFRL(i) = trapz(midsFRL(i,:), densFRL(i,:));
end

%% Plot Effectiveness Histograms with random Mid_Trials (discharge) bin values
p = 1;
figure
for j = 1:length(cvX1)
    for i = 1:length(beta)  
        eff = alpha.*midsCV(j,:).^beta(i).*densCV(j,:);
        Q = linspace(0, max(midsCV(j,:)), 200);
        
        subplot(length(cvX1)+length(FRL),length(beta),p)
        hbar = bar(midsCV(j,:),eff);
        hold on
        plot(Q, lognpdf( Q, meanY1(j), stdY1(j) ).*alpha.*Q.^beta(i),...
            'k-', 'LineWidth', 1.5);
        set(hbar,'FaceColor', [0 102/255 204/255],...
                 'EdgeColor', [0 102/255 204/255])
        set(gca, ...
              'Box','off',...
              'Layer','top',...
              'TickDir','out',...
              'TickLength',[.02 .02],...
              'LineWidth',1,...
              'FontSize',8)
        p = p + 1;
    end
end

% Add beta vals as titles to top row of histograms
for i = 1:length(beta)
    subplot(length(cvX1)+length(FRL),length(beta),i)
    title(['\bf\beta = ' num2str(beta(i))],'FontSize',10)
end

k = 1;
for i = 1:length(cvX1)
    subplot(length(cvX1)+length(FRL),length(beta),k)
    ylabel(['\bf\it Cv = \rm\bf' num2str(i)], 'FontSize',10)
    k = k + 4; % Plots only lefthandmost ylabels
end

m = p;
for j = 1:length(FRL)
    for i = 1:length(beta)  
        eff = alpha.*midsFRL(j,:).^beta(i).*densFRL(j,:);
        Q = linspace(0, max(midsFRL(j,:)), 200);
        
        subplot(length(cvX1)+length(FRL),length(beta),p)
        hbar = bar(midsFRL(j,:),eff);
        hold on
        plot(Q, lognpdf( Q, meanY2, stdY2 ).*alpha.*Q.^beta(i),...
            'k-', 'LineWidth', 1.5);
        set(hbar,'FaceColor', [0 102/255 204/255], ...
                 'EdgeColor', [0 102/255 204/255])
        set(gca, ...
              'Box','off',...
              'Layer','top',...
              'TickDir','out',...
              'TickLength',[.02 .02],...
              'LineWidth',1,...
              'FontSize',8)
        p = p + 1;
    end
end

for i = 1:length(FRL)
    subplot(length(FRL) + length(cvX1),length(beta),m)
    ylabel(['\bf' num2str(FRL(i)) ' yr'], 'FontSize',10)
    m = m + 4; % Plots only left-hand-most ylabels
end

%subplot(length(cvX_plot),length(beta_plot),7)
text(.01, .01, '\bf Sediment Discharge Density [tonnes/day/m^{3}/s]',...
    'rotation', 90, 'FontSize',11)

subplot(length(cvX1)+length(FRL),length(beta),23)
    xlabel('\bf Discharge [m^3/s]', 'FontSize',11)

