function Output = UBCRM_base(S, Q, D50, D84, Mu, H, Tstar, RL, TL)
% UBCRM_core is a function to isolate the maximum efficency channel using the
% UBCRMH model where bank strength is represented by modified friction
% angle, varying from the basic friction angle, Phi to just under 90
% degrees
%
%%%% INPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q       Formative discharge (cumecs)
%
%   S       Energy gradient (m/m)
%
%   D50     Characteristic grain size for sediment transport(mm)
%
%   D84     grain size for bed stability analysis and flow resistance (mm)
%
%   Mu      bank strength expressed relative to the bed entrainment
%           threshold (dimensionless) NOTE: if H > 0, then Mu should be set
%           to 1.0, so that bank strength is attributed to the cohesive
%           term, H
%
%   H       bank strength as effective rooting depth (m) (optional) NOTE
%           that H should be set to 0 if Mu > 1, so that bank strength is 
%           attributed to a modified friction angle determined by Mu
%
%   Tstar   critical Shield's number for D84 (optional)
%
%   RL      flow resistance law selected (optional)
%               1 = Ferguson (2007)
%               2 = Keulegan
%               n (you can also simply enter a value of Manning's n)
%
%   TL      transport law  to be selected (optional)
%               1 = Parker (1990) 
%               2 = Eaton and Church (2011)    
%
%%%% OUTPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1 W           solution water surface width [m]
%
% 2 d           mean hydraulic depth [m]
%
% 3 U           mean flow velocity [m]
%
% 4 QB          bed material transport rate [me/s]
%
% 5 Qcrit       critical flow for entrainment of surface D50 (assuming
%               Shields number of 0.06)
%
% 6 bed_str     shear stress acting on the bed [Pa]
%
% 7 bank_str    shear stres action on the banks [Pa]
%
% 8 thetaSI     side slope angle [degrees]
%
% 9 Y           height of the trapazoid portion of the channel cross section
%
%
% Copyright 2014 Brett Eaton
    % <brett.eaton@ubc.ca>
    % Department of Geography, University of British Columbia
    %
    % Modified: July 2014
    %
    % This file is part of UBCRM.
    %
    % UBCRM is free software: you can redistribute it and/or modify
    % it under the terms of the GNU General Public License as published by
    % the Free Software Foundation, either version 3 of the License, or
    % (at your option) any later version.
    %
    % UBCRM is distributed in the hope that it will be useful,
    % but WITHOUT ANY WARRANTY; without even the implied warranty of
    % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    % GNU General Public License for more details.
    %
    % You should have received a copy of the GNU General Public License
    % along with Foobar.  If not, see <http://www.gnu.org/licenses/>.


    
    Phi = 35;       %set Phi to default of 35 if no value selected

if nargin < 6
    H = 0;          %set H to 0 if no cohsive term is entered
end
    
if nargin < 7       %set Tstar to default of 0.02 if no value selected (following Andrews, 1983)
    Tstar = 0.02;
end

ModPhi = atand(Mu*tand(Phi));   %translate Mu to a friction angle

if nargin < 8       %set RL to default of 1 if no value selected
    RL = 1;
end

if nargin < 9       %set TL to default of 1 if no value selected
    TL = 1;
end

%%%%%%%%%% STEP 1: CHOOSE AN INITIAL VALUE AND PERFORM INITIAL CALCS %%%%%
    Ptest = 3.0*Q^0.5;  % use hydraulic geometry equation to guess the first bed width

% calculate the local gradient in sediment transport at Ptest
    Y_minus = FindTransp(Ptest*0.99999, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
    Y_plus = FindTransp(Ptest*1.00001, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
    dQ1(1) = Y_plus - Y_minus;

% now move in direction of gradient towards solution
    if dQ1 > 0
        Ptest = Ptest +0.25*Ptest;
    else
        Ptest = Ptest - 0.25*Ptest;
    end

% perform a second calculation
    Y_minus = FindTransp(Ptest*0.99999, S, Q, D50, D84, H,ModPhi,  Phi, Tstar, RL, TL);
    Y_plus = FindTransp(Ptest*1.00001, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
    dQ1(2) = Y_plus - Y_minus;

%%%%%%%%%% STEP 2: SET UP A WHILE LOOP TO QUICKLY FIND THE MAXIMUM %%%%%%%
    i=2; 
    while dQ1(i-1)/dQ1(i) > 0    % continue loop until we get a reversal of the signs
        if dQ1(i) > 0
            Ptest = Ptest +0.25*Ptest;
        else
            Ptest = Ptest - 0.25*Ptest;
        end
        i = i + 1;
        Y_minus = FindTransp(Ptest*0.99999, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
        Y_plus = FindTransp(Ptest*1.00001, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
        dQ1(i) = Y_plus - Y_minus;
    end

%%%%%%%%%%%% STEP 3: NOW THAT WE HAVE IDENTIFIED THE BOUNDS WITHIN WHICH
%%%%%%%%%%%% THE SOLUTION LIES, WE SET THE BOUNDS AND FIND THE SOLUTION

% set the gates where P_upper is above the solution and P_lower is below it
    if dQ1(i) < 0
        P_upper = Ptest;      
        P_lower = Ptest - 0.25*Ptest;%move back one step to other side of max
    else
        P_lower = Ptest;
        P_upper = Ptest + 0.25*Ptest; %move forward one step to other side 
    end

% seek the maximum efficiency solution using Newtonian method
    Tol = 0.001;    %set the convergence toleranace, as a proportion of the solution bottom width
    while abs((P_upper - P_lower)/Ptest) > Tol
        %calculate the transport rate for a stable channel on either side of
        %the test bed width
        Y_minus = FindTransp(Ptest*0.99999, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
        Y_plus = FindTransp(Ptest*1.00001, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
        dQ2 = Y_plus - Y_minus;

        %MOVE THE GATES BASED ON GRADIENT dQ/dP
        if dQ2 < 0
            P_upper = Ptest;
        else
            P_lower = Ptest;
        end
        %recalculate the test bottom width
        Ptest = 0.5*(P_upper + P_lower);
    end
    
%%%%%%%%%%% STEP 4: CALCULATE THE FINAL SOLUTION %%%%%%%%%%%%%%%%%%%%%%%%
    Pbed = Ptest;           %record the solution bed width
    
    [QB, thetaSI] = ...  %extract the solution transport rate and geometry
        FindTransp(Pbed, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
    
    [~, Y, bed_str, bank_str] = ...   %extract trap depth and stresses
        FindSS(thetaSI, Pbed, S, Q, D84, H, ModPhi, Phi, Tstar, RL);
    
    [~, W, d, U, ~, ~] = ... %extract channel dimensions
        FindQ(Y, H, Pbed, thetaSI, S, D84, RL);
    
    Ycrit = 1.65 .* 0.06 .* (D50/1000) ./ S; %estimate the flow at which the bed becomes entrained
    if Ycrit > Y
        Hcrit = Ycrit - Y;
        Ycrit = Y;
    else
        Hcrit = 0;
    end
    [Qcrit] = ... %extract channel dimensions
        FindQ(Ycrit, Hcrit, Pbed, thetaSI, S, D84, RL);
    
     %generate output including shear stress calculations
    Output = [W, d, U, QB, Qcrit, bed_str, bank_str, thetaSI, Y];     
end

%%%%%%%%%%%%%%% END OF MAIN FUNCTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% SUB FUNCTIONS LISTED BELOW %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%% SUB FUNCTION FINDTRANSP, CALLED BY MAIN %%%%%%%%%%%%%%%%%
function [QB, thetaSI] = FindTransp(Pbed, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL)
%FindTransp.m is a function that takes a given channel width, slope,
%discharge, bed textures, and bank strength parameters and then solves for
%the bank angle for a critically stable bank, then calculates the sediment
%transport rate associated with that geometry. The function calls FindSS.m
%to calculate the stable channel geometry, which in turn calls FindQ.m to
%ensure that all geometries considered are consistent with continuity
%
%%%% INPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Pbed    the bottom width of the trapezoidal channel [m]
%   S       energy gradient [m/m]
%   Q       discharge carried by the channel [m3/s]
%   D50     grain size used for sediment transport calculation [mm]
%   D84     grain size used for flow resistance and stability analysis [mm]
%   H       bank strength parameter, expressed as effective rooting depth [m] 
%           (ranges from phi to <90 degrees)
%   Phi     friction angle for bank and bed sediment
%   Tstar   Critical Shields number for D84
%   RL      choice of flow resistance laws (see FindQ.m for details)
%%%% OUTPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   QB          transport rate for specified conditions [m3/s]
%   thetaSI     angle of the channel banks at critical stability [degrees]
%   YQ          the trapezoid depth associated with specified Q [m]

%determine stable bank angle for the given bed width associated with continuity
Trange = [0.001 ModPhi - 0.001];
%Trange = 0.75*Phi;  %test alternate solving routine
options = optimset('TolX', 10e-5);
thetaSI = fzero(@(x) FindSS(x, Pbed, S, Q, D84, H, ModPhi, Phi, Tstar, RL) - 1.0,...
    Trange,options);

%extract the solution shear stress values and critical shear stresses
[~, Y, bed_str] = FindSS(thetaSI, Pbed, S, Q, D84, H, ModPhi, Phi, Tstar, RL);

%calculate the sediment transport rate using a sediment transport law
%specified below
[~,~,~,U] = FindQ(Y, H, Pbed, thetaSI, S, D84, RL);

QB = TranspRate(D50, bed_str, U, S, Pbed, TL);

end

%%%%%%%%%%%%%%% SUB FUNCTION FINDSS, CALLED BY MAIN, FINDTRANSP %%%%%%%%%%%
function [SI, YQ, bed_str, bank_str, bed_crit, bank_crit] = FindSS(theta, Pbed, S, Q, D84, H, ModPhi, Phi, Tstar, RL)
%FindSS is a function that calculates the ratio of the critical bank
%strength and the shear stress acting onthe bank for a trapezoidal channel
%with known bed width, slope and discharge
%
%%%% INPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   theta   side slope angle [degrees]
%   Pbed    the bottom width of the trapezoidal channel [m]
%   S       energy gradient [m/m]
%   Q       discharge carried by the channel [m3/s]
%   D84     grain size used for flow resistance and stability analysis [mm]
%   H       bank strength parameter, expressed as effective rooting depth
%   [m]
%           (ranges from phi to <90 degrees)
%   Phi     friction angle for bank and bed sediment
%   Tstar   Critical Shields number for D84
%   RL      choice of flow resistance laws (see FindQ.m for details)
%%%% OUTPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   SI          bank stability index (critical/applied stress
%   YQ          the trapezoid depth associated with S1 [m]
%   bed_str     shear stress acting on the bank [Pa]
%   bank_str    shear stress acting on the bed [Pa]
%   bed_crit    critical shear stress for bed entrainment [Pa]
%


%%%%%%%%%%%%%%% determine depth associated with continuity %%%%%%%%%%%%%%%%
%     Yrange = [(0.0001*Q^0.5 - H) (5.0*Q^0.5)]; %allow Y to approach -H,
%     options = optimset('TolX', 10e-5);
%     YQ = fzero(@(x) FindQ(x, H, Pbed, theta, S, D84, RL) - Q, Yrange, options);

    %this sets the range of depths within which the solution lies
    
    %try using the bisection method
%     Yu = 100;
%     Yl = - H;
%     Yt = (Yu + Yl) / 2;
% 
%     TestQt = (FindQ(Yt, H, Pbed, theta, S, D84, RL) - Q)/Q;
%     Tol = 10e-6;
% 
%     while abs(TestQt) > Tol
%         if TestQt > 0
%             Yu = Yt;
%         elseif TestQt < 0
%             Yl = Yt;
%         end
%         Yt = (Yu + Yl) / 2;
%         TestQt = (FindQ(Yt, H, Pbed, theta, S, D84, RL) - Q)/Q;
%     end
%     YQ = Yt;
    
    
    %try using Newton Raphson method
    Yt = 0.33*Q^(1/3);
    Tol = 10e-6;
    Qt = FindQ(Yt, H, Pbed, theta, S, D84, RL);
    TestQt = (Qt - Q) / Q;
    while abs(TestQt) > Tol
        GradQ = (FindQ(1.0001*Yt, H, Pbed, theta, S, D84, RL) - ...
                 FindQ(0.9999*Yt, H, Pbed, theta, S, D84, RL))...
                 /(1.0001*Yt - 0.9999*Yt);
        Intercept = Qt - GradQ*Yt; %estimate the linear equation at Yt,Qt
        
        Yt = (Q-Intercept)/GradQ;   %use the equation to recalculate Yt
        Qt = FindQ(Yt, H, Pbed, theta, S, D84, RL); %check the Q value there
        TestQt = (Qt - Q) / Q;
    end
    YQ = Yt;
    
%%%%%%%%%%%%%%% partitiion the shear stress %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %dimensions when YQ <= 0 are based on a rectangle
    if YQ <= 0
        Ht = H + YQ;
        Pbank =  2 * (Ht);
        W = Pbed;
    end
    
    %dimensions when YQ > 0 are a composite trapezoid and rectangle
    if YQ > 0   
        Ht = H + YQ;
        Pbank = 2*H + 2*YQ/sin(theta*pi/180);
        W = Pbed + 2*YQ/tan(theta*pi/180);
    end
    
    totstress = 9.8*1000*(Ht)*S;  
    SFbank = 10^(-1.4026*log10(Pbed/Pbank + 1.5) + 0.247);    
    %partioning equation
    bed_str = totstress*(1 - SFbank)*(W/(2*Pbed) + 0.5);        
    %stress acting on the bed
    bank_str = totstress*SFbank*(W + Pbed)*sin(theta*pi/180)/(4*(Ht));  
    %stress acting on the bank
    %this is a set of  empirical equations derrived by Knight and others
    
%%%%%%%%%%%%%%% calculate the critical shear stress values %%%%%%%%%%%%%%
    bed_crit = Tstar*9.8*1650*(D84/1000);
    bank_crit = 9.8*1650*(D84/1000)*...
        Tstar*(tan(ModPhi*pi/180) / tan(Phi*pi/180))*...
        (1 - (sin(theta*pi/180)^2 / sin(ModPhi*pi/180)^2))^0.5;

%%%%%%%%%%%%%%% calculate a bank stability index %%%%%%%%%%%%%%%%%%%%%%%%
%     if YQ <= 0  %in this case, root cohesion dominates the bank strength
%         
%         SI = 1; %and we therefore set SI = 1, regardless of theta
%         
%     else    %if YQ > 0, then proceed with calculation, per normal
        
        SI = bank_crit/ bank_str;
        
%     end
end

%%%%%%%%%%%%%%% SUB FUNCTION FINDQ, CALLED BY MAIN, FINDSS %%%%%%%%%%%%%%%
function [Q, W, d, U, P, R] = FindQ(Y, H, Pbed, Theta, S, D84, RL)
%FindQ is a function to solve for the discharge associated with a set of
%input geometry parameters, a bed surface grain size, and a choice of flow
%resistance law to use, all assuming a trapezoidal channel geometry
%
%%% INPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Y       flow depth [m]
%   P       bottom width [m]
%   Theta   angle of the channel banks [degrees]
%   S       energy gradient [m/m]
%   RL      choice of resistance law to use 
%           1 = Ferguson (2007)
%           2 = Keulegan
%
%%%% OUTPUT = Q [m3/s]  %%%%%%%%%%%%%%%%%%%%%%%%
    
    %apply default flow resistance law
    if nargin < 7
        RL = 1;
    end
    
    %convert to useful units
    Theta = Theta * pi / 180;   %convert to radians
    D84 = D84 / 1000;           %convert to m
    
    %define constants
    g = 9.81;
    a1 = 5.5;   %constant for RL = 1
    a2 = 2.2;   %constant for RL = 1
    
    %calculate geometric parameters
    
    %when Y is negative, then the cross section is rectangular
    if Y <= 0
        dW = 0;
        Area = (H+Y)*Pbed;      %cross sectional area of flow
        P = Pbed + 2*(H+Y);     %wetted perimeter
        R = Area / P;           %hydraulic radius
    
    %when Y is positive, add a trapezoid to the dimensions of the upper
    %rectangle of height H
    elseif Y > 0
        dW = Y / tan(Theta);        %width of the bank at the base
        dP = sqrt(Y^2 + dW^2);      %length of bank
        Area = Y*(Pbed + dW)...     %cross section of lower trapezod
            + H*(Pbed + 2*dW);      %cross sectional  of upper rectangle
        P = Pbed + 2*dP + 2*H;      %wetted perimeter
        R = Area / P;               %hydraulic radius
    end
        %choose the appropriate flow resistance law
    if RL == 1
        %Use Ferguson's continuously varying power law
        Res = a1*a2*(R/D84) / (a1^2 + a2^2*(R/D84)^(5/3))^(1/2);
    elseif RL ==2
        %Keulegan equation
        Res = (1/0.4)*log(12.2*R/(D84));
    else
        %here is a backdoor to use a value of n for RL
        Res = R^(1/6) / (RL * sqrt(g));
    end
    
    %calculate velocity, U and discharge, Q
    U = Res*sqrt(g*R*S);
    Q = Area * U;
    W = Pbed + 2*dW;
    d = Area / W;
   
end

%%%%%%%%%%%%%%% SUB FUNCTION TRANSPRATE, CALLED BY FINDTRANSP %%%%%%%%%%%%
    function Qb = TranspRate(D50, tau, U, S, W, TL)
    %D50    the surface D50, used as the characteristic grain size
    %tau    the bed shear stress
    %U      the mean flow velocity, used to estimate stream power
    %W      the width of the sediment transport zone
    
    if TL == 1; %calculate sediment transport using Parker 1990 equation
        %define constants used in the equation
        rho = 1000; 
        rho_sed = 2650; 
        g = 9.81;

        % calculate dimensionless shear stress
        Tstar = tau/(g*(rho_sed-rho)*(D50/1000));
        Tref=0.0386;
        phi=Tstar/Tref;

        %determine which equation to use
        if phi > 1.59
            G=5474*(1-0.853/phi)^4.5;
            %t=1
        elseif phi < 1
            G=phi^14.2;
            %t=2
        else 
            G=exp(14.2*(phi-1)-9.28*(phi-1).^2);
            %t=3
        end

        %Calculate the transport rate
        Qb=W*G*(0.0025*(tau/rho)^(3/2)/(g*(rho_sed/rho-1)));
    end
    
    
    if TL == 2; %calculate using Eaton and Church 2011 relation
        %define constants used in the equation
        rho = 1000; 
        rho_sed = 2650; 
        g = 9.81;
        R = (rho_sed - rho)/rho;
        Tcrit = 0.045;      %this is the critical shields number
        
        %calculate stream power
        omega = tau * U;
        om_star = omega / (rho*(g*R*(D50/1000))^(3/2));
        
        dcrit = Tcrit * R * (D50/1000) / S;
        
        %Use Ferguson's continuously varying power law
        a1 = 6.5;   %constant for RL = 1
        a2 = 2.5;   %constant for RL = 1
        Res = a1*a2*(dcrit/(D50/1000)) / (a1^2 + a2^2*(dcrit/(D50/1000))^(5/3))^(1/2);
        
        om_crit = Res * Tcrit^(3/2);
        
        %calculate the transport rate
        E_star = (0.92 - 0.25*sqrt(om_crit/om_star))^9;
        qb_star = E_star * om_star;
        Qb = W * qb_star * sqrt(R*g)*(D50/1000)^(3/2);
    end
        
    
    
    end