function [qiNiTi] = qI_NI(r)

% Calculates the collisional flux of fluctuation energy
% through a flat, frictional surface, counted positive
% toward the flow and
% dim(-) with the product of collisional normal stress
% and sqrt[3*granular_temperature_there] 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	-	qiNiTi      =   collisional flux of fluctuation energy divided by
%                               the product of collisional normal stress
%                               and sqrt[3*granular_temperature_there]
%
%   global variables
%   
%           -   muw         =   collisional friction
%           -   ew          =   collisional normal restitution
%           -   b0w         =   collisional tangential restitution
%
%               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
%
%           -   b1          =   constant in Louge's fit of qiNiTi at low r
%           -   b2          =   ditto
%           -   b3          =   ditto
%           -   b4          =   ditto
%           -   b5          =   ditto
%           -   b6          =   ditto
%           -   b7          =   ditto
%           -   b8          =   ditto
%           -   Bb1         =   empirical function in Louge's fit at low r
%           -   Bb2         =   ditto
%           -   Bb3         =   ditto
%           -   Bb4         =   ditto
%           -   Bb5         =   ditto
%           -   d1          =   constant in Louge's fit of qiNiTi at large r
%           -   d2          =   ditto
%           -   d3          =   ditto
%           -   qiNiTi0     =   empirical function for qiNiTi with r->0
%           -   qiNiTi      =   empirical function for qiNiTi with r->°
%

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

b1      =   0.02	;
b2      =   0.3     ;
b3      =   3.4     ;
b4      =   1.8     ;
b5      =   1.1     ;
b6      =   0.1     ;
b7      =   4.6     ;
b8      =   8.6     ;

Bb1     =   b1/(1+ew)                       ;
Bb2     =   b2*(1-ew)*(1-exp(-b3*muw))      ;
Bb4     =   b4*(1-ew)^2 + b5*(1-ew) + b6    ;
Bb5     =   b7*(1-ew)+b8                    ;
Bb3     =   Bb4*(1-exp(-Bb5*muw))           ;

qiNiTi0	=   Bb1*r^2 + Bb2*r - Bb3                           ;

d1      =   0.522   ;
d2      =   2.57    ;
d3      =   1.08    ;

qiNiTif	=   d1*(1+ew)*muw^2 + d2*(1-ew)*muw - d3*(1-ew)     ;

qiNiTi  =   min(qiNiTi0,qiNiTif)                            ;

%
% Louge [PF 6 2253-69 (94)] reports the fit above for qiNiTi 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
%
% qiNiTi    =   qiNiTi0 + (qiNiTif - qiNiTi0)*[1-exp(-0.24*r^1.38)],
%                   where qiNiTi0 Å -0.096 and qiNiTif Å -0.22
%                   indicate a negative flux for any r.
% 