%% Creates figures for Qs plus analysis

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

clear all; close all; 

cd('U:\Research\Qeff_Paper\MatLab\Record Length Results\Run3_50bins')

% Flow Statistics
meanX = 2;
% User defined CV of flows
cvX1 = 0.5:0.1:5;          
cvX2 = 0.5:0.2:5.2;
cvX_plot = [6 16 26];  % cv = 1, 2, 3
% Std Dev of daily flows (cms)
stdX1 = meanX*cvX1;
stdX2 = meanX*cvX2;

% 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;
meanY2 = (1/2)*log((meanX^2)./(1+cvX2.^2));
varY2 = log(1 + cvX2.^2);
stdY2 = varY2.^0.5;
% Sediment Rating Curve Parameters
beta1 = 0.5:0.1:5;
beta2 = 0.5:0.1:5.1;
beta_plot = [1 2 3 4];

%Record Length
RL = [10 100 1000 100000]; %Record Lengths: 5, 10, 50, 100, 500, 1000, 100,000
%RL = [5 25 50 100];
% Vogel et al 2003's fraction of sediment transported above Qeff function
f_LN = 1 - normcdf(-1.*sqrt(log(1 + cvX1.^2)));  


% Load Qs_plus and Qe results for each flow record length
for x = 1:size(RL,2)
    eval(['load ',num2str(RL(x)),'_yr_RL'])
    eval(['Qs_plus_',num2str(RL(x)),'= AnalysisResults.Qs_plus;'])
    eval(['Qe_',num2str(RL(x)),'= AnalysisResults.Qe;'])
    eval(['Mid_trials_',num2str(RL(x)),'= AnalysisResults.Mid_trials;'])
    eval(['Freq_trials_',num2str(RL(x)),'= AnalysisResults.Freq_trials;'])
end

alpha = AnalysisResults.alpha;

