function [Fres] = SFD_flat_chute(tanalphai,holdupi,ewi,muwi,b0wi,eei,muui,b00i,muEi,guess)

% Calculates and plots the Froude number and basal volume fraction 
% for steady, fully-developed (SFD) flow on a flat, frictional
% base in terms of the tangent of inclination and properties of
% the spheres and wall. 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).  Closure is
% achieved by assuming a steady balance of spin production
% and dissipation, see eq. (68) of Louge & Keast. Other model
% assumptions, see flat_chute.m
%
% Also plot the curves mueff vs Froude for general flows (SFD or not)
% as well as the data of Louge & Keast for comparison.
%
% Michel Louge 4/18/06.
%
%	input	-	** Note our convention that global variables are input with a name followed by the letter i **
%
%           -   tanalphai   =   line array of tangents of angle of inclination
%			-	holdupi     =	integral of volume fraction along z divided by D
%           -   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
%           -   guess       =   line vector of initial guess [r nu0]
%
%	output	-   Fres        =	array of results [Fr nu0 r mueff]
%           -                   with    Fr	= 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
%           -                           nu0	=	volume fraction in the basal layer
%           -                           r	=	contact point velocity/sqrt(3T)
%           -                           mueff	=	effective friction of the flat, frictional surface 
%                                                   on the granular flow (should equal tanalpha)
%   Example: the command...
%
%       [Fres] = SFD_flat_chute([.3],5,.8,.141,.35,.972,.177,.25,.593,[2 .5])
%
%   ... calculates Fr, nu0, r and mueff for holdup=5 at tanalpha = 0.3 for the particle 
%       properties of Louge & Keast's experiments; mueff should equate tanalpha.
%
%   Example: the command...
%
%       [Fres] = SFD_flat_chute([.278:.001:.364],5,.8,.141,.35,.972,.177,.25,.593,[2 .5])
%
%   ... calculates Fr, nu0, r and mueff for holdup=5 for values of tanalpha
%       from 0.278 to 0.364 in an increment of 0.001.
%
%   global variables
%   
%           -   tanalpha    =   tangent of angle of inclination
%           -   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)
%
%               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 holdup tanalpha

% 	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
%           -   res     =   line vector of [r nu0] minimizing error in find_rnu.m
%           -   cst     =   constant of o(1) providing best fit of Louge & Keast data to 
%                           predictions of this prgm
%           -   num     =   number of tanalpha considered
%           -   i       =   index
%           -   Fr
%           -   nu0
%           -   r
%           -   mueff
%           -   numm    =   number of points for plotting mueff vs. Fr
%           -   Frmin	=	minimum Froude number of interest
%           -   Frmax	=   maximum Froude number of interest
%           -   rguess	=   initial guess for the contact point velocity/sqrt(3T) at Frmin
%                               (This may require trial and error - use flat_chute.m to that end).
%           -   Frr  	=   array of Froude numbers > 0
%           -   rr      =   value of r from solution
%           -   muff	=   corresponding array of effective friction
%           -   muef    =   placeholder value of mueff
%           -   ii      =   index
%           -   data    =   line array of [tanalpha Fr*] containing Louge & Keast data
%                                   where Fr* is the Froude number at SFD
%                                   for the tangent of inclination tanalpha
%           -   nu00    =   array of values of nu0 for plotting mueff vs Fr curves
%

num         =   size(tanalphai,2)   ;

% 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

holdup      =   holdupi     ;   % holdup

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 guess(1) < 0 | guess(1) == 0
    disp(['warning: invalid guess for r'])
    return
end

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

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

% loop over all tanalpha considered

for i=1:num
    
tanalpha    =   tanalphai(1,i)	;	% tangent of inclination

% find the value of r and nu0 at SFD

options = optimset('TolFun',1e-9)                   ;
res     =   fminsearch('find_rnu',guess,options)	;
guess   =   res                                     ;   % reset the guess

r       =   res(1)                                  ;
nu0     =   res(2)                                  ;

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

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

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 interior pair distribution function

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)                ;

b00     =   xx*qI_NI(r)                                         ;  % using the simulations fit in Louge PF 94
sni     =   SI_NI(r)                                            ;  % using the simulations fit in Louge PF 94
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)          ;

cst     =   0.7839	;   % constant of o(1) providing best fit of Louge & Keast data to predictions of this prgm
        
Fr      =   sqrt( pi/3/nu0 * ( (muE+muu)/nuL+muu/holdup )/...
                (1+ee)/muu )*cst * sqrh/sqrt(3)/(r+spin(r))     ;
            
mueff   =   muE - nu0gew*(muE - sni)/sqrh^2*Fr^2                ;

Fres(i,1)   =   Fr      ;
Fres(i,2)   =   nu0     ;
Fres(i,3)   =   r       ;
Fres(i,4)   =   mueff   ;

end

if num ~= 1

% calculates and plot the curves of Fr vs mueff for three values of nu0

for ii=1:3
nu00(ii)        =   min(Fres(:,2))+(max(Fres(:,2))-min(Fres(:,2)))/2*(ii-1) ;

numm            =   100                     ;   % describe the curve with 100 points
Frmin           =   0.1                     ;   % sets the minimum Froude to 0.1
Frmax           =   max(10,max(Fres(:,1)))  ;   % maximum Froude is at least 10
rguess          =   0.4                     ;   % guess for r at Frmin

for i=1:numm
    Frr(i,ii)	=   Frmin+(Frmax-Frmin)/(numm-1)*(i-1)                                              ;
    [rr,muef]	=   flat_chute(Frr(i,ii),holdupi,nu00(ii),ewi,muwi,b0wi,eei,muui,b00i,muEi,rguess)	;
    muff(i,ii)	=   muef                                                                            ;
    rguess      =   rr                                                                              ;
end

end

% experimental data of Louge & Keast for [tanalpha   Fr*] plotted along for comparison:

data(1,:) = [0.363970234	4.230994989]    ;
data(2,:) = [0.344327613	3.040989095]    ;
data(3,:) = [0.324919696	2.724538784]    ;
data(4,:) = [0.315298789	2.495716711]    ;
data(5,:) = [0.305730681	2.213729173]    ;
data(6,:) = [0.296213495	2.031964493]    ;
data(7,:) = [0.286745386	1.981008271]    ;
data(8,:) = [0.277324544	1.722075365]    ;

scrsz = get(0,'ScreenSize');
figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/3])

subplot(1,2,1)
plot(Fres(:,1),Fres(:,4),'r',Frr(:,1),muff(:,1),'b',...
    Frr(:,2),muff(:,2),'g',Frr(:,3),muff(:,3),'m',data(:,2),data(:,1),'o')
xlabel('Fr')
ylabel('effective friction')
axis([0 1.1*Frmax 0 muE])
title(['nu0 = ' num2str(nu00(1),2) ' (blue) & '  ...
    num2str(nu00(2),2) ' (green) & ' ...
    num2str(nu00(3),2) ' (magenta)'])

subplot(1,2,2)
plot(Fres(:,4),Fres(:,2),'r')
xlabel('tanalpha')
ylabel('nu0')
title(['base layer volume fraction vs tanalpha at SFD'])

end

