function [r,mueff] = flat_chute(Fri,holdupi,nu0i,ewi,muwi,b0wi,eei,muui,b00i,muEi,rguess)

% Calculates the effective friction on a flat, frictional
% base in terms of the Froude number for a given value of the
% solid volume fraction in the basal layer. Model described in
% Louge, M.Y. and Keast S.C. 
% "On dense granular flows down flat frictional inclines",
% Physics of Fluids 15(5), pp. 1213-1233 (2001).
%
% Principal model assumptions:
%
%   spherical, monodisperse spheres
%   velocity parallel to base, flow parameters gradients perpendicular to base
%   static friction of flat surface = muE
%   collisional friction of flat surface = muI
%   flow has thin shear layer at base with dense passive overburden
%   granular temperature adjusts instantly to shear layer conditions
%   enduring and collisional stresses coexist and addup
%   spin of spheres colliding with the base given by simulations 
%       in terms of relative velocity at the contact point divided by
%       the square root of the granular temperature there
%   
%   Note:	For the calculation implemented here, the flow does NOT need to be
%           steady and fully-developed. A crucial input is the solid volume
%           fraction in the shear layer. Such volume fraction can be
%           calculated elsewhere by equating the production and
%           dissipation rates of angular momentum in the shear layer.
%
% Michel Louge 4/17/06
%
%	input	-	** Note our convention that global variables are input with a name followed by the letter i **
%
%           -   Fri         =	Froude number       = mfr/holdup^(3/2)/sqrt(cos(alpha))
%                                   where	mfr     = mass flow rate dim(-)
%                                                       with rho_s*sqrt(gD)*W*D
%                                           holdup  = integral of volume
%                                                       fraction along z divided by D
%                                           D       = sphere diameter
%                                           g       = acceleration of gravity
%                                           alpha   = angle of inclination
%                                           rho_s   = sphere material density
%                                           W       =   chute width
%                                       
%			-	holdupi     =	see above
%           -   nu0i        =   volume fraction in the basal layer
%           -   muwi        =   collisional friction with flat wall
%           -   ewi         =   collisional normal restitution with flat wall
%           -   b0wi        =   collisional tangential restitution with flat wall
%           -   muui        =   collisional binary friction
%           -   eei         =   collisional binary normal restitution
%           -   b00i        =   collisional binary tangential restitution
%           -   muEi        =   enduring friction of spheres with the base
%           -   rguess      =   guess for the contact point velocity/sqrt(3T)
%
%	output	-	mueff       =   effective friction of the flat, frictional surface on the granular flow
%           -   r           =   contact point velocity/sqrt(3T)
%
%   global variables
%   
%           -	Fr          =	Froude number
%           -   holdup      =   holdup
%           -   muw         =   collisional friction with flat wall
%           -   ew          =   collisional normal restitution with flat wall
%           -   b0w         =   collisional tangential restitution with flat wall
%           -   muu         =   collisional binary friction
%           -   ee          =   collisional binary normal restitution
%           -   b00         =   collisional binary tangential restitution
%           -   muE         =   enduring friction of spheres with the base
%           -   nuL         =   volume fraction at the erodible boundary
%                               (Louge & Keast assumed nuL = 0.55)
%           -   nu0         =   volume fraction in the basal layer
%
%               For the measurement technique of these collisional parameters, see
%               S. Foerster, M. Y. Louge, H. Chang & K. Allia
%               "Measuremenfs of the collision properties of
%               small spheres,Ó Phys. Fluids 6, 1108 (1994).
%               For the measurements of the specific spheres used by Louge
%               & Keast, see A. Lorenz, C. Tuozzolo, and M. Y. Louge, 
%               "Measurements of impact properties of small, nearly 
%               spherical particles," Experimental Mechanics 37, 292-298 (1997).
%           

global muw ew b0w muu ee b00 muE nuL nu0 Fr holdup

% 	other internal variables
%
%           -   mm  	=   1+9*pi/32           appears in the dense theory
%           -   jj      =   1+pi/12             ditto
%           -   xx      =   sqrt(3*pi)/2/mm     ditto
%           -   g12L    =   HAB sphere-flat wall pair distribution at contact
%           -   g12CS   =   C&S pair distribution at contact
%           -   nu0gew  =   parameter in wall eqn of state
%           -   gam     =   right-hand-side of eq.(46) in Louge & Keast
%           -   r       =   contact point velocity/sqrt(3 T)
%           -   sni     =   (SI/NI) ... calculated once only to save time
%           -   sni2    =   (SI/NI)^2 ... calculated to save time
%           -   b1      =   constant in eq. (37) of Louge & Keast
%           -   b00     =   flux-dependent part in eq. (41) of Louge et Keast
%           -   BB      =   
%           -   deltap  =  
%           -   eta     =   SE/SI =	enduring shear stress/collisional shear stress
%                               Note!! This definition ot eta is different
%                               than that of Louge, M.: "Model for dense
%                               granular flows down bumpy inclines", 
%                               Physical Review E 67, 061303 (2003).
%           -   b0p     =   eq. (41) in Louge & Keast
%           -   k2      =   defined in eq. (37) of Louge & Keast
%           -   thick	=   shear layer thickness divided by D e.g., eq (45) of Louge & Keast
%           -   sqrh    =   sqrt(h) in eq. (53,54) of Louge & Keast
%           -   etacalc =   calculated value of eta
%

