Index: sasmodels/models/_py_parallelepiped.py
===================================================================
--- sasmodels/models/_py_parallelepiped.py	(revision 597878a0f8b1b1b38cfce28323dfa620bdbaea50)
+++ sasmodels/models/_py_parallelepiped.py	(revision 597878a0f8b1b1b38cfce28323dfa620bdbaea50)
@@ -0,0 +1,316 @@
+# parallelepiped model
+# Note: model title and parameter table are inserted automatically
+r"""
+The form factor is normalized by the particle volume.
+For information about polarised and magnetic scattering, see
+the :ref:`magnetism` documentation.
+
+Definition
+----------
+
+ This model calculates the scattering from a rectangular parallelepiped
+ (\:numref:`parallelepiped-image`\).
+ If you need to apply polydispersity, see also :ref:`rectangular-prism`.
+
+.. _parallelepiped-image:
+
+
+.. figure:: img/parallelepiped_geometry.jpg
+
+   Parallelepiped with the corresponding definition of sides.
+
+.. note::
+
+   The edge of the solid used to have to satisfy the condition that $A < B < C$.
+   After some improvements to the effective radius calculation, used with
+   an S(Q), it is beleived that this is no longer the case.
+
+The 1D scattering intensity $I(q)$ is calculated as:
+
+.. Comment by Miguel Gonzalez:
+   I am modifying the original text because I find the notation a little bit
+   confusing. I think that in most textbooks/papers, the notation P(Q) is
+   used for the form factor (adim, P(Q=0)=1), although F(q) seems also to
+   be used. But here (as for many other models), P(q) is used to represent
+   the scattering intensity (in cm-1 normally). It would be good to agree on
+   a common notation.
+
+.. math::
+
+    I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2
+           \left< P(q, \alpha) \right> + \text{background}
+
+where the volume $V = A B C$, the contrast is defined as
+$\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$,
+$P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented
+at an angle $\alpha$ (angle between the long axis C and $\vec q$),
+and the averaging $\left<\ldots\right>$ is applied over all orientations.
+
+Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the
+form factor is given by (Mittelbach and Porod, 1961)
+
+.. math::
+
+    P(q, \alpha) = \int_0^1 \phi_Q\left(\mu \sqrt{1-\sigma^2},a\right)
+        \left[S(\mu c \sigma/2)\right]^2 d\sigma
+
+with
+
+.. math::
+
+    \phi_Q(\mu,a) &= \int_0^1
+        \left\{S\left[\frac{\mu}{2}\cos\left(\frac{\pi}{2}u\right)\right]
+               S\left[\frac{\mu a}{2}\sin\left(\frac{\pi}{2}u\right)\right]
+               \right\}^2 du
+
+    S(x) &= \frac{\sin x}{x}
+
+    \mu &= qB
+
+The scattering intensity per unit volume is returned in units of |cm^-1|.
+
+NB: The 2nd virial coefficient of the parallelepiped is calculated based on
+the averaged effective radius, after appropriately sorting the three
+dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and
+length $(= C)$ values, and used as the effective radius for
+$S(q)$ when $P(q) \cdot S(q)$ is applied.
+
+To provide easy access to the orientation of the parallelepiped, we define
+three angles $\theta$, $\phi$ and $\Psi$. The definition of $\theta$ and
+$\phi$ is the same as for the cylinder model (see also figures below).
+
+.. Comment by Miguel Gonzalez:
+   The following text has been commented because I think there are two
+   mistakes. Psi is the rotational angle around C (but I cannot understand
+   what it means against the q plane) and psi=0 corresponds to a||x and b||y.
+
+   The angle $\Psi$ is the rotational angle around the $C$ axis against
+   the $q$ plane. For example, $\Psi = 0$ when the $B$ axis is parallel
+   to the $x$-axis of the detector.
+
+The angle $\Psi$ is the rotational angle around the $C$ axis.
+For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the $B$ axis
+oriented parallel to the y-axis of the detector with $A$ along the z-axis.
+For other $\theta$, $\phi$ values, the parallelepiped has to be first rotated
+$\theta$ degrees around $z$ and $\phi$ degrees around $y$,
+before doing a final rotation of $\Psi$ degrees around the resulting $C$ to
+obtain the final orientation of the parallelepiped.
+For example, for $\theta = 0$ and $\phi = 90$, we have that $\Psi = 0$
+corresponds to $A$ along $x$ and $B$ along $y$,
+while for $\theta = 90$ and $\phi = 0$, $\Psi = 0$ corresponds to
+$A$ along $z$ and $B$ along $x$.
+
+.. _parallelepiped-orientation:
+
+.. figure:: img/parallelepiped_angle_definition.png
+
+    Definition of the angles for oriented parallelepiped, shown with $A<B<C$.
+
+.. figure:: img/parallelepiped_angle_projection.png
+
+    Examples of the angles for an oriented parallelepiped against the
+    detector plane.
+
+On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will
+appear. These are actually rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped, perpendicular to the $a$ x $c$ and $b$ x $c$ faces. 
+(When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.) The third orientation distribution, in $\psi$, is 
+about the $c$ axis of the particle, perpendicular to the $a$ x $b$ face. Some experimentation may be required to 
+understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation 
+distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.) 
+
+    
+For a given orientation of the parallelepiped, the 2D form factor is
+calculated as
+
+.. math::
+
+    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2
+                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2
+                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2
+
+with
+
+.. math::
+
+    \cos\alpha &= \hat A \cdot \hat q,
+
+    \cos\beta  &= \hat B \cdot \hat q,
+
+    \cos\gamma &= \hat C \cdot \hat q
+
+and the scattering intensity as:
+
+.. math::
+
+    I(q_x, q_y) = \frac{\text{scale}}{V} V^2 \Delta\rho^2 P(q_x, q_y)
+            + \text{background}
+
+.. Comment by Miguel Gonzalez:
+   This reflects the logic of the code, as in parallelepiped.c the call
+   to _pkernel returns $P(q_x, q_y)$ and then this is multiplied by
+   $V^2 * (\Delta \rho)^2$. And finally outside parallelepiped.c it will be
+   multiplied by scale, normalized by $V$ and the background added. But
+   mathematically it makes more sense to write
+   $I(q_x, q_y) = \text{scale} V \Delta\rho^2 P(q_x, q_y) + \text{background}$,
+   with scale being the volume fraction.
+
+
+Validation
+----------
+
+Validation of the code was done by comparing the output of the 1D calculation
+to the angular average of the output of a 2D calculation over all possible
+angles.
+
+
+References
+----------
+
+P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
+
+R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
+
+Authorship and Verification
+----------------------------
+
+* **Author:** This model is based on form factor calculations implemented
+    in a c-library provided by the NIST Center for Neutron Research (Kline, 2006).
+* **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017
+* **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017
+
+"""
+
+import numpy as np
+from numpy import pi, inf, sqrt, sin, cos, radians
+from sasmodels.special import sas_sinx_x, Gauss76Wt, Gauss76Z, square, ORIENT_ASYMMETRIC
+
+name = "parallelepiped"
+title = "Rectangular parallelepiped with uniform scattering length density."
+description = """
+    I(q)= scale*V*(sld - sld_solvent)^2*P(q,alpha)+background
+        P(q,alpha) = integral from 0 to 1 of ...
+           phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma
+        with
+            phi(mu,a) = integral from 0 to 1 of ..
+            (S((/mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du
+            S(x) = sin(x)/x
+            mu = q*B
+        V: Volume of the rectangular parallelepiped
+        alpha: angle between the long axis of the
+            parallelepiped and the q-vector for 1D
+"""
+category = "shape:parallelepiped"
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
+               "Parallelepiped scattering length density"],
+              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
+               "Solvent scattering length density"],
+              ["length_a", "Ang", 35, [0, inf], "volume",
+               "Shorter side of the parallelepiped"],
+              ["length_b", "Ang", 75, [0, inf], "volume",
+               "Second side of the parallelepiped"],
+              ["length_c", "Ang", 400, [0, inf], "volume",
+               "Larger side of the parallelepiped"],
+              ["theta", "degrees", 60, [-360, 360], "orientation",
+               "c axis to beam angle"],
+              ["phi", "degrees", 60, [-360, 360], "orientation",
+               "rotation about beam"],
+              ["psi", "degrees", 60, [-360, 360], "orientation",
+               "rotation about c axis"],
+             ]
+
+def form_volume(length_a, length_b, length_c):
+    return length_a * length_b * length_c
+             
+def Iq(q,
+        sld, solvent_sld,
+        length_a, length_b, length_c):
+
+    mu = 0.5 * q * length_b
+    
+    # Scale sides by B
+    a_scaled = length_a / length_b
+    c_scaled = length_c / length_b
+        
+    # outer integral (with gauss points), integration limits = 0, 1
+    outer_total = 0 #initialize integral
+
+    for i in range(76):
+        sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
+        mu_proj = mu * sqrt(1.0-sigma*sigma)
+
+        # inner integral (with gauss points), integration limits = 0, 1
+        # corresponding to angles from 0 to pi/2.
+        inner_total = 0.0
+        for j in range(76):
+            uu = 0.5 * ( Gauss76Z[j] + 1.0 );
+            sin_uu, cos_uu = sin(0.5*pi*uu), cos(0.5*pi*uu)
+            si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled)
+            si2 = sas_sinx_x(mu_proj * cos_uu)
+            inner_total += Gauss76Wt[j] * square(si1 * si2)
+        inner_total *= 0.5
+
+        si = sas_sinx_x(mu * c_scaled * sigma)
+        outer_total += Gauss76Wt[i] * inner_total * si * si
+    outer_total *= 0.5
+
+    # Multiply by contrast^2 and convert from [1e-12 A-1] to [cm-1]
+    V = form_volume(length_a, length_b, length_c)
+    drho = (sld-solvent_sld)
+    return 1.0e-4 * square(drho * V) * outer_total
+Iq.vectorized = True
+
+def Iqxy(qx, qy,
+        sld, solvent_sld,
+        length_a, length_b, length_c,
+        theta, phi, psi):
+    q, xhat, yhat, zhat = ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi)
+    siA = sas_sinx_x(0.5*length_a*q*xhat)
+    siB = sas_sinx_x(0.5*length_b*q*yhat)
+    siC = sas_sinx_x(0.5*length_c*q*zhat)
+    V = form_volume(length_a, length_b, length_c)
+    drho = (sld - solvent_sld)
+    form = V * drho * siA * siB * siC
+    # Square and convert from [1e-12 A-1] to [cm-1]
+    return 1.0e-4 * form * form
+Iqxy.vectorized = True
+
+
+def ER(length_a, length_b, length_c):
+    """
+    Return effective radius (ER) for P(q)*S(q)
+    """
+    # now that axes can be in any size order, need to sort a,b,c where a~b and c is either much smaller
+    # or much larger
+    abc = np.vstack((length_a, length_b, length_c))
+    abc = np.sort(abc, axis=0)
+    selector = (abc[1] - abc[0]) > (abc[2] - abc[1])
+    length = np.where(selector, abc[0], abc[2])
+    # surface average radius (rough approximation)
+    radius = np.sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi)
+
+    ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius))
+    return 0.5 * (ddd) ** (1. / 3.)
+
+# VR defaults to 1.0
+
+# parameters for demo
+demo = dict(scale=1, background=0,
+            sld=6.3, sld_solvent=1.0,
+            length_a=35, length_b=75, length_c=400,
+            theta=45, phi=30, psi=15,
+            length_a_pd=0.1, length_a_pd_n=10,
+            length_b_pd=0.1, length_b_pd_n=1,
+            length_c_pd=0.1, length_c_pd_n=1,
+            theta_pd=10, theta_pd_n=1,
+            phi_pd=10, phi_pd_n=1,
+            psi_pd=10, psi_pd_n=10)
+# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models
+qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
+tests = [[{}, 0.2, 0.17758004974],
+         [{}, [0.2], [0.17758004974]],
+         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0089517140475],
+         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0089517140475]],
+        ]
+del qx, qy  # not necessary to delete, but cleaner
Index: sasmodels/special.py
===================================================================
--- sasmodels/special.py	(revision 597878a0f8b1b1b38cfce28323dfa620bdbaea50)
+++ sasmodels/special.py	(revision 597878a0f8b1b1b38cfce28323dfa620bdbaea50)
@@ -0,0 +1,805 @@
+"""
+Special Functions
+.................
+
+The C code follows the C99 standard, with the usual math functions,
+as defined in
+`OpenCL <https://www.khronos.org/registry/cl/sdk/1.1/docs/man/xhtml/mathFunctions.html>`_.
+This includes the following:
+
+    M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E:
+        $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$
+    exp, log, pow(x,y), expm1, sqrt:
+        Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\sqrt{x}$.
+        The function expm1(x) is accurate across all $x$, including $x$
+        very close to zero.
+    sin, cos, tan, asin, acos, atan:
+        Trigonometry functions and inverses, operating on radians.
+    sinh, cosh, tanh, asinh, acosh, atanh:
+        Hyperbolic trigonometry functions.
+    atan2(y,x):
+        Angle from the $x$\ -axis to the point $(x,y)$, which is equal to
+        $\tan^{-1}(y/x)$ corrected for quadrant.  That is, if $x$ and $y$ are
+        both negative, then atan2(y,x) returns a value in quadrant III where
+        atan(y/x) would return a value in quadrant I. Similarly for
+        quadrants II and IV when $x$ and $y$ have opposite sign.
+    fmin(x,y), fmax(x,y), trunc, rint:
+        Floating point functions.  rint(x) returns the nearest integer.
+    NAN:
+        NaN, Not a Number, $0/0$.  Use isnan(x) to test for NaN.  Note that
+        you cannot use :code:`x == NAN` to test for NaN values since that
+        will always return false.  NAN does not equal NAN!
+    INFINITY:
+        $\infty, 1/0$.  Use isinf(x) to test for infinity, or isfinite(x)
+        to test for finite and not NaN.
+    erf, erfc, tgamma, lgamma:  **do not use**
+        Special functions that should be part of the standard, but are missing
+        or inaccurate on some platforms. Use sas_erf, sas_erfc and sas_gamma
+        instead (see below). Note: lgamma(x) has not yet been tested.
+
+Some non-standard constants and functions are also provided:
+
+    M_PI_180, M_4PI_3:
+        $\frac{\pi}{180}$, $\frac{4\pi}{3}$
+    SINCOS(x, s, c):
+        Macro which sets s=sin(x) and c=cos(x). The variables *c* and *s*
+        must be declared first.
+    square(x):
+        $x^2$
+    cube(x):
+        $x^3$
+    sas_sinx_x(x):
+        $\sin(x)/x$, with limit $\sin(0)/0 = 1$.
+    powr(x, y):
+        $x^y$ for $x \ge 0$; this is faster than general $x^y$ on some GPUs.
+    pown(x, n):
+        $x^n$ for $n$ integer; this is faster than general $x^n$ on some GPUs.
+    FLOAT_SIZE:
+        The number of bytes in a floating point value.  Even though all
+        variables are declared double, they may be converted to single
+        precision float before running. If your algorithm depends on
+        precision (which is not uncommon for numerical algorithms), use
+        the following::
+
+            #if FLOAT_SIZE>4
+            ... code for double precision ...
+            #else
+            ... code for single precision ...
+            #endif
+    SAS_DOUBLE:
+        A replacement for :code:`double` so that the declared variable will
+        stay double precision; this should generally not be used since some
+        graphics cards do not support double precision.  There is no provision
+        for forcing a constant to stay double precision.
+
+The following special functions and scattering calculations are defined in
+`sasmodels/models/lib <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib>`_.
+These functions have been tuned to be fast and numerically stable down
+to $q=0$ even in single precision.  In some cases they work around bugs
+which appear on some platforms but not others, so use them where needed.
+Add the files listed in :code:`source = ["lib/file.c", ...]` to your *model.py*
+file in the order given, otherwise these functions will not be available.
+
+    polevl(x, c, n):
+        Polynomial evaluation $p(x) = \sum_{i=0}^n c_i x^i$ using Horner's
+        method so it is faster and more accurate.
+
+        $c = \{c_n, c_{n-1}, \ldots, c_0 \}$ is the table of coefficients,
+        sorted from highest to lowest.
+
+        :code:`source = ["lib/polevl.c", ...]` (`link to code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/polevl.c>`_)
+
+    p1evl(x, c, n):
+        Evaluation of normalized polynomial $p(x) = x^n + \sum_{i=0}^{n-1} c_i x^i$
+        using Horner's method so it is faster and more accurate.
+
+        $c = \{c_{n-1}, c_{n-2} \ldots, c_0 \}$ is the table of coefficients,
+        sorted from highest to lowest.
+
+        :code:`source = ["lib/polevl.c", ...]`
+        (`link to code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/polevl.c>`_)
+
+    sas_gamma(x):
+        Gamma function $\text{sas_gamma}(x) = \Gamma(x)$.
+
+        The standard math function, tgamma(x) is unstable for $x < 1$
+        on some platforms.
+
+        :code:`source = ["lib/sasgamma.c", ...]`
+        (`link to code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gamma.c>`_)
+
+    sas_erf(x), sas_erfc(x):
+        Error function
+        $\text{sas_erf}(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt$
+        and complementary error function
+        $\text{sas_erfc}(x) = \frac{2}{\sqrt\pi}\int_x^{\infty} e^{-t^2}\,dt$.
+
+        The standard math functions erf(x) and erfc(x) are slower and broken
+        on some platforms.
+
+        :code:`source = ["lib/polevl.c", "lib/sas_erf.c", ...]`
+        (`link to error functions' code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_erf.c>`_)
+
+    sas_J0(x):
+        Bessel function of the first kind $\text{sas_J0}(x)=J_0(x)$ where
+        $J_0(x) = \frac{1}{\pi}\int_0^\pi \cos(x\sin(\tau))\,d\tau$.
+
+        The standard math function j0(x) is not available on all platforms.
+
+        :code:`source = ["lib/polevl.c", "lib/sas_J0.c", ...]`
+        (`link to Bessel function's code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_J0.c>`_)
+
+    sas_J1(x):
+        Bessel function of the first kind  $\text{sas_J1}(x)=J_1(x)$ where
+        $J_1(x) = \frac{1}{\pi}\int_0^\pi \cos(\tau - x\sin(\tau))\,d\tau$.
+
+        The standard math function j1(x) is not available on all platforms.
+
+        :code:`source = ["lib/polevl.c", "lib/sas_J1.c", ...]`
+        (`link to Bessel function's code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_J1.c>`_)
+
+    sas_JN(n, x):
+        Bessel function of the first kind and integer order $n$:
+        $\text{sas_JN}(n, x)=J_n(x)$ where
+        $J_n(x) = \frac{1}{\pi}\int_0^\pi \cos(n\tau - x\sin(\tau))\,d\tau$.
+        If $n$ = 0 or 1, it uses sas_J0(x) or sas_J1(x), respectively.
+
+        The standard math function jn(n, x) is not available on all platforms.
+
+        :code:`source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c", ...]`
+        (`link to Bessel function's code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_JN.c>`_)
+
+    sas_Si(x):
+        Sine integral $\text{Si}(x) = \int_0^x \tfrac{\sin t}{t}\,dt$.
+
+        This function uses Taylor series for small and large arguments:
+
+        For large arguments,
+
+        .. math::
+
+             \text{Si}(x) \sim \frac{\pi}{2}
+             - \frac{\cos(x)}{x}\left(1 - \frac{2!}{x^2} + \frac{4!}{x^4} - \frac{6!}{x^6} \right)
+             - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right)
+
+        For small arguments,
+
+        .. math::
+
+           \text{Si}(x) \sim x
+           - \frac{x^3}{3\times 3!} + \frac{x^5}{5 \times 5!} - \frac{x^7}{7 \times 7!}
+           + \frac{x^9}{9\times 9!} - \frac{x^{11}}{11\times 11!}
+
+        :code:`source = ["lib/Si.c", ...]`
+        (`link to code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/Si.c>`_)
+
+    sas_3j1x_x(x):
+        Spherical Bessel form
+        $\text{sph_j1c}(x) = 3 j_1(x)/x = 3 (\sin(x) - x \cos(x))/x^3$,
+        with a limiting value of 1 at $x=0$, where $j_1(x)$ is the spherical
+        Bessel function of the first kind and first order.
+
+        This function uses a Taylor series for small $x$ for numerical accuracy.
+
+        :code:`source = ["lib/sas_3j1x_x.c", ...]`
+        (`link to code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_3j1x_x.c>`_)
+
+
+    sas_2J1x_x(x):
+        Bessel form $\text{sas_J1c}(x) = 2 J_1(x)/x$, with a limiting value
+        of 1 at $x=0$, where $J_1(x)$ is the Bessel function of first kind
+        and first order.
+
+        :code:`source = ["lib/polevl.c", "lib/sas_J1.c", ...]`
+        (`link to Bessel form's code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_J1.c>`_)
+
+
+    Gauss76Z[i], Gauss76Wt[i]:
+        Points $z_i$ and weights $w_i$ for 76-point Gaussian quadrature, respectively,
+        computing $\int_{-1}^1 f(z)\,dz \approx \sum_{i=1}^{76} w_i\,f(z_i)$.
+
+        Similar arrays are available in :code:`gauss20.c` for 20-point
+        quadrature and in :code:`gauss150.c` for 150-point quadrature.
+
+        :code:`source = ["lib/gauss76.c", ...]`
+        (`link to code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/gauss76.c>`_)
+
+
+"""
+import numpy as np
+import scipy.special
+
+# Functions to add to our standard set
+from numpy import degrees, radians
+
+# C99 standard math library functions
+M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E = np.pi, np.pi/2, np.pi/4, np.sqrt(0.5), np.e
+from numpy import exp, log, power as pow, expm1, sqrt
+from numpy import sin, cos, tan, arcsin as asin, arccos as acos, arctan as atan
+from numpy import sinh, cosh, tanh, arcsinh as asinh, arccosh as acosh, arctanh as atanh
+from numpy import arctan2 as atan2
+from numpy import fmin, fmax, trunc, rint
+from numpy import NAN, inf as INFINITY
+
+# erf, erfc, tgamma, lgamma  **do not use**
+
+# Rotations of q
+def ORIENT_SYMMETRIC(qx, qy, theta, phi):
+    q = sqrt(qx*qx + qy*qy)
+    q = sqrt(qx*qx + qy*qy)
+    sin_phi, cos_phi = sin(radians(phi)), cos(radians(phi))
+    cn = (cn*qx + sn*qy)/q * sin(radians(theta))
+    cn[q==0.] = 1.
+    sn = sqrt(1 - cn**2)
+    return q, cn, sn
+
+def ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi):
+    q = sqrt(qx*qx + qy*qy)
+    qxhat = qx/q
+    qyhat = qy/q
+    sin_theta, cos_theta = sin(radians(theta)), cos(radians(theta))
+    sin_phi, cos_phi = sin(radians(phi)), cos(radians(phi))
+    sin_psi, cos_psi = sin(radians(psi)), cos(radians(psi))
+
+    xhat = (qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi)
+          + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi))
+    yhat = (qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi)
+          + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi))
+    zhat = (qxhat*(-sin_theta*cos_phi)
+          + qyhat*(-sin_theta*sin_phi))
+    return q, xhat, yhat, zhat
+
+
+# non-standard constants and functions
+
+M_PI_180, M_4PI_3 = M_PI/180, 4*M_PI/3
+
+# can't do SINCOS in python; use "s, c = SINCOS(x)" instead
+def SINCOS(x): return sin(x), cos(x)
+
+def square(x): return x*x
+
+def cube(x): return x*x*x
+
+from numpy import sinc as _sinc
+def sas_sinx_x(x): return _sinc(x/M_PI)
+def powr(x, y): return x**y
+def pown(x, n): return x**n
+
+FLOAT_SIZE = 8
+
+def polevl(x, c, n): return np.polyval(c[:n], x)
+def p1evl(x, c, n): return np.polyval(np.hstack(([1.], c))[:n], x)
+
+from scipy.special import gamma as sas_gamma
+from scipy.special import erf as sas_erf
+from scipy.special import erfc as sas_erfc
+from scipy.special import j0 as sas_J0
+from scipy.special import j1 as sas_J1
+from scipy.special import jn as sas_JN
+
+# missing sas_Si
+def sas_Si(x):
+    return scipy.special.sici(x)[0]
+
+def sas_3j1x_x(x):
+    retvalue = 3*sas_j1(x)/x
+    retvalue[x==0] = 1.
+    
+def sas_2J1x_x(x):
+    retvalue = 2*sas_J1(x)/x
+    retvalue[x==0] = 1.
+
+
+# Gaussians
+
+Gauss20Wt = [
+	.0176140071391521,
+	.0406014298003869,
+	.0626720483341091,
+	.0832767415767047,
+	.10193011981724,
+	.118194531961518,
+	.131688638449177,
+	.142096109318382,
+	.149172986472604,
+	.152753387130726,
+	.152753387130726,
+	.149172986472604,
+	.142096109318382,
+	.131688638449177,
+	.118194531961518,
+	.10193011981724,
+	.0832767415767047,
+	.0626720483341091,
+	.0406014298003869,
+	.0176140071391521
+]
+
+Gauss20Z = [
+	-.993128599185095,
+	-.963971927277914,
+	-.912234428251326,
+	-.839116971822219,
+	-.746331906460151,
+	-.636053680726515,
+	-.510867001950827,
+	-.37370608871542,
+	-.227785851141645,
+	-.076526521133497,
+	.0765265211334973,
+	.227785851141645,
+	.37370608871542,
+	.510867001950827,
+	.636053680726515,
+	.746331906460151,
+	.839116971822219,
+	.912234428251326,
+	.963971927277914,
+	.993128599185095
+]
+
+Gauss76Wt = [
+	.00126779163408536,		#0
+	.00294910295364247,
+	.00462793522803742,
+	.00629918049732845,
+	.00795984747723973,
+	.00960710541471375,
+	.0112381685696677,
+	.0128502838475101,
+	.0144407317482767,
+	.0160068299122486,
+	.0175459372914742,		#10
+	.0190554584671906,
+	.020532847967908,
+	.0219756145344162,
+	.0233813253070112,
+	.0247476099206597,
+	.026072164497986,
+	.0273527555318275,
+	.028587223650054,
+	.029773487255905,
+	.0309095460374916,		#20
+	.0319934843404216,
+	.0330234743977917,
+	.0339977794120564,
+	.0349147564835508,
+	.0357728593807139,
+	.0365706411473296,
+	.0373067565423816,
+	.0379799643084053,
+	.0385891292645067,
+	.0391332242205184,		#30
+	.0396113317090621,
+	.0400226455325968,
+	.040366472122844,
+	.0406422317102947,
+	.0408494593018285,
+	.040987805464794,
+	.0410570369162294,
+	.0410570369162294,
+	.040987805464794,
+	.0408494593018285,		#40
+	.0406422317102947,
+	.040366472122844,
+	.0400226455325968,
+	.0396113317090621,
+	.0391332242205184,
+	.0385891292645067,
+	.0379799643084053,
+	.0373067565423816,
+	.0365706411473296,
+	.0357728593807139,		#50
+	.0349147564835508,
+	.0339977794120564,
+	.0330234743977917,
+	.0319934843404216,
+	.0309095460374916,
+	.029773487255905,
+	.028587223650054,
+	.0273527555318275,
+	.026072164497986,
+	.0247476099206597,		#60
+	.0233813253070112,
+	.0219756145344162,
+	.020532847967908,
+	.0190554584671906,
+	.0175459372914742,
+	.0160068299122486,
+	.0144407317482767,
+	.0128502838475101,
+	.0112381685696677,
+	.00960710541471375,		#70
+	.00795984747723973,
+	.00629918049732845,
+	.00462793522803742,
+	.00294910295364247,
+	.00126779163408536		#75 (indexed from 0)
+]
+
+Gauss76Z = [
+	-.999505948362153,		#0
+	-.997397786355355,
+	-.993608772723527,
+	-.988144453359837,
+	-.981013938975656,
+	-.972229228520377,
+	-.961805126758768,
+	-.949759207710896,
+	-.936111781934811,
+	-.92088586125215,
+	-.904107119545567,		#10
+	-.885803849292083,
+	-.866006913771982,
+	-.844749694983342,
+	-.822068037328975,
+	-.7980001871612,
+	-.77258672828181,
+	-.74587051350361,
+	-.717896592387704,
+	-.688712135277641,
+	-.658366353758143,		#20
+	-.626910417672267,
+	-.594397368836793,
+	-.560882031601237,
+	-.526420920401243,
+	-.491072144462194,
+	-.454895309813726,
+	-.417951418780327,
+	-.380302767117504,
+	-.342012838966962,
+	-.303146199807908,		#30
+	-.263768387584994,
+	-.223945802196474,
+	-.183745593528914,
+	-.143235548227268,
+	-.102483975391227,
+	-.0615595913906112,
+	-.0205314039939986,
+	.0205314039939986,
+	.0615595913906112,
+	.102483975391227,			#40
+	.143235548227268,
+	.183745593528914,
+	.223945802196474,
+	.263768387584994,
+	.303146199807908,
+	.342012838966962,
+	.380302767117504,
+	.417951418780327,
+	.454895309813726,
+	.491072144462194,		#50
+	.526420920401243,
+	.560882031601237,
+	.594397368836793,
+	.626910417672267,
+	.658366353758143,
+	.688712135277641,
+	.717896592387704,
+	.74587051350361,
+	.77258672828181,
+	.7980001871612,	#60
+	.822068037328975,
+	.844749694983342,
+	.866006913771982,
+	.885803849292083,
+	.904107119545567,
+	.92088586125215,
+	.936111781934811,
+	.949759207710896,
+	.961805126758768,
+	.972229228520377,		#70
+	.981013938975656,
+	.988144453359837,
+	.993608772723527,
+	.997397786355355,
+	.999505948362153		#75
+]
+
+Gauss150Z = [
+  	-0.9998723404457334,
+  	-0.9993274305065947,
+  	-0.9983473449340834,
+  	-0.9969322929775997,
+  	-0.9950828645255290,
+  	-0.9927998590434373,
+  	-0.9900842691660192,
+  	-0.9869372772712794,
+  	-0.9833602541697529,
+  	-0.9793547582425894,
+  	-0.9749225346595943,
+  	-0.9700655145738374,
+  	-0.9647858142586956,
+  	-0.9590857341746905,
+  	-0.9529677579610971,
+  	-0.9464345513503147,
+  	-0.9394889610042837,
+  	-0.9321340132728527,
+  	-0.9243729128743136,
+  	-0.9162090414984952,
+  	-0.9076459563329236,
+  	-0.8986873885126239,
+  	-0.8893372414942055,
+  	-0.8795995893549102,
+  	-0.8694786750173527,
+  	-0.8589789084007133,
+  	-0.8481048644991847,
+  	-0.8368612813885015,
+  	-0.8252530581614230,
+  	-0.8132852527930605,
+  	-0.8009630799369827,
+  	-0.7882919086530552,
+  	-0.7752772600680049,
+  	-0.7619248049697269,
+  	-0.7482403613363824,
+  	-0.7342298918013638,
+  	-0.7198995010552305,
+  	-0.7052554331857488,
+  	-0.6903040689571928,
+  	-0.6750519230300931,
+  	-0.6595056411226444,
+  	-0.6436719971150083,
+  	-0.6275578900977726,
+  	-0.6111703413658551,
+  	-0.5945164913591590,
+  	-0.5776035965513142,
+  	-0.5604390262878617,
+  	-0.5430302595752546,
+  	-0.5253848818220803,
+  	-0.5075105815339176,
+  	-0.4894151469632753,
+  	-0.4711064627160663,
+  	-0.4525925063160997,
+  	-0.4338813447290861,
+  	-0.4149811308476706,
+  	-0.3959000999390257,
+  	-0.3766465660565522,
+  	-0.3572289184172501,
+  	-0.3376556177463400,
+  	-0.3179351925907259,
+  	-0.2980762356029071,
+  	-0.2780873997969574,
+  	-0.2579773947782034,
+  	-0.2377549829482451,
+  	-0.2174289756869712,
+  	-0.1970082295132342,
+  	-0.1765016422258567,
+  	-0.1559181490266516,
+  	-0.1352667186271445,
+  	-0.1145563493406956,
+  	-0.0937960651617229,
+  	-0.0729949118337358,
+  	-0.0521619529078925,
+  	-0.0313062657937972,
+  	-0.0104369378042598,
+  	0.0104369378042598,
+  	0.0313062657937972,
+  	0.0521619529078925,
+  	0.0729949118337358,
+  	0.0937960651617229,
+  	0.1145563493406956,
+  	0.1352667186271445,
+  	0.1559181490266516,
+  	0.1765016422258567,
+  	0.1970082295132342,
+  	0.2174289756869712,
+  	0.2377549829482451,
+  	0.2579773947782034,
+  	0.2780873997969574,
+  	0.2980762356029071,
+  	0.3179351925907259,
+  	0.3376556177463400,
+  	0.3572289184172501,
+  	0.3766465660565522,
+  	0.3959000999390257,
+  	0.4149811308476706,
+  	0.4338813447290861,
+  	0.4525925063160997,
+  	0.4711064627160663,
+  	0.4894151469632753,
+  	0.5075105815339176,
+  	0.5253848818220803,
+  	0.5430302595752546,
+  	0.5604390262878617,
+  	0.5776035965513142,
+  	0.5945164913591590,
+  	0.6111703413658551,
+  	0.6275578900977726,
+  	0.6436719971150083,
+  	0.6595056411226444,
+  	0.6750519230300931,
+  	0.6903040689571928,
+  	0.7052554331857488,
+  	0.7198995010552305,
+  	0.7342298918013638,
+  	0.7482403613363824,
+  	0.7619248049697269,
+  	0.7752772600680049,
+  	0.7882919086530552,
+  	0.8009630799369827,
+  	0.8132852527930605,
+  	0.8252530581614230,
+  	0.8368612813885015,
+  	0.8481048644991847,
+  	0.8589789084007133,
+  	0.8694786750173527,
+  	0.8795995893549102,
+  	0.8893372414942055,
+  	0.8986873885126239,
+  	0.9076459563329236,
+  	0.9162090414984952,
+  	0.9243729128743136,
+  	0.9321340132728527,
+  	0.9394889610042837,
+  	0.9464345513503147,
+  	0.9529677579610971,
+  	0.9590857341746905,
+  	0.9647858142586956,
+  	0.9700655145738374,
+  	0.9749225346595943,
+  	0.9793547582425894,
+  	0.9833602541697529,
+  	0.9869372772712794,
+  	0.9900842691660192,
+  	0.9927998590434373,
+  	0.9950828645255290,
+  	0.9969322929775997,
+  	0.9983473449340834,
+  	0.9993274305065947,
+  	0.9998723404457334
+]
+
+Gauss150Wt = [
+  	0.0003276086705538,
+  	0.0007624720924706,
+  	0.0011976474864367,
+  	0.0016323569986067,
+  	0.0020663664924131,
+  	0.0024994789888943,
+  	0.0029315036836558,
+  	0.0033622516236779,
+  	0.0037915348363451,
+  	0.0042191661429919,
+  	0.0046449591497966,
+  	0.0050687282939456,
+  	0.0054902889094487,
+  	0.0059094573005900,
+  	0.0063260508184704,
+  	0.0067398879387430,
+  	0.0071507883396855,
+  	0.0075585729801782,
+  	0.0079630641773633,
+  	0.0083640856838475,
+  	0.0087614627643580,
+  	0.0091550222717888,
+  	0.0095445927225849,
+  	0.0099300043714212,
+  	0.0103110892851360,
+  	0.0106876814158841,
+  	0.0110596166734735,
+  	0.0114267329968529,
+  	0.0117888704247183,
+  	0.0121458711652067,
+  	0.0124975796646449,
+  	0.0128438426753249,
+  	0.0131845093222756,
+  	0.0135194311690004,
+  	0.0138484622795371,
+  	0.0141714592928592,
+  	0.0144882814685445,
+  	0.0147987907597169,
+  	0.0151028518701744,
+  	0.0154003323133401,
+  	0.0156911024699895,
+  	0.0159750356447283,
+  	0.0162520081211971,
+  	0.0165218992159766,
+  	0.0167845913311726,
+  	0.0170399700056559,
+  	0.0172879239649355,
+  	0.0175283451696437,
+  	0.0177611288626114,
+  	0.0179861736145128,
+  	0.0182033813680609,
+  	0.0184126574807331,
+  	0.0186139107660094,
+  	0.0188070535331042,
+  	0.0189920016251754,
+  	0.0191686744559934,
+  	0.0193369950450545,
+  	0.0194968900511231,
+  	0.0196482898041878,
+  	0.0197911283358190,
+  	0.0199253434079123,
+  	0.0200508765398072,
+  	0.0201676730337687,
+  	0.0202756819988200,
+  	0.0203748563729175,
+  	0.0204651529434560,
+  	0.0205465323660984,
+  	0.0206189591819181,
+  	0.0206824018328499,
+  	0.0207368326754401,
+  	0.0207822279928917,
+  	0.0208185680053983,
+  	0.0208458368787627,
+  	0.0208640227312962,
+  	0.0208731176389954,
+  	0.0208731176389954,
+  	0.0208640227312962,
+  	0.0208458368787627,
+  	0.0208185680053983,
+  	0.0207822279928917,
+  	0.0207368326754401,
+  	0.0206824018328499,
+  	0.0206189591819181,
+  	0.0205465323660984,
+  	0.0204651529434560,
+  	0.0203748563729175,
+  	0.0202756819988200,
+  	0.0201676730337687,
+  	0.0200508765398072,
+  	0.0199253434079123,
+  	0.0197911283358190,
+  	0.0196482898041878,
+  	0.0194968900511231,
+  	0.0193369950450545,
+  	0.0191686744559934,
+  	0.0189920016251754,
+  	0.0188070535331042,
+  	0.0186139107660094,
+  	0.0184126574807331,
+  	0.0182033813680609,
+  	0.0179861736145128,
+  	0.0177611288626114,
+  	0.0175283451696437,
+  	0.0172879239649355,
+  	0.0170399700056559,
+  	0.0167845913311726,
+  	0.0165218992159766,
+  	0.0162520081211971,
+  	0.0159750356447283,
+  	0.0156911024699895,
+  	0.0154003323133401,
+  	0.0151028518701744,
+  	0.0147987907597169,
+  	0.0144882814685445,
+  	0.0141714592928592,
+  	0.0138484622795371,
+  	0.0135194311690004,
+  	0.0131845093222756,
+  	0.0128438426753249,
+  	0.0124975796646449,
+  	0.0121458711652067,
+  	0.0117888704247183,
+  	0.0114267329968529,
+  	0.0110596166734735,
+  	0.0106876814158841,
+  	0.0103110892851360,
+  	0.0099300043714212,
+  	0.0095445927225849,
+  	0.0091550222717888,
+  	0.0087614627643580,
+  	0.0083640856838475,
+  	0.0079630641773633,
+  	0.0075585729801782,
+  	0.0071507883396855,
+  	0.0067398879387430,
+  	0.0063260508184704,
+  	0.0059094573005900,
+  	0.0054902889094487,
+  	0.0050687282939456,
+  	0.0046449591497966,
+  	0.0042191661429919,
+  	0.0037915348363451,
+  	0.0033622516236779,
+  	0.0029315036836558,
+  	0.0024994789888943,
+  	0.0020663664924131,
+  	0.0016323569986067,
+  	0.0011976474864367,
+  	0.0007624720924706,
+  	0.0003276086705538
+]
+