% Load flow record data results for each flow record length and convert to 
% effectiveness histograms for illustrative purposes
% Here I want to plot a series of Eff histograms featuring (median value of
% mid_trials ? or simply one example? Try both. For each CVX, I will show
% how the Eff histogram changes with beta. Show using 100 yr and 100,000 yr
% FRL (or all?) also, show a series for one or more FRL and one or more
% beta values I will show how things change with increasing CVX.

% Medians and Standard deviations for trials of Qs_plus and Qe
for x = 1:size(RL,2)
    eval(['Median_Qs_plus_',num2str(RL(x)),' = median(Qs_plus_',num2str(RL(x)),',1);'])
%     eval(['Std_Qs_plus_',num2str(RL(x)),' = std(Qs_plus_',num2str(RL(x)),');'])
    eval(['Median_Qe_',num2str(RL(x)),' = median(Qe_',num2str(RL(x)),',1);'])
%     eval(['Std_Qe_',num2str(RL(x)),' = std(Qe_',num2str(RL(x)),');'])
end

% Convert to MxNxP to MxN Matrices with CV in the row dir and Beta in the
% column dir
%
% NOTE: for RL = 10, 100, and 1,000 we have 46x46 matrices, but for 
% RL = 100k, we have 24x24's
for x = 1:size(RL,2)
    if x < 6  % 4 for RL > 1000
        for i = 1:length(cvX1)
            eval(['Qsp_med_' num2str(RL(x)) '(' num2str(i) ',:) = Median_Qs_plus_' ...
                num2str(RL(x)) '(:,:,' num2str(i) ');'])
            eval(['Qe_med_' num2str(RL(x)) '(' num2str(i) ',:) = Median_Qe_' ...
                num2str(RL(x)) '(:,:,' num2str(i) ');'])

    %         %Create Median histograms of flow data for each CV
    %         for j = 1:size(Mid_trials_10,2)
    %             eval(['Med_mid_trials_' num2str(RL(x)) '(i,j) = median(Mid_trials_'...
    %                 num2str(RL(x)) '(:,j,i),1);'])
    %             eval(['Med_freq_trials_' num2str(RL(x)) '(i,j) = median(Freq_trials_'...
    %                 num2str(RL(x)) '(:,j,i),1);'])
        end
    else
        for i = 1:length(cvX2)
            eval(['Qsp_med_' num2str(RL(x)) '(' num2str(i) ',:) = Median_Qs_plus_' ...
                num2str(RL(x)) '(:,:,' num2str(i) ');'])
            eval(['Qe_med_' num2str(RL(x)) '(' num2str(i) ',:) = Median_Qe_' ...
                num2str(RL(x)) '(:,:,' num2str(i) ');'])
        end
        
    end
end

% Create matrices with PofE, RI of Qeff, and logRI of Qeff
for x = 1:size(RL,2)
    if x < 6% 4 for RL > 1000
        for i = 1:length(cvX1)
            for j = 1:length(beta1)
            eval(['Qe_PofE_' num2str(RL(x)) '(i,j) = ' ...
                  '1-logncdf(Qe_med_' num2str(RL(x)) '(i,j),meanY1(i),stdY1(i));'])
            eval(['Qe_RI_' num2str(RL(x)) '(i,j) = '...
                  '1/365/(1-normcdf((log(Qe_med_' num2str(RL(x)) ...
                  '(i,j))-meanY1(i))/stdY1(i)));'])
            eval(['logQe_RI_' num2str(RL(x)) '(i,j) = '... 
                  'log10(Qe_RI_' num2str(RL(x)) '(i,j));'])
            end
        end
    else
        for i = 1:length(cvX2)
            for j = 1:length(beta2)
            eval(['Qe_PofE_' num2str(RL(x)) '(i,j) = ' ...
                  '1-logncdf(Qe_med_' num2str(RL(x)) '(i,j),meanY2(i),stdY2(i));'])
            eval(['Qe_RI_' num2str(RL(x)) '(i,j) = '...
                  '1/365/(1-normcdf((log(Qe_med_' num2str(RL(x)) ...
                  '(i,j))-meanY2(i))/stdY2(i)));'])
            eval(['logQe_RI_' num2str(RL(x)) '(i,j) = '... 
                  'log10(Qe_RI_' num2str(RL(x)) '(i,j));'])
            end
        end
    end
end
% 
% % Create Median Effectiveness Histograms
% for x = 1:size(RL,2)
%     for i = 1:length(cvX)
%         for j = 1:length(beta_plot)
%             eval(['Eff_med_' num2str(RL(x)) '(i,:,j) = alpha.*Med_mid_trials_' ...
%             num2str(RL(x)) '(i,:).^beta_plot(j).*Med_freq_trials_' ...
%             num2str(RL(x)) '(i,:);'])
%         end
%     end
% end

% %% Plot Discharge CDF by FRL (Random Selection) - Try median value again!
% % Each plot is of a different flow record length with values of CVx from
% % 1:4 in it
% NEEDS WORK STILL !!!
% trial = round(100*rand(1));
% figure
% for i = 1:length(cvX_plot)
%     p = 1;
%     for j = 1:length(RL)
%         subplot(2,2,p)
%         eval(['x = Med_mid_trials_' num2str(RL(j)) ...
%               '(cvX_plot(i),:);'])
%         eval(['y = cumsum(Med_freq_trials_' num2str(RL(j)) ...
%               '(cvX_plot(i),:));'])
%         plot(x,y)
%         hold on
%         p = p + 1;
%         axis([0 50 0 1])
%         axis square
%     end
% 
% end
% 
%    
% for i = 1:length(RL)
%     subplot(length(cvX_plot),length(RL),i)
%     title(['Record Length = ' num2str(RL(i)) ' yr'])
% end
% 
% p = 1;
% for i = 1:length(cvX_plot)
%     subplot(length(cvX_plot),length(RL),p)
%     ylabel(['\it CV = \rm\bf' num2str(i)])
%     p = p + length(RL);
% end
% 
% %subplot(length(cvX_plot),length(beta_plot),7)
% text(.01, .01, '\bf Sediment Discharge [mass/day]','rotation', 90)

%% Plot Effectiveness Histograms with random Mid_Trials (discharge) bin values
FRL = 100;   %Choose flow record length
trial = round(100*rand(1));
figure
p = 1;
% need to convert relative frequency histo values to density by dividing by
% bin width
for j = 1:length(cvX_plot) % check to see if 1-16
    for i = 1:length(beta_plot)  
        subplot(length(cvX_plot),length(beta_plot),p)
        eval(['width = Mid_trials_' num2str(FRL) '(trial,2,cvX_plot(j)) - ' ...
              'Mid_trials_' num2str(FRL) '(trial,1,cvX_plot(j));'])
        eval(['eff = alpha.*Mid_trials_' num2str(FRL) '(trial,:,cvX_plot(j)).^' ...
              'beta_plot(i).*Freq_trials_' num2str(FRL) '(trial,:,cvX_plot(j)) ./ width;'])
        eval(['hbar = bar(Mid_trials_' num2str(FRL) '(trial,:,cvX_plot(j)),' ...
              'eff);'])
        eval(['Q = linspace(0, max(Mid_trials_' num2str(FRL) ...
             '(trial,:,cvX_plot(j))), 200);'])         
        hold on
        plot(Q, lognpdf( Q, meanY1(j), stdY1(j) ).*alpha.*Q.^beta_plot(i),...
            'k-', 'LineWidth', 1.5);
        set(hbar,'FaceColor', 'b', 'EdgeColor', 'b')
        set(gca, ...
              'Box','off',...
              'Layer','top',...
              'TickDir','out',...
              'TickLength',[.02 .02],...
              'LineWidth',1,...
              'FontSize',8)
        p = p + 1;
    end
end

for i = 1:length(beta_plot)
    subplot(length(cvX_plot),length(beta_plot),i)
    title(['\bf\beta = ' num2str(beta_plot(i))],'FontSize',11)
end

p = 1;
for i = 1:length(cvX_plot)
    subplot(length(cvX_plot),length(beta_plot),p)
    ylabel(['\bf\it Cv = \rm\bf' num2str(i)], 'FontSize',11)
    p = p + 4; % Plots only lefthandmost 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)
