Index: sasmodels/models/HayterMSAsq.py
===================================================================
--- sasmodels/models/HayterMSAsq.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
+++ sasmodels/models/HayterMSAsq.py	(revision ec2ca99b9e3d90e4b91919c523a94d9855ca9c3c)
@@ -52,23 +52,22 @@
     [Hayter-Penfold MSA charged sphere interparticle S(Q) structure factor]
         Interparticle structure factor S(Q)for a charged hard spheres.
-        Routine takes absolute value of charge, use HardSphere if charge goes to zero.
+        Routine takes absolute value of charge, use HardSphere if charge
+        goes to zero.
         In sasview the effective radius will be calculated from the
         parameters used in P(Q).
 """
 single = False  # double precision only for now
+
+# pylint: disable=bad-whitespace, line-too-long
 #             [ "name", "units", default, [lower, upper], "type", "description" ],
-parameters = [["effect_radius", "Ang", 20.75, [0, inf], "volume",
-               "effective radius of hard sphere"],
-              ["charge", "e", 19.0, [0, inf], "",
-               "charge on sphere (in electrons)"],
-              ["volfraction", "", 0.0192, [0, 0.74], "",
-               "volume fraction of spheres"],
-              ["temperature", "K", 318.16, [0, inf], "",
-               "temperature, in Kelvin, for Debye length calculation"],
-              ["saltconc", "M", 0.0, [-inf, inf], "",
-               "conc of salt, 1:1 electolyte, for Debye length"],
-              ["dielectconst", "", 71.08, [-inf, inf], "",
-               "dielectric constant of solvent (default water), for Debye length"],
-             ]
+parameters = [
+    ["effect_radius", "Ang", 20.75,   [0, inf],    "volume", "effective radius of hard sphere"],
+    ["charge",        "e",   19.0,    [0, inf],    "", "charge on sphere (in electrons)"],
+    ["volfraction",   "",     0.0192, [0, 0.74],   "", "volume fraction of spheres"],
+    ["temperature",   "K",  318.16,   [0, inf],    "", "temperature, in Kelvin, for Debye length calculation"],
+    ["saltconc",      "M",    0.0,    [-inf, inf], "", "conc of salt, 1:1 electolyte, for Debye length"],
+    ["dielectconst",  "",    71.08,   [-inf, inf], "", "dielectric constant of solvent (default water), for Debye length"],
+    ]
+# pylint: enable=bad-whitespace, line-too-long
 category = "structure-factor"
 
@@ -88,12 +87,28 @@
 oldpars = dict()
 # default parameter set,  use  compare.sh -midQ -linear
-# note the calculation varies in different limiting cases so a wide range of parameters will be required for a thorough test!
+# note the calculation varies in different limiting cases so a wide range of
+# parameters will be required for a thorough test!
 # odd that the default st has saltconc zero
-demo = dict(effect_radius = 20.75,charge=19.0,volfraction = 0.0192,temperature=318.16,saltconc=0.05,dielectconst=71.08,effect_radius_pd = 0.1,effect_radius_pd_n = 40)
+demo = dict(effect_radius=20.75,
+            charge=19.0,
+            volfraction=0.0192,
+            temperature=318.16,
+            saltconc=0.05,
+            dielectconst=71.08,
+            effect_radius_pd=0.1,
+            effect_radius_pd_n=40)
 #
 # attempt to use same values as old sasview unit test
 tests = [
-        [ {'scale': 1.0, 'background' : 0.0, 'effect_radius' : 20.75, 'charge' : 19.0, 'volfraction' : 0.0192, 'temperature' : 298.0,
-          'saltconc' : 0,'dielectconst' : 78.0, 'effect_radius_pd' : 0}, [0.0010], [0.0712928]]
-        ]
+    [{'scale': 1.0,
+      'background': 0.0,
+      'effect_radius': 20.75,
+      'charge': 19.0,
+      'volfraction': 0.0192,
+      'temperature': 298.0,
+      'saltconc': 0,
+      'dielectconst': 78.0,
+      'effect_radius_pd': 0},
+     [0.0010], [0.0712928]]
+    ]
 
Index: sasmodels/models/bcc.c
===================================================================
--- sasmodels/models/bcc.c	(revision 7ed702f2d0d4ecbcb0115970cdabf0d05ac08de0)
+++ sasmodels/models/bcc.c	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
@@ -40,12 +40,4 @@
 }
 
-double _sphereform(double q, double radius, double sld, double solvent_sld){
-    const double qr = q*radius;
-    double sn, cn;
-    SINCOS(qr, sn, cn);
-    const double bes = (qr == 0.0 ? 1.0 : 3.0*(sn-qr*cn)/(qr*qr*qr));
-    const double fq = bes * (sld - solvent_sld)*form_volume(radius);
-    return 1.0e-4*fq*fq;
-}
 
 double form_volume(double radius){
@@ -87,5 +79,5 @@
 
 	answer = (vb-va)/2.0*summ;
-	answer = answer*_sphereform(q,radius,sld,solvent_sld)*latticescale;
+	answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale;
 
     return answer;
@@ -174,5 +166,5 @@
 
   // Use SphereForm directly from libigor
-  answer = _sphereform(q,radius,sld,solvent_sld)*Zq*latticescale;
+  answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
 
   return answer;
Index: sasmodels/models/bcc.py
===================================================================
--- sasmodels/models/bcc.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
+++ sasmodels/models/bcc.py	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
@@ -132,5 +132,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/J1.c", "lib/gauss150.c", "bcc.c"]
+source = ["lib/J1.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc.c"]
 
 # parameters for demo
Index: sasmodels/models/core_shell_sphere.c
===================================================================
--- sasmodels/models/core_shell_sphere.c	(revision edc9f8d2c5449514c0206c9ed35db99d865b35b6)
+++ sasmodels/models/core_shell_sphere.c	(revision edc9f8d2c5449514c0206c9ed35db99d865b35b6)
@@ -0,0 +1,31 @@
+double form_volume(double radius, double thickness);
+double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld);
+double Iqxy(double qx, double qy, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld);
+
+
+double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld) {
+    // Core first, then add in shell
+    const double core_qr = q * radius;
+    const double core_contrast = core_sld - shell_sld;
+    const double core_bes = sph_j1c(core_qr);
+    const double core_volume = 4.0 * M_PI / 3.0 * radius * radius * radius;
+    double f = core_volume * core_bes * core_contrast;
+
+    // Now the shell
+    const double shell_qr = q * (radius + thickness);
+    const double shell_contrast = shell_sld - solvent_sld;
+    const double shell_bes = sph_j1c(shell_qr);
+    const double shell_volume = 4.0 * M_PI / 3.0 * pow((radius + thickness), 3);
+    f += shell_volume * shell_bes * shell_contrast;
+    return f * f * 1.0e-4;
+}
+
+double Iqxy(double qx, double qy, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld) {
+    const double q = sqrt(qx*qx+qy*qy);
+    return Iq(q, radius, thickness, core_sld, shell_sld, solvent_sld);
+}
+
+double form_volume(double radius, double thickness)
+{
+    return 4.0 * M_PI / 3.0 * pow((radius + thickness), 3);
+}
Index: sasmodels/models/core_shell_sphere.py
===================================================================
--- sasmodels/models/core_shell_sphere.py	(revision 8c9dbc9161c033ba16aba0b54213deb314c953a0)
+++ sasmodels/models/core_shell_sphere.py	(revision 8c9dbc9161c033ba16aba0b54213deb314c953a0)
@@ -0,0 +1,110 @@
+r"""
+This model provides the form factor, $P(q)$, for a spherical particle with a core-shell structure.
+The form factor is normalized by the particle volume.
+
+Definition
+----------
+
+The 1D scattering intensity is calculated in the following way (Guinier, 1955)
+
+.. math::
+
+    P(q) = \frac{\text{scale}}{V} F^2(q) + \text{background}
+
+where
+
+.. math::
+    F^2(q)=\frac{3}{V_s}\left[V_c(\rho_c-\rho_s)\frac{\sin(qr_c)-qr_c\cos(qr_c)}{(qr_c)^3}+
+    V_s(\rho_s-\rho_{solv})\frac{\sin(qr_s)-qr_s\cos(qr_s)}{(qr_s)^3}\right]
+
+where $V_s$ is the volume of the outer shell, $V_c$ is
+the volume of the core, $r_s$ is the radius of the shell, $r_c$ is the radius of the
+core, $\rho_c$ is the scattering length density of the core, $\rho_s$ is the scattering length
+density of the shell, $\rho_{solv}$ is the scattering length density of the solvent.
+
+The 2D scattering intensity is the same as $P(q)$ above, regardless of the
+orientation of the $q$ vector.
+
+NB: The outer most radius (ie, = radius + thickness) is used as the
+effective radius for $S(Q)$ when $P(Q) \cdot S(Q)$ is applied.
+
+Reference
+---------
+
+A Guinier and G Fournet, *Small-Angle Scattering of X-Rays*, John Wiley and Sons, New York, (1955)
+
+Validation
+----------
+
+Validation of our code was done by comparing the output of the 1D model to the output of
+the software provided by NIST (Kline, 2006). Figure 1 shows a comparison of the output of
+our model and the output of the NIST software.
+
+.. image:: img/core_shell_sphere_1d.jpg
+
+    Figure 1: Comparison of the SasView scattering intensity for a core-shell sphere with
+    the output of the NIST SANS analysis software. The parameters were set to:
+    *scale* = 1.0, *radius* = 60 , *contrast* = 1e-6 |Ang^-2|, and
+    *background* = 0.001 |cm^-1|.
+"""
+
+from numpy import pi, inf
+
+name = "core_shell_sphere"
+title = "Form factor for a monodisperse spherical particle with particle with a core-shell structure."
+description = """
+    F^2(q) = 3/V_s [V_c (core_sld-shell_sld) (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3 
+                   + V_s (shell_sld-solvent_sld) (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3]
+
+            V_s: Volume of the sphere shell
+            V_c: Volume of the sphere core
+            r_s: Shell radius = radius + thickness
+"""
+category = "shape:sphere"
+
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["radius",      "Ang",        60.0, [0, inf],    "volume", "Sphere core radius"],
+              ["thickness",   "Ang",        10.0, [0, inf],    "volume", "Sphere shell thickness"],
+              ["core_sld",    "1e-6/Ang^2", 1.0,  [-inf, inf], "",       "Sphere core scattering length density"],
+              ["shell_sld",   "1e-6/Ang^2", 2.0,  [-inf, inf], "",       "Sphere shell scattering length density"],
+              ["solvent_sld", "1e-6/Ang^2", 3.0,  [-inf, inf],  "",      "Solvent scattering length density"]]
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lib/sph_j1c.c", "core_shell_sphere.c"]
+
+demo = dict(scale=1, background=0, radius=60, thickness=10,
+            core_sld=1.0, shell_sld=2.0, solvent_sld=0.0)
+
+oldname = 'CoreShellModel'
+oldpars = {}
+
+def ER(radius, thickness):
+    """
+        Equivalent radius
+        @param radius: core radius
+        @param thickness: shell thickness
+    """
+    return radius + thickness
+
+def VR(radius, thickness):
+    """
+        Volume ratio
+        @param radius: core radius
+        @param thickness: shell thickness
+    """
+    whole = 4.0 * pi / 3.0 * pow((radius + thickness), 3)
+    core = 4.0 * pi / 3.0 * radius * radius * radius
+    return whole, whole - core
+
+tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0],
+         [{'radius': 20.0, 'thickness': 10.0}, 'VR', 0.703703704],
+
+         # The SasView test result was 0.00169, with a background of 0.001
+         [{'radius': 60.0,
+           'thickness': 10.0,
+           'core_sld': 1.0,
+           'shell_sld':2.0,
+           'solvent_sld':3.0,
+           'background':0.0
+          }, 0.4, 0.000698838]]
Index: sasmodels/models/elliptical_cylinder.c
===================================================================
--- sasmodels/models/elliptical_cylinder.c	(revision a8b3cdb60df4f2e22c9b3c5c3110bb51b7dba2c6)
+++ sasmodels/models/elliptical_cylinder.c	(revision a8b3cdb60df4f2e22c9b3c5c3110bb51b7dba2c6)
@@ -0,0 +1,170 @@
+double form_volume(double r_minor, double r_ratio, double length);
+double Iq(double q, double r_minor, double r_ratio, double length,
+          double sld, double solvent_sld);
+double Iqxy(double qx, double qy, double r_minor, double r_ratio, double length,
+            double sld, double solvent_sld, double theta, double phi, double psi);
+
+
+double _elliptical_cylinder_kernel(double q, double r_minor, double r_ratio, double theta);
+
+double _elliptical_cylinder_kernel(double q, double r_minor, double r_ratio, double theta)
+{
+    // This is the function LAMBDA1^2 in Feigin's notation
+    // q is the q-value for the calculation (1/A)
+    // r_minor is the transformed radius"a" in Feigin's notation
+    // r_ratio is the ratio (major radius)/(minor radius) of the Ellipsoid [=] ---
+    // theta is the dummy variable of the integration
+
+    double retval,arg;
+
+    arg = q*r_minor*sqrt((1.0+r_ratio*r_ratio)/2+(1.0-r_ratio*r_ratio)*cos(theta)/2);
+    if (arg == 0.0){
+        retval = 1.0;
+    }else{
+        retval = 2.0*NR_BessJ1(arg)/arg;
+    }
+    return retval*retval ;
+}
+
+
+double form_volume(double r_minor, double r_ratio, double length)
+{
+    return M_PI * r_minor * r_minor * r_ratio * length;
+}
+
+double Iq(double q, double r_minor, double r_ratio, double length,
+          double sld, double solvent_sld) {
+
+    const int nordi=76; //order of integration
+    const int nordj=20;
+    double va,vb;       //upper and lower integration limits
+    double summ,zi,yyy,answer;         //running tally of integration
+    double summj,vaj,vbj,zij,arg,si;            //for the inner integration
+
+    // orientational average limits
+    va = 0.0;
+    vb = 1.0;
+    // inner integral limits
+    vaj=0.0;
+    vbj=M_PI;
+
+    //initialize integral
+    summ = 0.0;
+
+    const double delrho = sld - solvent_sld;
+
+    for(int i=0;i<nordi;i++) {
+        //setup inner integral over the ellipsoidal cross-section
+        summj=0;
+        zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;     //the "x" dummy
+        arg = r_minor*sqrt(1.0-zi*zi);
+        for(int j=0;j<nordj;j++) {
+            //20 gauss points for the inner integral
+            zij = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0;        //the "y" dummy
+            yyy = Gauss20Wt[j] * _elliptical_cylinder_kernel(q, arg, r_ratio, zij);
+            summj += yyy;
+        }
+        //now calculate the value of the inner integral
+        answer = (vbj-vaj)/2.0*summj;
+        //divide integral by Pi
+        answer /= M_PI;
+
+        //now calculate outer integral
+        arg = q*length*zi/2.0;
+        if (arg == 0.0){
+            si = 1.0;
+        }else{
+            si = sin(arg) * sin(arg) / arg / arg;
+        }
+        yyy = Gauss76Wt[i] * answer * si;
+        summ += yyy;
+    }
+
+    answer = (vb-va)/2.0*summ;
+    // Multiply by contrast^2
+    answer *= delrho*delrho;
+
+    const double vol = form_volume(r_minor, r_ratio, length);
+    return answer*vol*vol*1e-4;
+}
+
+
+double Iqxy(double qx, double qy, double r_minor, double r_ratio, double length,
+            double sld, double solvent_sld, double theta, double phi, double psi) {
+    const double _theta = theta * M_PI / 180.0;
+    const double _phi = phi * M_PI / 180.0;
+    const double _psi = psi * M_PI / 180.0;
+    const double q = sqrt(qx*qx+qy*qy);
+    const double q_x = qx/q;
+    const double q_y = qy/q;
+
+    //Cylinder orientation
+    double cyl_x = cos(_theta) * cos(_phi);
+    double cyl_y = sin(_theta);
+
+    //cyl_z = -cos(_theta) * sin(_phi);
+
+    // q vector
+    //q_z = 0;
+
+    // Note: cos(alpha) = 0 and 1 will get an
+    // undefined value from CylKernel
+    //alpha = acos( cos_val );
+
+    //ellipse orientation:
+    // the elliptical corss section was transformed and projected
+    // into the detector plane already through sin(alpha)and furthermore psi remains as same
+    // on the detector plane.
+    // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt
+    // the wave vector q.
+
+    //x- y- component on the detector plane.
+    const double ella_x =  -cos(_phi)*sin(_psi) * sin(_theta)+sin(_phi)*cos(_psi);
+    const double ella_y =  sin(_psi)*cos(_theta);
+    const double ellb_x =  -sin(_theta)*cos(_psi)*cos(_phi)-sin(_psi)*sin(_phi);
+    const double ellb_y =  cos(_theta)*cos(_psi);
+
+    // Compute the angle btw vector q and the
+    // axis of the cylinder
+    double cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
+
+    // calculate the axis of the ellipse wrt q-coord.
+    double cos_nu = ella_x*q_x + ella_y*q_y;
+    double cos_mu = ellb_x*q_x + ellb_y*q_y;
+
+    // The following test should always pass
+    if (fabs(cos_val)>1.0) {
+      //printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
+      cos_val = 1.0;
+    }
+    if (fabs(cos_nu)>1.0) {
+      //printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n");
+      cos_nu = 1.0;
+    }
+    if (fabs(cos_mu)>1.0) {
+      //printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n");
+      cos_mu = 1.0;
+    }
+
+    const double r_major = r_ratio * r_minor;
+    const double qr = q*sqrt( r_major*r_major*cos_nu*cos_nu + r_minor*r_minor*cos_mu*cos_mu );
+    const double qL = q*length*cos_val/2.0;
+
+    double Be;
+    if (qr==0){
+      Be = 0.5;
+    }else{
+      Be = NR_BessJ1(qr)/qr;
+    }
+
+    double Si;
+    if (qL==0){
+      Si = 1.0;
+    }else{
+      Si = sin(qL)/qL;
+    }
+
+    const double k = 2.0 * Be * Si;
+    const double vol = form_volume(r_minor, r_ratio, length);
+    return (sld - solvent_sld) * (sld - solvent_sld) * k * k *vol*vol*1.0e-4;
+}
Index: sasmodels/models/elliptical_cylinder.py
===================================================================
--- sasmodels/models/elliptical_cylinder.py	(revision b7c2fcea5b553bdc935287b253bc4a360bf6cf0c)
+++ sasmodels/models/elliptical_cylinder.py	(revision b7c2fcea5b553bdc935287b253bc4a360bf6cf0c)
@@ -0,0 +1,146 @@
+# pylint: disable=line-too-long
+r"""
+This function calculates the scattering from an elliptical cylinder.
+
+Definition for 2D (orientated system)
+-------------------------------------
+
+The angles |theta| and |phi| define the orientation of the axis of the cylinder. The angle |bigpsi| is defined as the
+orientation of the major axis of the ellipse with respect to the vector *Q*\ . A gaussian polydispersity can be added
+to any of the orientation angles, and also for the minor radius and the ratio of the ellipse radii.
+
+.. image:: img/elliptical_cylinder_geometry.gif
+
+    *Figure.* *a* = *r_minor* and |nu|\ :sub:`n` = $r_ratio$ (i.e., $r_major / r_minor$).
+
+The function calculated is
+
+.. math::
+    I(\mathbf{q})=\frac{1}{V_{cyl}}\int{d\psi}\int{d\phi}\int{p(\theta,\phi,\psi)F^2(\mathbf{q},\alpha,\psi)\sin(\theta)d\theta}
+
+with the functions
+
+.. math::
+    F(\mathbf{q},\alpha,\psi)=2\frac{J_1(a)\sin(b)}{ab}
+    \\
+    a = \mathbf{q}\sin(\alpha)\left[ r^2_{major}\sin^2(\psi)+r^2_{minor}\cos(\psi) \right]^{1/2}
+    \\
+    b=\mathbf{q}\frac{L}{2}\cos(\alpha)
+
+and the angle |bigpsi| is defined as the orientation of the major axis of the ellipse with respect to the vector $\vec q$ .
+The angle $\alpha$ is the angle between the axis of the cylinder and $\vec q$.
+
+
+Definition for 1D (no preferred orientation)
+--------------------------------------------
+
+The form factor is averaged over all possible orientation before normalized by the particle volume
+
+.. math::
+    P(q) = scale  <F^2> / V
+
+The returned value is scaled to units of |cm^-1|.
+
+To provide easy access to the orientation of the elliptical cylinder, we define the axis of the cylinder using two
+angles |theta|, |phi| and |bigpsi|. As for the case of the cylinder, the angles |theta| and |phi| are defined on
+Figure 2 of CylinderModel. The angle |bigpsi| is the rotational angle around its own long_c axis against the *q* plane.
+For example, |bigpsi| = 0 when the *r_minor* axis is parallel to the *x*\ -axis of the detector.
+
+All angle parameters are valid and given only for 2D calculation; ie, an oriented system.
+
+.. image:: img/elliptical_cylinder_geometry_2d.jpg
+
+    *Figure. Definition of angles for 2D*
+
+.. image:: img/core_shell_bicelle_fig2.jpg
+
+    *Figure. Examples of the angles for oriented elliptical cylinders against the detector plane.*
+
+NB: The 2nd virial coefficient of the cylinder is calculated based on the averaged radius (= sqrt(*r_minor*\ :sup:`2` \* *r_ratio*))
+and length values, and used as the effective radius for *S(Q)* when *P(Q)* \* *S(Q)* is applied.
+
+
+.. image:: img/elliptical_cylinder_comparison_1d.jpg
+
+    *Figure. 1D plot using the default values (w/1000 data point).*
+
+Validation
+----------
+
+Validation of our code was done by comparing the output of the 1D calculation to the angular average of the output of
+the 2D calculation over all possible angles. The figure below shows the comparison where the solid dot refers to
+averaged 2D values while the line represents the result of the 1D calculation (for the 2D averaging, values of 76, 180,
+and 76 degrees are taken for the angles of |theta|, |phi|, and |bigpsi| respectively).
+
+.. image:: img/elliptical_cylinder_validation_1d.gif
+
+    *Figure. Comparison between 1D and averaged 2D.*
+
+In the 2D average, more binning in the angle |phi| is necessary to get the proper result. The following figure shows
+the results of the averaging by varying the number of angular bins.
+
+.. image:: img/elliptical_cylinder_averaging.gif
+
+    *Figure. The intensities averaged from 2D over different numbers of bins and angles.*
+
+Reference
+---------
+
+L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum,
+New York, (1987)
+"""
+
+import math
+from numpy import pi, inf
+
+name = "elliptical_cylinder"
+title = "Form factor for an elliptical cylinder."
+description = """
+    Form factor for an elliptical cylinder.
+    See L A Feigin and D I Svergun, Structure Analysis by Small-Angle X-Ray and Neutron Scattering, Plenum, New York, (1987).
+"""
+category = "shape:cylinder"
+
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["r_minor",     "Ang",        20.0,  [0, inf],    "volume",      "Ellipse minor radius"],
+              ["r_ratio",     "",           1.5,   [1, inf],    "volume",      "Ratio of major radius over minor radius"],
+              ["length",      "Ang",        400.0, [1, inf],    "volume",      "Length of the cylinder"],
+              ["sld",         "1e-6/Ang^2", 4.0,   [-inf, inf], "",            "Cylinder scattering length density"],
+              ["solvent_sld", "1e-6/Ang^2", 1.0,   [-inf, inf], "",            "Solvent scattering length density"],
+              ["theta",       "degrees",    90.0,  [-360, 360], "orientation", "In plane angle"],
+              ["phi",         "degrees",    0,     [-360, 360], "orientation", "Out of plane angle"],
+              ["psi",         "degrees",    0,     [-360, 360], "orientation", "Major axis angle relative to Q"]]
+
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lib/nr_bess_j1.c", "lib/gauss76.c", "lib/gauss20.c", "elliptical_cylinder.c"]
+
+demo = dict(scale=1, background=0, r_minor=100, r_ratio=1.5, length=400.0,
+            sld=4.0, solvent_sld=1.0, theta=10.0, phi=20, psi=30, theta_pd=10, phi_pd=2, psi_pd=3)
+
+oldname = 'EllipticalCylinderModel'
+oldpars = dict(theta='cyl_theta', phi='cyl_phi', psi='cyl_psi', sld='sldCyl', solvent_sld='sldSolv')
+
+def ER(r_minor, r_ratio, length):
+    """
+        Equivalent radius
+        @param r_minor: Ellipse minor radius
+        @param r_ratio: Ratio of major radius over minor radius
+        @param length: Length of the cylinder
+    """
+    radius = math.sqrt(r_minor * r_minor * r_ratio)
+    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius))
+    return 0.5 * (ddd) ** (1. / 3.)
+
+tests = [[{'r_minor': 20.0, 'r_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024],
+         [{'r_minor': 20.0, 'r_ratio': 1.2, 'length':300.0}, 'VR', 1],
+
+         # The SasView test result was 0.00169, with a background of 0.001
+         [{'r_minor': 20.0,
+           'r_ratio': 1.5,
+           'sld': 4.0,
+           'length':400.0,
+           'solvent_sld':1.0,
+           'background':0.0
+          }, 0.001, 675.504402]]
Index: sasmodels/models/fcc.c
===================================================================
--- sasmodels/models/fcc.c	(revision eeb8bacf06af69341f80d5e76d7b10a7239c4e67)
+++ sasmodels/models/fcc.c	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
@@ -7,5 +7,4 @@
 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
 double _FCCeval(double Theta, double Phi, double temp1, double temp3);
-double _sphereform(double q, double radius, double sld, double solvent_sld);
 
 
@@ -39,13 +38,4 @@
 
 	return (result);
-}
-
-double _sphereform(double q, double radius, double sld, double solvent_sld){
-    const double qr = q*radius;
-    double sn, cn;
-    SINCOS(qr, sn, cn);
-    const double bes = (qr == 0.0 ? 1.0 : 3.0*(sn-qr*cn)/(qr*qr*qr));
-    const double fq = bes * (sld - solvent_sld)*form_volume(radius);
-    return 1.0e-4*fq*fq;
 }
 
@@ -88,5 +78,5 @@
 
 	answer = (vb-va)/2.0*summ;
-	answer = answer*_sphereform(q,radius,sld,solvent_sld)*latticescale;
+	answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale;
 
     return answer;
@@ -175,5 +165,5 @@
 
   // Use SphereForm directly from libigor
-  answer = _sphereform(q,radius,sld,solvent_sld)*Zq*latticescale;
+  answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
 
   return answer;
Index: sasmodels/models/fcc.py
===================================================================
--- sasmodels/models/fcc.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
+++ sasmodels/models/fcc.py	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
@@ -125,5 +125,5 @@
              ]
 
-source = ["lib/J1.c", "lib/gauss150.c", "fcc.c"]
+source = ["lib/J1.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc.c"]
 
 # parameters for demo
Index: sasmodels/models/hardsphere.py
===================================================================
--- sasmodels/models/hardsphere.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
+++ sasmodels/models/hardsphere.py	(revision 093f7545303dc8a4e5074424203e44456c438f42)
@@ -35,6 +35,6 @@
 from numpy import inf
 
-name = "hardsphere"
-title = "Hard sphere structure factor, with Percus-Yevick closure"
+name = "hardsphere_fish"
+title = "Hard sphere structure factor from FISH, with Percus-Yevick closure"
 description = """\
     [Hard sphere structure factor, with Percus-Yevick closure]
@@ -55,5 +55,4 @@
                "volume fraction of hard spheres"],
              ]
