function [SiNi] = SI_NI(r)

% Calculates the ratio of collisional shear and normal stresses
% on a flat, frictional wall using an empirical expression fit 
% from the numerical simmulations described in Louge, M.Y. 
% "Computer simulations of rapid granular flows of spheres 
% interacting with a flat, frictional boundary", 
% Physics of Fluids 6(7), pp. 2253-2269 (1994).
%
% Michel Louge 4/11/06
%
%	input	-   r           =	relative velocity at the contact point
%                                   divided by sqrt[3*granular_temperature_there]
%
%	output	-	SiNi        =   ratio of collisional shear and normal stresses
%
%   global variables
%   
%           -   muw         =   collisional friction with flat wall
%           -   ew          =   collisional normal restitution with flat wall
%           -   b0w         =   collisional tangential restitution with flat wall
%
%               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

% 	other internal variables
%
%           -   a1          =   constant in Louge's fit of SiNi
%           -   a2          =   ditto
%           -   a3          =   ditto
%           -   a4          =   ditto
%           -   a5          =   ditto
%           -   a           =   empirical coefficient
%           -   aprime      =   empirical exponent
%

if muw < 0 | muw == 0
    disp(['impossible collisional friction = ' num2str(muw)])
    pause
end

if ew < 0 | ew > 1
    disp(['impossible normal restitution = ' num2str(ew)])
    pause
end

if b0w < 0 | b0w > 1
    disp(['impossible tangential restitution = ' num2str(b0w)])
    pause
end

a1      =   0.3537  ;
a2      =   1.042   ;
a3      =   4.453   ;
a4      =   2.068   ;
a5      =   0.8468  ;

aprime  =   a4*(1-exp(-a5*muw))                                 ;
a       =   a1*exp(-a2*(1-ew))*(1-exp(-a3*muw))*(1+b0w)^aprime  ;

SiNi    =   muw*(1-exp(-a*r/muw))                               ;

%
% Louge [PF 6 2253-69 (94)] reports the fit above for SiNi vs. r
% for arbitrary values of collisional parameters.
%
% However, note that he implemented a different, closer fit to his 
% simulations carried out with the specific collisional properties 
% of his spent glass beads in Louge & Keast PF 15, 1213-33 (01).
% This latter fit is
% SiNi	=	muw*(1-1/(1+r^0.95785*(1.28816*r +0.898728)))
% 