%    'VerticalAlignment','middle','HorizontalAlignment','left');
% 
% str = {['\bf Random Effectiveness Histograms for ' num2str(FRL) 'yr Flow Record']};
% text(7,220,str)
subplot(length(RL),length(beta_plot),11)
    xlabel('\bf Discharge [m^3/s]', 'FontSize',11)

%% Plot Effectiveness Histograms with random Mid_Trials (discharge) bin values
%  for CV = 1, range of Beta values, and vary FRL in the row direction
RL = [10 100 100000];
trial = round(100*rand(1));
figure
p = 1;
for j = 1:length(RL)-1
    for i = 1:length(beta_plot)  
        subplot(length(RL),length(beta_plot),p)
        eval(['eff = alpha.*Mid_trials_' num2str(RL(j)) '(trial,:,6).^' ...
              'beta_plot(i).*Freq_trials_' num2str(RL(j)) '(trial,:,6);'])
        eval(['hbar = bar(Mid_trials_' num2str(RL(j)) '(trial,:,6),' ...
              'eff);'])
        eval(['Q = linspace(0, max(Mid_trials_' num2str(RL(j)) ...
            '(trial,:,6)), 200);'])
        hold on
        plot(Q, lognpdf( Q, meanY1(cvX_plot(j)), stdY1(cvX_plot(j)) )...
            .*alpha.*Q.^beta_plot(i),'k-', 'LineWidth', 1.5);
        x_max = max(eff);
        set(hbar,'FaceColor', 'b', 'EdgeColor', 'b')
        set(gca, ...
              'Box','off',...
              'Layer','top',...
              'TickDir','out',...
              'TickLength',[.02 .02],...
              'LineWidth',1,...
              'FontSize',8)
        p = p + 1;
    end