-single = False
 
 # No volume normalization despite having a volume parameter
@@ -64,33 +63,49 @@
 
 Iq = """
-    double denom,dnum,alpha,beta,gamm,a,asq,ath,afor,rca,rsa;
-    double calp,cbeta,cgam,prefac,c,vstruc;
-    double struc;
+      double D,A,B,G,X,X2,X4,S,C,FF,HARDSPH;
 
-    //  compute constants
-    denom = pow((1.0-volfraction),4);
-    dnum = pow((1.0 + 2.0*volfraction),2);
-    alpha = dnum/denom;
-    beta = -6.0*volfraction*pow((1.0 + volfraction/2.0),2)/denom;
-    gamm = 0.50*volfraction*dnum/denom;
-    //
-    //  calculate the structure factor
-    //
-    a = 2.0*q*effect_radius;
-    asq = a*a;
-    ath = asq*a;
-    afor = ath*a;
-    SINCOS(a,rsa,rca);
-    //rca = cos(a);
-    //rsa = sin(a);
-    calp = alpha*(rsa/asq - rca/a);
-    cbeta = beta*(2.0*rsa/asq - (asq - 2.0)*rca/ath - 2.0/ath);
-    cgam = gamm*(-rca/a + (4.0/a)*((3.0*asq - 6.0)*rca/afor + (asq - 6.0)*rsa/ath + 6.0/afor));
-    prefac = -24.0*volfraction/a;
-    c = prefac*(calp + cbeta + cgam);
-    vstruc = 1.0/(1.0-c);
-    struc = vstruc;
+      if(fabs(effect_radius) < 1.E-12) {
+               HARDSPH=1.0;
+               return(HARDSPH);
+      }
+      D=pow((1.-volfraction),2);
+      A=pow((1.+2.*volfraction)/D, 2);
+      X=fabs(q*effect_radius*2.0);
 
-    return(struc);
+      if(X < 5.E-06) {
+                 HARDSPH=1./A;
+                 return(HARDSPH);
+      }
+      X2=pow(X,2);
+      X4=pow(X2,2);
+      B=-6.*volfraction* pow((1.+0.5*volfraction)/D ,2);
+      G=0.5*volfraction*A;
+
+      if(X < 0.2) {
+      // use Taylor series expansion for small X, IT IS VERY PICKY ABOUT THE X CUT OFF VALUE, ought to be lower in double. 
+      // No obvious way to rearrange the equations to avoid needing a very high number of significant figures. 
+      // Series expansion found using Mathematica software. Numerical test in .xls showed terms to X^2 are sufficient 
+      // for 5 or 6 significant figures but I put the X^4 one in anyway 
+            FF = 8*A +6*B + 4*G - (0.8*A +2.0*B/3.0 +0.5*G)*X2 +(A/35. +B/40. +G/50.)*X4;
+            // combining the terms makes things worse at smallest Q in single precision
+            //FF = (8-0.8*X2)*A +(3.0-X2/3.)*2*B + (4+0.5*X2)*G +(A/35. +B/40. +G/50.)*X4;
+            // note that G = -volfraction*A/2, combining this makes no further difference at smallest Q
+            //FF = (8 +2.*volfraction + ( volfraction/4. -0.8 +(volfraction/100. -1./35.)*X2 )*X2 )*A  + (3.0 -X2/3. +X4/40)*2*B;
+            HARDSPH= 1./(1. + volfraction*FF );
+            return(HARDSPH);
+      }
+      SINCOS(X,S,C);
+
+// RKH Feb 2016, use version from FISH code as it is better than original sasview one at small Q in single precision
+      FF=A*(S-X*C)/X + B*(2.*X*S -(X2-2.)*C -2.)/X2 + G*( (4.*X2*X -24.*X)*S -(X4 -12.*X2 +24.)*C +24. )/X4;
+      HARDSPH= 1./(1. + 24.*volfraction*FF/X2 );
+
+// rearrange the terms, is now about same as sasmodels
+//     FF=A*(S/X-C) + B*(2.*S/X - C +2.0*(C-1.0)/X2) + G*( (4./X -24./X3)*S -(1.0 -12./X2 +24./X4)*C +24./X4 );
+//     HARDSPH= 1./(1. + 24.*volfraction*FF/X2 );
+// remove 1/X2 from final line, take more powers of X inside the brackets, stil bad
+//      FF=A*(S/X3-C/X2) + B*(2.*S/X3 - C/X2 +2.0*(C-1.0)/X4) + G*( (4./X -24./X3)*S -(1.0 -12./X2 +24./X4)*C +24./X4 )/X2;
+//      HARDSPH= 1./(1. + 24.*volfraction*FF );
+      return(HARDSPH);
    """
 