% initialization of global variables

muw     =   muwi        ;   % wall collisions
ew      =   ewi         ;
b0w     =   b0wi        ;

muu     =   muui        ;   % binary collisions
ee      =   eei         ;
b00     =   b00i        ;

muE     =   muEi        ;   % static friction of the base

Fr      =   Fri         ;   % Froude
holdup  =   holdupi     ;   % holdup

nu0     =   nu0i        ;   % volume fraction in the basal layer

nuL     =   0.55        ;   % volume fraction at the erodible boundary between shear layer and overburden

% check validity of input parameters

if muw == 0 | muw < 0
    disp(['warning: muw ² 0'])
    return
end

if ew < 0 | ew > 1
    disp(['warning: invalid ew'])
    return
end

if b0w < 0 | b0w > 1
    disp(['warning: invalid b0w'])
    return
end

if muu == 0 | muu < 0
    disp(['warning: mu ² 0'])
    return
end

if ee < 0 | ee > 1
    disp(['warning: invalid e'])
    return
end

if b00 < 0 | b00 > 1
    disp(['warning: invalid b0'])
    return
end

if muE < 0
    disp(['warning: invalid muE'])
    return
end

if Fr < 0 | Fr == 0
    disp(['warning: invalid Fr'])
    return
end

if holdup < 0
    disp(['warning: invalid holdup'])
    return
end

if nu0 < 0 | nu0 > 0.64
    disp(['warning: invalid basal volume fraction'])
    return
end

% calculate parameters

Eeff    =   ee - pi/2*muu                   ;   % effective restitution

mm      =   1+9*pi/32                       ;   % various constants in the dense theory
jj      =   1+pi/12                         ;
xx      =   sqrt(3*pi)/2/mm                 ;
rr      =   5*pi/4/mm/jj                    ;
pj      =   5*sqrt(pi)/2/jj                 ;

g12L    =   (1+2*nu0)/(1-nu0)^2             ;   % wall pair-distribution function
g12CS   =   (2-nu0)/2/(1-nu0)^3             ;   % Carnahan & Starling pair distribution function in the shear layer

nu0gew  =   nu0*g12L*(1+ew)/2                   ;   % various parameters involved in the solution
gam     =   nu0gew*(1-nuL)^3/2/nuL^2/(2-nuL)    ;

b1      =   -sqrt(3*(1-Eeff)/mm)                ;

% find the value of r

options = optimset('TolFun',1e-9)       ;
r =   fzero('find_r',rguess,options)         ;

% calculate results

b00     =   xx*qI_NI(r)                                         ;
sni     =   SI_NI(r)                                            ;
sni2    =   sni^2                                               ;
BB      =   sni2*rr*(1-gam)+2*b00*xx*r*sni                      ;
deltap  =   BB^2+4*xx^2*r^2*sni2*(b1^2-b00^2-rr*sni2*(1-gam))	;
eta     =   (-BB+sqrt(deltap))/2/xx^2/r^2/sni2                  ;
b0p     =   b00+xx*r*sni*eta                                    ;
k2      =   (b0p^2-gam*b1^2)/(1-gam)                            ;
if k2 > 0
    thick   =   (1/sqrt(k2))*atanh(sqrt(b0p^2-gam*b1^2)*...
                                    sqrt(1-gam)/(b0p-gam*b1))   ;
else
    if k2 == 0
        thick   =   (1/b0p)+(1/b1)                              ;
    else
        thick   =   (1/sqrt(-k2))*atan(sqrt(-b0p^2+gam*b1^2)*...
            sqrt(1-gam)/(b0p-gam*b1))                           ;
    end
end

sqrh    =    pj/holdup*sni/k2*(-b0p*nu0/2+nu0*...
            ( sqrt( (b0p^2-k2)/(b1^2-k2) ) - 1)+...
            b1*(holdup-nu0*(thick+0.5))*...
            sqrt( (b0p^2-k2)/(b1^2-k2) )+...
            b0p*holdup ) + sqrt(3)*spin(r) + r*sqrt(3)          ;
        
% check the existence of a solution, see Louge & Keast section E p. 1226:

if k2 < 0
    if b0p > abs(b1) | gam < b0p/b1
            disp(['no solution with k^2 < 0'])
            pause
    end
else
    if b0p < 0 | b0p > abs(b1)
             disp(['no solution with k^2 > 0'])
             pause
    end
end

if thick > holdup/nu0 - 0.5
    disp(['warning: excessive basal layer thickness'])
    pause
end

if sqrh < 0
    disp(['warning: sqrth < 0'])
    pause
end


mueff   =   muE - nu0gew*(muE - sni)/sqrh^2*Fr^2                ;