end
% 1000 yr record length has different number of CVx values (23 instead of
% 46)
for j = length(RL)
    for i = 1:length(beta_plot)  
        subplot(length(RL),length(beta_plot),p)
        eval(['eff = alpha.*Mid_trials_' num2str(RL(j)) '(trial,:,3).^' ...
              'beta_plot(i).*Freq_trials_' num2str(RL(j)) '(trial,:,3);'])
        eval(['hbar = bar(Mid_trials_' num2str(RL(j)) '(trial,:,3),' ...
              'eff);'])
        eval(['Q = linspace(0, max(Mid_trials_' num2str(RL(j)) ...
            '(trial,:,3)), 200);'])
        hold on
        plot(Q, lognpdf( Q, meanY2(cvX_plot(j)), stdY2(cvX_plot(j)) )...
            .*alpha.*Q.^beta_plot(i), 'k-', 'LineWidth', 1.5);
        x_max = max(eff);
        set(hbar,'FaceColor', 'b', 'EdgeColor', 'b')
        set(gca, ...
              'Box','off',...
              'Layer','top',...
              'TickDir','out',...
              'TickLength',[.02 .02],...
              'LineWidth',1,...
              'FontSize',8)
        p = p + 1;
    end
end

for i = 1:length(beta_plot)
    subplot(length(RL),length(beta_plot),i)
    title(['\bf\beta = ' num2str(beta_plot(i))], 'FontSize',11)
end

p = 1;
for i = 1:length(RL)
    subplot(length(RL),length(beta_plot),p)
    ylabel(['\bf' num2str(RL(i)) ' yr'], 'FontSize',11)
    p = p + 4; % Plots only left-hand-most ylabels
end

subplot(length(RL),length(beta_plot),11)
    xlabel('\bf Discharge [m^3/s]', 'FontSize',11)

text(.01, .01, '\bf Sediment Discharge Density [tonnes/day/m^{3}/s]','rotation', 90)

%% Contour Plots for Qs_plus
% I need two contour commands: one for the 46x46 plots and one for the
% 24x24 plots
figure
%RL = [100000 1000 100 10];
RL = [100 50 25 5];
x = length(RL);
i = 1;
while x > 1 
        subplot(1,4,i)
        contourfm(beta1,cvX1,eval(['transpose(Qsp_med_' num2str(RL(x)) ')']),10,...
                  'LineColor','none')
        colormap(flipud(gray)) % Change to gray if you want BW
        title(['Record Length = ' num2str(RL(x)) ' yr'])
        axis square
        axis([1 4 0.5 5])
        set(gca,...
            'Box', 'on',...
            'Layer', 'top',...
            'TickDir', 'out',...
            'TickLength', [0.02 0.02],...
            'XTickmode'   , 'manual',...
            'YTick', [0.5:0.5:5],...
            'YTickLabel',[' ';'1';' ';'2';' ';'3';' ';'4';' ';'5'],...
            'XTick',[0.5:0.5:5],...
            'XTickLabel',[' ';'1';' ';'2';' ';'3';' ';'4';' ';'5']);
        caxis([0 1])        % Sets contour color range to be standard
        x = x - 1;
        i = i + 1;
end

% subplot(1,4,4)
% contourfm(beta2,cvX2,transpose(Qsp_med_100000),10,...
%           'LineColor','none')
% colormap(flipud(gray)) % Change to gray if you want BW
% title('Record Length = 100,000 yr')
% axis square
% axis([0.5 5 0.5 5])
%         set(gca,...
%             'Box', 'on',...
%             'Layer', 'top',...
%             'TickDir', 'out',...
%             'TickLength', [0.02 0.02],...
%             'XTickmode'   , 'manual',...
%             'YTick', [0.5:0.5:5],...
%             'YTickLabel',[' ';'1';' ';'2';' ';'3';' ';'4';' ';'5'],...
%             'XTick',[0.5:0.5:5],...
%             'XTickLabel',[' ';'1';' ';'2';' ';'3';' ';'4';' ';'5']);
%         caxis([0 1])        % Sets contour color range to be standard
        
subplot(1,4,1)
    ylabel('Sediment Rating Curve Exponent, \beta','FontSize',12)
subplot(1,4,2)
    text(2,5,'Fraction of Sediment Transported above Q_{eff} , f_+','FontSize',12)
subplot(1,4,3)
    xlabel('Daily Flow Coefficient of Variation','FontSize',12)
subplot(1,4,4)
    colorbar
    