@@ -106,8 +121,8 @@
 oldname = 'HardsphereStructure'
 oldpars = dict()
-
+# Q=0.001 is in the Taylor series, low Q part, so add Q=0.1, assuming double precision sasview is correct
 tests = [
         [ {'scale': 1.0, 'background' : 0.0, 'effect_radius' : 50.0, 'volfraction' : 0.2,
-           'effect_radius_pd' : 0}, [0.001], [0.209128]]
+           'effect_radius_pd' : 0}, [0.001,0.1], [0.209128,0.930587]]
         ]
 
Index: sasmodels/models/hollow_cylinder.py
===================================================================
--- sasmodels/models/hollow_cylinder.py	(revision 0420af743cbd4738603a3f63dde3bffd78ebfdcd)
+++ sasmodels/models/hollow_cylinder.py	(revision ec2ca99b9e3d90e4b91919c523a94d9855ca9c3c)
@@ -78,19 +78,26 @@
 """
 category = "shape:cylinder"
-
-#             ["name", "units", default, [lower, upper], "type","description"],
+# pylint: disable=bad-whitespace, line-too-long
+#   ["name", "units", default, [lower, upper], "type","description"],
 parameters = [
-              ["radius", "Ang", 30.0, [0, inf], "volume", "Cylinder radius"],
-              ["core_radius", "Ang", 20.0, [0, inf], "volume", "Hollow core radius"],
-              ["length", "Ang", 400.0, [0, inf], "volume", "Cylinder length"],
-              ["sld", "1/Ang^2", 6.3, [-inf, inf], "", "Cylinder sld"],
-              ["solvent_sld", "1/Ang^2", 1, [-inf, inf], "", "Solvent sld"],
-              ["theta", "degrees", 90, [-360, 360], "orientation", "Theta angle"],
-              ["phi", "degrees", 0, [-360, 360], "orientation", "Phi angle"],
-              ]
+    ["radius",      "Ang",     30.0, [0, inf],    "volume",      "Cylinder radius"],
+    ["core_radius", "Ang",     20.0, [0, inf],    "volume",      "Hollow core radius"],
+    ["length",      "Ang",    400.0, [0, inf],    "volume",      "Cylinder length"],
+    ["sld",         "1/Ang^2",  6.3, [-inf, inf], "",            "Cylinder sld"],
+    ["solvent_sld", "1/Ang^2",  1,   [-inf, inf], "",            "Solvent sld"],
+    ["theta",       "degrees", 90,   [-360, 360], "orientation", "Theta angle"],
+    ["phi",         "degrees",  0,   [-360, 360], "orientation", "Phi angle"],
+    ]
+# pylint: enable=bad-whitespace, line-too-long
 
 source = ["lib/J1.c", "lib/gauss76.c", "hollow_cylinder.c"]
 
 def ER(radius, core_radius, length):
+    """
+    :param radius:      Cylinder radius
+    :param core_radius: Hollow core radius, UNUSED
+    :param length:      Cylinder length
+    :return:            Effective radius
+    """
     if radius == 0 or length == 0:
         return 0.0
@@ -104,4 +111,10 @@
 
 def VR(radius, core_radius, length):
+    """
+    :param radius:      Cylinder radius
+    :param core_radius: Hollow core radius
+    :param length:      Cylinder length
+    :return:            Volf ratio for P(q)*S(q)
+    """
     vol_core = pi*core_radius*core_radius*length
     vol_total = pi*radius*radius*length
@@ -110,23 +123,23 @@
 
 # parameters for demo
-demo = dict(scale=1.0,background=0.0,length=400.0,radius=30.0,core_radius=20.0,
-            sld=6.3,solvent_sld=1,theta=90,phi=0,
+demo = dict(scale=1.0, background=0.0, length=400.0, radius=30.0,
+            core_radius=20.0, sld=6.3, solvent_sld=1, theta=90, phi=0,
             radius_pd=.2, radius_pd_n=9,
             length_pd=.2, length_pd_n=10,
             core_radius_pd=.2, core_radius_pd_n=9,
             theta_pd=10, theta_pd_n=5,
-            )
+           )
 
 # For testing against the old sasview models, include the converted parameter
 # names and the target sasview model name.
 oldname = 'HollowCylinderModel'
-oldpars = dict(scale='scale',background='background',radius='radius',
-               core_radius='core_radius',sld='sldCyl',length='length',
-               solvent_sld='sldSolv',phi='axis_phi',theta='axis_theta')
+oldpars = dict(scale='scale', background='background', radius='radius',
+               core_radius='core_radius', sld='sldCyl', length='length',
+               solvent_sld='sldSolv', phi='axis_phi', theta='axis_theta')
 
 # Parameters for unit tests
 tests = [
-         [{"radius" : 30.0},0.00005,1764.926],
-         [{},'VR',1.8],
-         [{},0.001,1756.76]
-         ]
+    [{"radius": 30.0}, 0.00005, 1764.926],
+    [{}, 'VR', 1.8],
+    [{}, 0.001, 1756.76]
+    ]
Index: sasmodels/models/lamellarCaille.py
===================================================================
--- sasmodels/models/lamellarCaille.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
+++ sasmodels/models/lamellarCaille.py	(revision ec2ca99b9e3d90e4b91919c523a94d9855ca9c3c)
@@ -1,3 +1,2 @@
-# Note: model title and parameter table are inserted automatically
 r"""
 This model provides the scattering intensity, $I(q) = P(q) S(q)$, for a
