sticky hard sphere fails for q<epsilon

[pkienzle]

Looking at the sticky hard sphere code it should converge in the limit of q → 0

Here is the calculation:

S(q) = 1 / (A^2 + B^2) => A^-2

for

A = 1 + 12 eta [A1 + A2 - A3] => 1 + 12 eta [alpha/3 + beta/2 - lambda/12]
B = 12 eta [B1 + B2 - B3] => 0

A1 = alpha [(sin(qr) - qr*cos)/(qr)^3] = alpha/3 3j1(qr)/qr => alpha/3
A2 = beta [(1 - cos(qr))/(qr)^2] => beta/2
A3 = lambda/12 [sin(qr)/qr] => lambda/12

B1 = alpha [ 1/(2qr) - sin(qr)/(qr)^2 + (1-cos(qr))/(qr)^3] => 0
B2 = beta [ 1 / qr - sin(qr)/(qr)^2 ] => 0
B3 = lambda/12 [ (1 - cos(qr))/qr ] => 0

where

alpha = (1 + 2 eta - mu)/(1 - eta)^2
beta = (mu - 3 eta)/(2 (1-eta))
mu = lambda eta (1-eta)
lambda = (-b - sqrt(b^2 - 4 a c))/2a
a = eta/6
b = tau + eta/(1-eta)
c = (1 + eta/2)/(1 - eta)^2
eta = Vf/(1-epsilon)^3

and

r = 2 Re / (1 - epsilon)

using the input parameters:

Re is radius_effective,
Vf is volume fraction, 
epsilon is perturbation and 
tau is stickiness.

We already have an accurate calculator for 3 j1(x)/x and sin(x)/x

Need to check accuracy of the following for low q, possibly replacing them with the Taylor expansion

    [1-cos(x)]/x^2 = 1/2 - x^2/24 + x^4/720 - x^6/40320 + ...
    1/(2x) - sin(x)/x^2 + (1 - cos(x))/x^3 = x/8 - x^3/144 + x^5/5760 - x^7/403200 + ...
    1/x - sin(x)/x^2 = x/6 - x^3/120 + x^5/5040 - x^7/362880 + ...