%% Filled Contour Plots for Qe RI using Log10 intervals
figure
x = length(RL);
i = 1;
% Creates Blue Colormap 
% Otherwise, change to grey below
while x > 1 
        subplot(1,4,i)
        v = [-3:.5:6];
        contourfm(beta1,cvX1,eval(['transpose(logQe_RI_' num2str(RL(x)) ')']),v,...
                  'LineColor','none')
        colormap(flipud(gray)) % Change to gray if you want BW
        title(['Record Length = ' num2str(RL(x)) ' yr'])
        axis square
        axis([0.5 5 0.5 5])
        set(gca,...
            'Box', 'on',...
            'Layer', 'top',...
            'TickDir', 'out',...
            'TickLength', [0.02 0.02],...
            'XTickmode'   , 'manual',...
            'XTick', [0.5:0.5:5],...
            'XTickLabel',[' ';'1';' ';'2';' ';'3';' ';'4';' ';'5'],...
            'YTick',[0.5:0.5:5],...
            'YTickLabel',[' ';'1';' ';'2';' ';'3';' ';'4';' ';'5']);
 
        caxis([-3 6])        % Sets contour color range to be standard
        x = x - 1;
        i = i + 1;
end

        subplot(1,4,4)
        v = [-2.5:.5:6];
        contourfm(beta2,cvX2,transpose(logQe_RI_100000),v,...
                  'LineColor','none')
        colormap(flipud(gray)) % Change to gray if you want BW
        title('Record Length = 100,000 yr')
        axis square
        axis([0.5 5 0.5 5])
        set(gca,...
            'Box', 'on',...
            'Layer', 'top',...
            'TickDir', 'out',...
            'TickLength', [0.02 0.02],...
            'XTickmode'   , 'manual',...
            'XTick', [0.5:0.5:5],...
            'XTickLabel',[' ';'1';' ';'2';' ';'3';' ';'4';' ';'5'],...
            'YTick',[0.5:0.5:5],...
            'YTickLabel',[' ';'1';' ';'2';' ';'3';' ';'4';' ';'5']);
 
        caxis([-3 6])        % Sets contour color range to be standard
        

subplot(1,4,1)
    ylabel('Sediment Rating Curve Exponent, \beta','FontSize',12)
subplot(1,4,2)
    text(5,6,'Return interval of \itQ_{EFF} \rm [yr]','FontSize',12)
subplot(1,4,3)
    xlabel('Daily Flow Coefficient of Variation','FontSize',12)
subplot(1,4,4)
    colorbar('YTickLabel',{'1E-2','1E0','1E2','1E4','1E6'},...
             'TickDir','out')
    