@@ -77,26 +76,26 @@
 title = "Random lamellar sheet with Caille structure factor"
 description = """\
-	[Random lamellar phase with Caille  structure factor]
-        randomly oriented stacks of infinite sheets
-		with Caille S(Q), having polydisperse spacing.
-    	sld = sheet scattering length density
-		sld_solvent = solvent scattering length density
-		background = incoherent background
-		scale = scale factor
+    [Random lamellar phase with Caille  structure factor]
+    randomly oriented stacks of infinite sheets
+    with Caille S(Q), having polydisperse spacing.
+    sld = sheet scattering length density
+    sld_solvent = solvent scattering length density
+    background = incoherent background
+    scale = scale factor
 """
 category = "shape:lamellae"
 
 single = False
-
+# pylint: disable=bad-whitespace, line-too-long
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["thickness", "Ang",  30.0, [0, inf], "volume", "sheet thickness"],
-              ["Nlayers", "",  20, [0, inf], "", "Number of layers"],
-              ["spacing", "Ang", 400., [0.0,inf], "volume", "d-spacing of Caille S(Q)"],
-              ["Caille_parameter", "1/Ang^2", 0.1, [0.0,0.8], "", "Caille parameter"],
-              ["sld", "1e-6/Ang^2", 6.3, [-inf,inf], "",
-               "layer scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2", 1.0, [-inf,inf], "",
-               "Solvent scattering length density"],
-             ]
+parameters = [
+    ["thickness",        "Ang",      30.0,  [0, inf],   "volume", "sheet thickness"],
+    ["Nlayers",          "",          20,   [0, inf],   "",       "Number of layers"],
+    ["spacing",          "Ang",      400.,  [0.0,inf],  "volume", "d-spacing of Caille S(Q)"],
+    ["Caille_parameter", "1/Ang^2",    0.1, [0.0,0.8],  "",       "Caille parameter"],
+    ["sld",              "1e-6/Ang^2", 6.3, [-inf,inf], "",       "layer scattering length density"],
+    ["solvent_sld",      "1e-6/Ang^2", 1.0, [-inf,inf], "",       "Solvent scattering length density"],
+    ]
+# pylint: enable=bad-whitespace, line-too-long
 
 source = ["lamellarCaille_kernel.c"]
