function [error]  =   find_r(r)

% Helps find the dim(-) contact point velocity
%
% Michel Louge 4/14/06
%
%   input   -   r       =   contact point velocity/sqrt(3 T)
%
%   output  -   error	=   error to annul using fzero('find_r',rguess,options)
%
%   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
%           -   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
%

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)                                         ;
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)          ;
        
etacalc =   muE/sni*(sqrh^2/nu0gew/Fr^2 - 1)                    ;

error   =   real(etacalc - eta)                                 ;   % use real part to make more robust