% %% Plot Qs_plus vs Beta
% lines = cellstr(char('-',':','-.','--','-',':','-.','--','-',':'));
% markers = cellstr(char('s','^','d','o','v','<','>','.','*','+','x'));
% % colors = {'m','c','r','g','b','k','m','c','r','g','b'};
% colors = summer(length(RL)+2); %summer % gray
% figure
% x = 1;
% while x <= 4 
%     for z = [1 2 6 16]
%         subplot(2,2,x)
%         for i = 1:size(RL,2)
%             eval(['Qs_plot1 = Median_Qs_plus_',num2str(RL(i)),'(1,:,z);'])
%             plot(beta,Qs_plot1,[markers{i}],'Color',colors(i+1,:),...
%                 'Linestyle','-','MarkerFaceColor',colors(i+1,:),...
%                 'MarkerSize',5)
%             hold on
%         end
%     title(['CV = ',num2str(cvX(z))])
% %     xlabel('beta')
% %     ylabel('Qs plus')
%     ylim([0 1])
%     axis square
% %     legend('100000', '1000', '100', '10','Location','EastOutside')
%     x = x + 1;
%     end
% 
% end
% 
% %% Plot Qs_plus vs CV_X
% figure
% b = 1;
% while b <= 6
%     for x = [1 2 4 6 8 10] % hard coding indices of beta values to pull
%         subplot(3,2,b)
%         for i = 1:size(RL,2)
%             eval(['Qs_plot2(1,:) = Median_Qs_plus_',num2str(RL(i)),'(1,x,:);'])
%             plot(cvX,Qs_plot2,[markers{i}],'Color',colors(i+1,:),...
%                 'Linestyle','-','MarkerFaceColor',colors(i+1,:),...
%                 'MarkerSize',5)
%             hold on
%         end
%         plot(cvX,f_LN,'k-')
%         title(['\beta = ', num2str(beta(x))])
%     %     xlabel('cvX')
%     %     ylabel('Qs plus')
%         ylim([0 1])
%         axis square
% %         legend('100000', '1000', '100', '10','Location','EastOutside')
%         b = b + 1;
%     end
% end
% 
% 
% %% Plot Qs Plus vs Beta for one FRL and all CVx
% figure
% %RL = [1000000 100000 1000 100 10];
% RL = [10 100 1000 100000];
% colors = flipud(gray(length(beta)+3)); %summer % gray
% for b = 1:length(RL)
%     i = 1;
%     subplot(2,2,b)
%     for x = [1 2 4 6 8 10]
%         eval(['Qs_plot3(1,:) = Median_Qs_plus_',num2str(RL(b)),'(1,x,:);'])
%         plot(cvX,Qs_plot3,[markers{i}],'Color',colors(x+2,:),...
%             'Linestyle','-','MarkerFaceColor',colors(x+2,:),...
%             'MarkerSize',5,'LineWidth',1.5)
%         hold on
%         i = i + 1;
%     end
%     plot(cvX,f_LN,'color',[0.3 0.3 0.3],'LineWidth',1.5,'LineStyle','--')
%     title([num2str(RL(b)) ' yr Flow Record'])
%     axis square
% end
% 
% subplot(2,2,1)
% ylabel('\bf Percent of Sediment Transported Above \itQ_{EFF}', 'FontSize',12)
% subplot(2,2,4)
% xlabel('\bf Daily Flow Coefficient of Variation','FontSize',12)
% legend('\beta = 0.5', '\beta = 1', '\beta = 2', '\beta = 3', ...
%                    '\beta = 4', '\beta = 5', 'Analytical Relation',...
%                    'Location','EastOutside')
% 
% %% Plot Qs Plus vs CVx for one FRL and all Beta
% figure
% %RL = [1000000 100000 1000 100 10];
% markers = cellstr(char('s','^','d','o','v','<','>','.','*','+','x'));
% RL = [10 100 1000 100000];
% colors = flipud(gray(length(cvX_plot)+1)); %summer % gray
% for b = 1:length(RL)
%     i = 1;
%     subplot(2,2,b)
%     for x = [1 3 6 16] % indices of cvX to plot
%         eval(['Qs_plot3(1,:) = Median_Qs_plus_',num2str(RL(b)),'(1,:,x);'])
%         plot(beta,Qs_plot3,[markers{i}],'Color',colors(i+1,:),...
%             'Linestyle','-','MarkerFaceColor',colors(i+1,:),...
%             'MarkerSize',5,'LineWidth',1.5)
%         hold on
%         i = i + 1;
%     end
%     title([num2str(RL(b)) ' yr Flow Record'])
%     axis square
%     axis([0 5 0 1])
% end
% 
% subplot(2,2,1)
% ylabel('\bf Percent of Sediment Transported Above \itQ_{EFF}', 'FontSize',12)
% subplot(2,2,4)
% xlabel('\bf Sediment Rating Curve Exponent Value, \beta', 'FontSize',12)
% legend('CV = 1.0','CV = 1.4','CV = 2.0','CV = 4.0',...
%        'Location','EastOutside')

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