@@ -116,16 +115,18 @@
 
 demo = dict(scale=1, background=0,
-            thickness=67.,Nlayers=3.75,spacing=200.,
-            Caille_parameter=0.268,sld=1.0, solvent_sld=6.34,
-            thickness_pd= 0.1, thickness_pd_n=100,
-            spacing_pd= 0.05, spacing_pd_n=40)
+            thickness=67., Nlayers=3.75, spacing=200.,
+            Caille_parameter=0.268, sld=1.0, solvent_sld=6.34,
+            thickness_pd=0.1, thickness_pd_n=100,
+            spacing_pd=0.05, spacing_pd_n=40)
 
 oldname = 'LamellarPSModel'
-oldpars = dict(thickness='delta', Nlayers='N_plates', Caille_parameter='caille',
-               sld='sld_bi',solvent_sld='sld_sol')
+oldpars = dict(thickness='delta', Nlayers='N_plates',
+               Caille_parameter='caille',
+               sld='sld_bi', solvent_sld='sld_sol')
 #
 tests = [
-        [ {'scale': 1.0, 'background' : 0.0, 'thickness' : 30.,'Nlayers' : 20.0, 'spacing' : 400.,
-            'Caille_parameter' : 0.1, 'sld' : 6.3, 'solvent_sld' : 1.0,
-            'thickness_pd' : 0.0, 'spacing_pd' : 0.0 }, [0.001], [28895.13397]]
-        ]
+    [{'scale': 1.0, 'background': 0.0, 'thickness': 30., 'Nlayers': 20.0,
+      'spacing': 400., 'Caille_parameter': 0.1, 'sld': 6.3,
+      'solvent_sld': 1.0, 'thickness_pd': 0.0, 'spacing_pd': 0.0},
+     [0.001], [28895.13397]]
+    ]
Index: sasmodels/models/lib/nr_bess_j1.c
===================================================================
--- sasmodels/models/lib/nr_bess_j1.c	(revision a8b3cdb60df4f2e22c9b3c5c3110bb51b7dba2c6)
+++ sasmodels/models/lib/nr_bess_j1.c	(revision a8b3cdb60df4f2e22c9b3c5c3110bb51b7dba2c6)
@@ -0,0 +1,28 @@
+double NR_BessJ1(double x);
+double NR_BessJ1(double x)
+{
+    double ax,z;
+    double xx,y,ans,ans1,ans2;
+
+    if ((ax=fabs(x)) < 8.0) {
+        y=x*x;
+        ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
+                                                  +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
+        ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
+                                               +y*(99447.43394+y*(376.9991397+y*1.0))));
+        ans=ans1/ans2;
+    } else {
+        z=8.0/ax;
+        y=z*z;
+        xx=ax-2.356194491;
+        ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
+                                   +y*(0.2457520174e-5+y*(-0.240337019e-6))));
+        ans2=0.04687499995+y*(-0.2002690873e-3
+                              +y*(0.8449199096e-5+y*(-0.88228987e-6
+                                                     +y*0.105787412e-6)));
+        ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
+        if (x < 0.0) ans = -ans;
+    }
+
+    return(ans);
+}
Index: sasmodels/models/lib/sphere_form.c
===================================================================
--- sasmodels/models/lib/sphere_form.c	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
+++ sasmodels/models/lib/sphere_form.c	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
@@ -0,0 +1,26 @@
+double
+sphere_form(double q, double radius, double sld, double solvent_sld)
+{
+    const double qr = q*radius;
+    double sn, cn;
+    SINCOS(qr, sn, cn);
+    const double vol = 1.333333333333333*M_PI*radius*radius*radius;
+    const double bes = (qr == 0.0 ? 1.0 : 3.0*(sn-qr*cn)/(qr*qr*qr));
+    const double fq = bes * (sld - solvent_sld)*vol;
+    return 1.0e-4*fq*fq;
+}
+
+// Do we even need the paracrystal form?
+// The only difference is the return value at r=0
+double
+sphere_form_paracrystal(double q, double radius, double delrho)
+{
+    const double qr = q*radius;
+    const double vol = 1.333333333333333*M_PI*radius*radius*radius;
+	const double f = vol*delrho*delrho*1.0e-4;
+    double sn, cn;
+    SINCOS(qr, sn, cn);
+	const double bes = (qr == 0.0 ? f : 3.0*(sn-qr*cn)/(qr*qr*qr));
+	const double fq = bes*delrho*vol;
+	return fq*fq * 1.0e-4;
+}
Index: sasmodels/models/porod.py
===================================================================
--- sasmodels/models/porod.py	(revision 16bb4337d95ce5f999aec885781f06ccb80925c0)
+++ sasmodels/models/porod.py	(revision 16bb4337d95ce5f999aec885781f06ccb80925c0)
@@ -0,0 +1,58 @@
+r"""
+This model fits the Porod function
+
+.. math::
+
+    I(q) = C/q^4
+    \\
+    C = 2\pi (\Delta\rho)^2 S_v
+
+to the data directly without any need for linearisation (cf. Log I(q) vs Log q).
+
+Here $C$ is the scale factor and $S_v$ is the specific surface area (ie, surface area / volume)
+of the sample, and $\Delta\rho$ is the contrast factor.
+
+For 2D data: The 2D scattering intensity is calculated in the same way as 1D, 
+where the q vector is defined as
+
+.. math::
+    q = \sqrt{q_x^2+q_y^2}
+
+"""
+
+from numpy import sqrt, power
+
+name = "porod"
+title = "Porod function"
+description = """\
+          I(q) = scale/q^4 + background
+"""
+
+category = "shape-independent"
+
+parameters = []
+
+def Iq(q):
+    """
+    @param q: Input q-value
+    """
+    return 1.0/power(q, 4)
+
+Iq.vectorized = True  # Iq accepts an array of q values
+
+def Iqxy(qx, qy, *args):
+    """
+    @param qx:   Input q_x-value
+    @param qy:   Input q_y-value
+    @param args: Remaining arguments
+    """
+    return Iq(sqrt(qx ** 2 + qy ** 2), *args)
+
+Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
+
+demo = dict(scale=1.5, background=0.5)
+
+oldname = "PorodModel"
+oldpars = dict(scale='scale', background='background')
+
+tests = [[{'scale': 0.00001, 'background':0.01}, 0.04, 3.916250]]
Index: sasmodels/models/sc_crystal.c
===================================================================
--- sasmodels/models/sc_crystal.c	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
+++ sasmodels/models/sc_crystal.c	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
@@ -0,0 +1,228 @@
+double form_volume(double radius);
+
+double Iq(double q,
+          double dnn,
+          double d_factor,
+          double radius,
+          double sphere_sld,
+          double solvent_sld);
+
+double Iqxy(double qx, double qy,
+            double dnn,
+            double d_factor,
+            double radius,
+            double sphere_sld,
+            double solvent_sld,
+            double theta,
+            double phi,
+            double psi);
+
+double form_volume(double radius)
+{
+    return 1.333333333333333*M_PI*radius*radius*radius;
+}
+
+static double
+sc_eval(double theta, double phi, double temp3, double temp4, double temp5)
+{
+    double cnt, snt;
+    SINCOS(theta, cnt, snt);
+
+    double cnp, snp;
+    SINCOS(phi, cnp, snp);
+
+	double temp6 = snt;
+	double temp7 = -1.0*temp3*snt*cnp;
+	double temp8 = temp3*snt*snp;
+	double temp9 = temp3*cnt;
+	double result = temp6/((1.0-temp4*cos((temp7))+temp5)*
+	                       (1.0-temp4*cos((temp8))+temp5)*
+	                       (1.0-temp4*cos((temp9))+temp5));
+	return (result);
+}
+
+static double
+sc_integrand(double dnn, double d_factor, double qq, double xx, double yy)
+{
+    //Function to calculate integrand values for simple cubic structure
+
+	double da = d_factor*dnn;
+	double temp1 = qq*qq*da*da;
+	double temp2 = pow( 1.0-exp(-1.0*temp1) ,3);
+	double temp3 = qq*dnn;
+	double temp4 = 2.0*exp(-0.5*temp1);
+	double temp5 = exp(-1.0*temp1);
+
+	double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5);
+	integrand *= 2.0/M_PI;
+
+	return(integrand);
+}
+
+static
+double sc_crystal_kernel(double q,
+          double dnn,
+          double d_factor,
+          double radius,
+          double sphere_sld,
+          double solvent_sld)
+{
+	const double va = 0.0;
+	const double vb = M_PI/2.0; //orientation average, outer integral
+
+    double summ=0.0;
+    double answer=0.0;
+
+	for(int i=0;i<150;i++) {
+		//setup inner integral over the ellipsoidal cross-section
+		double summj=0.0;
+		double zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;
+		for(int j=0;j<150;j++) {
+			//150 gauss points for the inner integral
+			double zij = ( Gauss150Z[j]*(vb-va) + va + vb )/2.0;
+			double tmp = Gauss150Wt[j] * sc_integrand(dnn,d_factor,q,zi,zij);
+			summj += tmp;
+		}
+		//now calculate the value of the inner integral
+		answer = (vb-va)/2.0*summj;
+
+		//now calculate outer integral
+		double tmp = Gauss150Wt[i] * answer;
+		summ += tmp;
+	}		//final scaling is done at the end of the function, after the NT_FP64 case
+
+	answer = (vb-va)/2.0*summ;
+
+	//Volume fraction calculated from lattice symmetry and sphere radius
+	const double latticeScale = (4.0/3.0)*M_PI*(radius*radius*radius)/pow(dnn,3);
+
+	//answer *= sphere_form_paracrystal(q, radius,contrast)*latticeScale;
+	answer *= sphere_form(q, radius, sphere_sld, solvent_sld)*latticeScale;
+
+	return answer;
+}
+
+static
+double sc_crystal_kernel_2d(double q, double q_x, double q_y,
+          double dnn,
+          double d_factor,
+          double radius,
+          double sphere_sld,
+          double solvent_sld,
+          double theta,
+          double phi,
+          double psi)
+{
+    //convert angle degree to radian
+    theta = theta * M_PI_180;
+    phi = phi * M_PI_180;
+    psi = psi * M_PI_180;
+
+    const double qda_2 = pow(q*d_factor*dnn,2.0);
+
+    double snt, cnt;
+    SINCOS(theta, snt, cnt);
+
+    double snp, cnp;
+    SINCOS(phi, snp, cnp);
+
+    double sns, cns;
+    SINCOS(psi, sns, cns);
+
+    /// Angles here are respect to detector coordinate instead of against
+    //  q coordinate(PRB 36, 3, 1754)
+    // a3 axis orientation
+
+    const double a3_x = cnt * cnp;
+    const double a3_y = snt;
+
+    // Compute the angle btw vector q and the a3 axis
+    double cos_val_a3 = a3_x*q_x + a3_y*q_y;
+
+    // a1 axis orientation
+    const double a1_x = -cnp*sns * snt+snp*cns;
+    const double a1_y = sns*cnt;
+
+    double cos_val_a1 = a1_x*q_x + a1_y*q_y;
+
+    // a2 axis orientation
+    const double a2_x = -snt*cns*cnp-sns*snp;
+    const double a2_y = cnt*cns;
+
+    // a2 axis
+    const double cos_val_a2 =  a2_x*q_x + a2_y*q_y;
+
+    // The following test should always pass
+    if (fabs(cos_val_a3)>1.0) {
+        //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+        cos_val_a3 = 1.0;
+    }
+    if (fabs(cos_val_a1)>1.0) {
+        //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+        cos_val_a1 = 1.0;
+    }
+    if (fabs(cos_val_a2)>1.0) {
+        //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+        cos_val_a3 = 1.0;
+    }
+
+    const double a3_dot_q = dnn*q*cos_val_a3;
+    const double a1_dot_q = dnn*q*cos_val_a1;
+    const double a2_dot_q = dnn*q*cos_val_a2;
+
+    // Call Zq=Z1*Z2*Z3
+    double Zq = (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a1_dot_q)+exp(-qda_2));
+    Zq *= (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a2_dot_q)+exp(-qda_2));
+    Zq *= (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a3_dot_q)+exp(-qda_2));
+
+    // Use SphereForm directly from libigor
+    double answer = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq;
+
+    //consider scales
+    const double latticeScale = (4.0/3.0)*M_PI*(radius*radius*radius)/pow(dnn,3.0);
+    answer *= latticeScale;
+
+    return answer;
+}
+
+double Iq(double q,
+          double dnn,
+          double d_factor,
+          double radius,
+          double sphere_sld,
+          double solvent_sld)
+{
+    return sc_crystal_kernel(q,
+              dnn,
+              d_factor,
+              radius,
+              sphere_sld,
+              solvent_sld);
+}
+
+// Iqxy is never called since no orientation or magnetic parameters.
+double Iqxy(double qx, double qy,
+            double dnn,
+            double d_factor,
+            double radius,
+            double sphere_sld,
+            double solvent_sld,
+            double theta,
+            double phi,
+            double psi)
+{
+    double q = sqrt(qx*qx + qy*qy);
+
+
+    return sc_crystal_kernel_2d(q, qx/q, qy/q,
+                  dnn,
+                  d_factor,
+                  radius,
+                  sphere_sld,
+                  solvent_sld,
+                  theta,
+                  phi,
+                  psi);
+
+}
+
Index: sasmodels/models/sc_crystal.py
===================================================================
--- sasmodels/models/sc_crystal.py	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
+++ sasmodels/models/sc_crystal.py	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
@@ -0,0 +1,168 @@
+r"""
+Calculates the scattering from a **simple cubic lattice** with
+paracrystalline distortion. Thermal vibrations are considered to be
+negligible, and the size of the paracrystal is infinitely large.
+Paracrystalline distortion is assumed to be isotropic and characterized
+by a Gaussian distribution.
+
+Definition
+----------
+
+The scattering intensity $I(q)$ is calculated as
+
+.. math::
+
+    I(q) = \frac{scale}{V_p}V_{lattice}P(q)Z(q)
+
+where scale is the volume fraction of spheres, $V_p$ is the volume of
+the primary particle, $V_{lattice}$ is a volume correction for the crystal
+structure, $P(q)$ is the form factor of the sphere (normalized), and
+$Z(q)$ is the paracrystalline structure factor for a simple cubic structure.
+
+Equation (16) of the 1987 reference is used to calculate $Z(q)$, using
+equations (13)-(15) from the 1987 paper for Z1, Z2, and Z3.
+
+The lattice correction (the occupied volume of the lattice) for a simple
+cubic structure of particles of radius *R* and nearest neighbor separation *D* is
+
+.. math::
+
+    V_{lattice}=\frac{4\pi}{3}\frac{R^3}{D^3}
+
+The distortion factor (one standard deviation) of the paracrystal is included
+in the calculation of $Z(q)$
+
+.. math::
+
+    \Delta a = gD
+
+where *g* is a fractional distortion based on the nearest neighbor distance.
+
+The simple cubic lattice is
+
+.. figure:: img/sc_crystal_fig1.jpg
+
+For a crystal, diffraction peaks appear at reduced q-values given by
+
+.. math::
+
+    \frac{qD}{2\pi} = \sqrt{h^2+k^2+l^2}
+
+where for a simple cubic lattice any h, k, l are allowed and none are
+forbidden. Thus the peak positions correspond to (just the first 5)
+
+.. math::
+
+    \begin{align*}
+    q/q_0 \quad & \quad 1
+                & \sqrt{2} \quad
+                & \quad  \sqrt{3} \quad
+                & \sqrt{4} \quad
+                & \quad \sqrt{5}\quad \\
+    Indices \quad & (100)
+                  & \quad (110) \quad
+                  & \quad (111)
+                  & (200) \quad
+                  & \quad (210)
+    \end{align*}
+
+.. note::
+
+    The calculation of *Z(q)* is a double numerical integral that must be
+    carried out with a high density of points to properly capture the sharp
+    peaks of the paracrystalline scattering.
+    So be warned that the calculation is SLOW. Go get some coffee.
+    Fitting of any experimental data must be resolution smeared for any
+    meaningful fit. This makes a triple integral. Very, very slow.
+    Go get lunch!
+
+This example dataset is produced using 200 data points,
+$q_{min} = 0.01A^{-1}, q_{max} = 0.1A^{-1}$ and the above default values.
+
+.. figure:: img/sc_crystal_1d.jpg
+
+    1D plot in the linear scale using the default values (w/200 data point).
+
+The 2D (Anisotropic model) is based on the reference below where *I(q)* is
+approximated for 1d scattering. Thus the scattering pattern for 2D may not
+be accurate. Note that we are not responsible for any incorrectness of the 2D
+model computation.
+
+.. figure:: img/sc_crystal_fig2.jpg
+.. figure:: img/sc_crystal_fig3.jpg
+
+    2D plot using the default values (w/200X200 pixels).
+
+Reference
+---------
+Hideki Matsuoka et. al. *Physical Review B,* 36 (1987) 1754-1765
+(Original Paper)
+
+Hideki Matsuoka et. al. *Physical Review B,* 41 (1990) 3854 -3856
+(Corrections to FCC and BCC lattice structure calculation)
+
+"""
+
+from numpy import inf
+
+name = "sc_crystal"
+title = "Simple cubic lattice with paracrystalline distortion"
+description = """
+        P(q)=(scale/Vp)*V_lattice*P(q)*Z(q)+bkg where scale is the volume
+        fraction of sphere,
+        Vp = volume of the primary particle,
+        V_lattice = volume correction for
+        for the crystal structure,
+        P(q)= form factor of the sphere (normalized),
+        Z(q)= paracrystalline structure factor
+        for a simple cubic structure.
+        [Simple Cubic ParaCrystal Model]
+        Parameters;
+        scale: volume fraction of spheres
+        bkg:background, R: radius of sphere
+        dnn: Nearest neighbor distance
+        d_factor: Paracrystal distortion factor
+        radius: radius of the spheres
+        sldSph: SLD of the sphere
+        sldSolv: SLD of the solvent
+        """
+category = "shape:paracrystal"
+single = False
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["dnn",         "Ang",       220.0,  [0.0, inf],  "",            "Nearest neighbor distance"],
+              ["d_factor",    "",            0.06, [-inf, inf], "",            "Paracrystal distortion factor"],
+              ["radius",      "Ang",        40.0,  [0.0, inf],  "volume",      "Radius of sphere"],
+              ["sphere_sld",  "1e-6/Ang^2",  3.0,  [0.0, inf],  "",            "Sphere scattering length density"],
+              ["solvent_sld", "1e-6/Ang^2",  6.3,  [0.0, inf],  "",            "Solvent scattering length density"],
+              ["theta",       "degrees",     0.0,  [-inf, inf], "orientation", "Orientation of the a1 axis w/respect incoming beam"],
+              ["phi",         "degrees",     0.0,  [-inf, inf], "orientation", "Orientation of the a2 in the plane of the detector"],
+              ["psi",         "degrees",     0.0,  [-inf, inf], "orientation", "Orientation of the a3 in the plane of the detector"],
+             ]
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lib/sphere_form.c", "lib/gauss150.c", "sc_crystal.c"]
+
+demo = dict(scale=1, background=0,
+            dnn=220.0,
+            d_factor=0.06,
+            radius=40.0,
+            sphere_sld=3.0,
+            solvent_sld=6.3,
+            theta=0.0,
+            phi=0.0,
+            psi=0.0)
+
+oldname = 'SCCrystalModel'
+
+oldpars = dict(sphere_sld='sldSph',
+               solvent_sld='sldSolv')
+
+tests = [
+    # Accuracy tests based on content in test/utest_extra_models.py
+    [{}, 0.001, 10.3038],
+    [{}, 0.215268, 0.00714889],
+    [{}, (0.414467), 0.000313289]
+    ]
+
+
