function feta = fetaJR(phi, e_R)
% This function calculates the normalized shear viscosity eta from the 
% Jenkins & Richman (1985) theory.
% The normalized viscosity is defined as
% 
% feta = eta / [(5*pi^(1/2)/96)* rho_s * d * T^(1/2)]
%
% syntax:
%   feta = fetaJR(phi, e_R)
%
% inputs: volume fraction phi and the coefficient of restitution e_R, 
%
% output: normalized viscosity feta, its size is [length(phi), length(e_R)]
%
Nphi = length(phi);
Ne = length(e_R);
feta = zeros(Nphi, Ne);

chi = chiCS(phi);

for i=1:Ne
    ei = e_R(i);
    E1 = 4 / ((1+ei)*(3-ei));
    E2 = 8*(1+3*ei)/(5*(3-ei));
    E3 = 32*(1+ei)/(25*(3-ei)) * ((3*ei-1)+(12/pi)*(3-ei));
    
    feta(:,i) = E1./chi + E2*phi + E3*(phi.^2).*chi;
end