% %% Effectiveness Histograms with median Mid_Trials (discharge) bin values
% NEEDS UPDATING: NOT SELECTING APPROPRIATE CVX VALUES
%
% % Chose Flow Record Length
% FRL = 100000;
% figure
% p = 1;
% for j = 1:length(cvX_plot)
%     for i = 1:length(beta_plot)
%         subplot(length(cvX_plot),length(beta_plot),p)
%         eval(['hbar = bar(Med_mid_trials_' num2str(FRL) '(j,:),Eff_med_' ...
%             num2str(FRL) '(j,:,i));'])
%         set(hbar,'FaceColor', 'k', 'EdgeColor', 'k');
%         set(gca, ...
%               'Box','on',...
%               'Layer','top',...
%               'TickDir','in',...
%               'TickLength',[.02 .02],...
%               'LineWidth',1)
%         p = p + 1;
%     end
% end
% for i = 1:length(beta_plot)
%     subplot(length(cvX_plot),length(beta_plot),i)
%     title(['\bf \beta = ' num2str(beta_plot(i))])
% end
% 
% p = 1;
% for i = 1:length(cvX_plot)
%     subplot(length(cvX_plot),length(beta_plot),p)
%     ylabel(['\bf\it CV_x = \rm\bf' num2str(cvX_plot(i))])
%     p = p+6;
% % subplot(length(cvX_plot),length(beta_plot),3)
% % title(['\bf Median Effectiveness Curves for ' num2str(FRL) ' yr Flow Record'])
% end
% hbar(Med_mid_trials_10(1,:),Eff_med_10(1,:))
% text(2,3, 'Sediment Discharge [mass/day]')
% h = text(2,4, 'MyStr3');
% set(h, 'rotation', 60)

RL = [10 100 100000];
trial = round(100*rand(1));
figure
p = 1;
for j = 1:length(RL)-1
    for i = 1:length(beta_plot)  
        subplot(length(RL),length(beta_plot),p)
        eval(['freq = Freq_trials_' num2str(RL(j)) '(trial,:,6);'])
        eval(['width = Mid_trials_' num2str(RL(j)) '(trial,2,6) -' ...
            'Mid_trials_' num2str(RL(j)) '(trial,1,6);'])
        density = freq./width;
        eval(['hbar = bar(Mid_trials_' num2str(RL(j)) '(trial,:,6),' ...
              'density);'])
        eval(['Q = linspace(0, max(Mid_trials_' num2str(RL(j)) ...
            '(trial,:,6)), 200);'])
        hold on
        plot(Q, lognpdf( Q, meanY1(cvX_plot(j)), stdY1(cvX_plot(j)) ),...
            'k-', 'LineWidth', 1.5);
        x_max = max(density);
        set(hbar,'FaceColor', 'b', 'EdgeColor', 'b')
        set(gca, ...
              'Box','off',...
              'Layer','top',...
              'TickDir','out',...
              'TickLength',[.02 .02],...
              'LineWidth',1,...
              'FontSize',8)
        p = p + 1;
    end
end
% 1000 yr record length has different number of CVx values (23 instead of
% 46)
for j = length(RL)
    for i = 1:length(beta_plot)  
        subplot(length(RL),length(beta_plot),p)
        eval(['eff = Freq_trials_' num2str(RL(j)) '(trial,:,3);'])
        eval(['hbar = bar(Mid_trials_' num2str(RL(j)) '(trial,:,3),' ...
              'eff);'])
        eval(['Q = linspace(0, max(Mid_trials_' num2str(RL(j)) ...
            '(trial,:,3)), 200);'])
        hold on
        plot(Q, lognpdf( Q, meanY2(3), stdY2(3) ),...
            'k-', 'LineWidth', 1.5);
        x_max = max(eff);
        set(hbar,'FaceColor', 'b', 'EdgeColor', 'b')
        set(gca, ...
              'Box','off',...
              'Layer','top',...
              'TickDir','out',...
              'TickLength',[.02 .02],...
              'LineWidth',1,...
              'FontSize',8)
        p = p + 1;
    end
end

y = Freq_trials_10(1,:,3);
x = Mid_trials_10(1,:,3);
width = x(2) - x(1);
Y = cumsum(y);
d_y = gradient(Y);
dy = y/width;
Q = linspace(0, max(Mid_trials_10(1,:,3)), 200);
ln_dy = lognpdf(Q,meanY2(3), stdY2(3));

FRL = 10;
ProbDaily = rand(1,FRL*365);
flows = logninv(ProbDaily,meanY2(3),stdY2(3));

[freq mid] = hist(flows, 50);

figure
plot(x,dy)
hold on
plot(Q,ln_dy, 'k--')

trapz(x,dy)
figure
plot(x,Y)

