Index: sasmodels/models/HayterMSAsq.py
===================================================================
--- sasmodels/models/HayterMSAsq.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/HayterMSAsq.py	(revision 7f47777b1828197366e17d650a92f794a3433fcc)
@@ -90,3 +90,9 @@
 # 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)
+#
+# 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]]
+        ]
 
Index: sasmodels/models/barbell.py
===================================================================
--- sasmodels/models/barbell.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/barbell.py	(revision dcdf29d606cfaa2f3bb9d578ac394c6f167fa360)
@@ -100,20 +100,24 @@
 title = "Cylinder with spherical end caps"
 description = """
-        Calculates the scattering from a barbell-shaped cylinder. That is a sphereocylinder with spherical end caps that have a radius larger than that of the cylinder and the center of the end cap
-        radius lies outside of the cylinder.
-        Note: As the length of cylinder(bar) -->0,it becomes a dumbbell. And when rad_bar = rad_bell, it is a spherocylinder.
-        It must be that rad_bar <(=) rad_bell.
+    Calculates the scattering from a barbell-shaped cylinder.
+    That is a sphereocylinder with spherical end caps that have a radius larger
+    than that of the cylinder and the center of the end cap radius lies outside
+    of the cylinder.
+    Note: As the length of cylinder(bar) -->0,it becomes a dumbbell. And when
+    rad_bar = rad_bell, it is a spherocylinder.
+    It must be that rad_bar <(=) rad_bell.
 """
 category = "shape:cylinder"
-
+# pylint: disable=bad-whitespace, line-too-long
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["sld", "4e-6/Ang^2", 4, [-inf, inf], "", "Barbell scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "", "Solvent scattering length density"],
-              ["bell_radius", "Ang", 40, [0, inf], "volume", "Spherical bell radius"],
-              ["radius", "Ang", 20, [0, inf], "volume", "Cylindrical bar radius"],
-              ["length", "Ang", 400, [0, inf], "volume", "Cylinder bar length"],
-              ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"],
-              ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"],
+parameters = [["sld",         "4e-6/Ang^2",   4, [-inf, inf], "",            "Barbell scattering length density"],
+              ["solvent_sld", "1e-6/Ang^2",   1, [-inf, inf], "",            "Solvent scattering length density"],
+              ["bell_radius", "Ang",         40, [0, inf],    "volume",      "Spherical bell radius"],
+              ["radius",      "Ang",         20, [0, inf],    "volume",      "Cylindrical bar radius"],
+              ["length",      "Ang",        400, [0, inf],    "volume",      "Cylinder bar length"],
+              ["theta",       "degrees",     60, [-inf, inf], "orientation", "In plane angle"],
+              ["phi",         "degrees",     60, [-inf, inf], "orientation", "Out of plane angle"],
              ]
+# pylint: enable=bad-whitespace, line-too-long
 
 source = ["lib/J1.c", "lib/gauss76.c", "barbell.c"]
Index: sasmodels/models/bcc.py
===================================================================
--- sasmodels/models/bcc.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/bcc.py	(revision dcdf29d606cfaa2f3bb9d578ac394c6f167fa360)
@@ -109,20 +109,23 @@
 title = "Body-centred cubic lattic with paracrystalline distortion"
 description = """
-    Calculates the scattering from a **body-centered 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.
+    Calculates the scattering from a **body-centered 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.
     """
 category = "shape:paracrystal"
-
+# pylint: disable=bad-whitespace, line-too-long
 #             ["name", "units", default, [lower, upper], "type","description" ],
-parameters = [["dnn", "Ang", 220, [-inf, inf], "", "Nearest neighbour distance"],
-              ["d_factor", "", 0.06, [-inf, inf], "", "Paracrystal distortion factor"],
-              ["radius", "Ang", 40, [0, inf], "volume", "Particle radius"],
-              ["sld", "1e-6/Ang^2", 4, [-inf, inf], "", "Particle scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "", "Solvent scattering length density"],
-              ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"],
-              ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"],
-              ["psi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"]
+parameters = [["dnn",         "Ang",       220,    [-inf, inf], "",            "Nearest neighbour distance"],
+              ["d_factor",    "",            0.06, [-inf, inf], "",            "Paracrystal distortion factor"],
+              ["radius",      "Ang",        40,    [0, inf],    "volume",      "Particle radius"],
+              ["sld",         "1e-6/Ang^2",  4,    [-inf, inf], "",            "Particle scattering length density"],
+              ["solvent_sld", "1e-6/Ang^2",  1,    [-inf, inf], "",            "Solvent scattering length density"],
+              ["theta",       "degrees",    60,    [-inf, inf], "orientation", "In plane angle"],
+              ["phi",         "degrees",    60,    [-inf, inf], "orientation", "Out of plane angle"],
+              ["psi",         "degrees",    60,    [-inf, inf], "orientation", "Out of plane angle"]
              ]
+# pylint: enable=bad-whitespace, line-too-long
 
 source = ["lib/J1.c", "lib/gauss150.c", "bcc.c"]
Index: sasmodels/models/broad_peak.py
===================================================================
--- sasmodels/models/broad_peak.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/broad_peak.py	(revision dcdf29d606cfaa2f3bb9d578ac394c6f167fa360)
@@ -60,14 +60,33 @@
 category = "shape-independent"
 
+# pylint: disable=bad-whitespace, line-too-long
 #             ["name", "units", default, [lower, upper], "type", "description"],
-parameters = [["porod_scale", "", 1.0e-05, [-inf, inf], "", "Power law scale factor"],
-              ["porod_exp", "", 3.0, [-inf, inf], "", "Exponent of power law"],
-              ["lorentz_scale", "", 10.0, [-inf, inf], "", "Scale factor for broad Lorentzian peak"],
-              ["lorentz_length", "Ang", 50.0, [-inf, inf], "", "Lorentzian screening length"],
-              ["peak_pos", "1/Ang", 0.1, [-inf, inf], "", "Peak postion in q"],
-              ["lorentz_exp", "", 2.0, [-inf, inf], "", "exponent of Lorentz function"],
+parameters = [["porod_scale",    "",  1.0e-05, [-inf, inf], "", "Power law scale factor"],
+              ["porod_exp",      "",      3.0, [-inf, inf], "", "Exponent of power law"],
+              ["lorentz_scale",  "",     10.0, [-inf, inf], "", "Scale factor for broad Lorentzian peak"],
+              ["lorentz_length", "Ang",  50.0, [-inf, inf], "", "Lorentzian screening length"],
+              ["peak_pos",       "1/Ang", 0.1, [-inf, inf], "", "Peak position in q"],
+              ["lorentz_exp",    "",      2.0, [-inf, inf], "", "Exponent of Lorentz function"],
              ]
+# pylint: enable=bad-whitespace, line-too-long
 
-def Iq(q, porod_scale, porod_exp, lorentz_scale, lorentz_length, peak_pos, lorentz_exp):
+def Iq(q,
+       porod_scale=1.0e-5,
+       porod_exp=3.0,
+       lorentz_scale=10.0,
+       lorentz_length=50.0,
+       peak_pos=0.1,
+       lorentz_exp=2.0):
+    """
+    :param q:              Input q-value
+    :param porod_scale:    Power law scale factor
+    :param porod_exp:      Exponent of power law
+    :param lorentz_scale:  Scale factor for broad Lorentzian peak
+    :param lorentz_length: Lorentzian screening length
+    :param peak_pos:       Peak position in q
+    :param lorentz_exp:    Exponent of Lorentz function
+    :return:               Calculated intensity
+    """
+
     inten = (porod_scale / q ** porod_exp + lorentz_scale
              / (1.0 + (abs(q - peak_pos) * lorentz_length) ** lorentz_exp))
@@ -76,7 +95,13 @@
 
 def Iqxy(qx, qy, *args):
+    """
+    :param qx:   Input q_x-value
+    :param qy:   Input q_y-value
+    :param args: Remaining arguments
+    :return:     2D-Intensity
+    """
     return Iq(sqrt(qx ** 2 + qy ** 2), *args)
+
 Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
-
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/capped_cylinder.py
===================================================================
--- sasmodels/models/capped_cylinder.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/capped_cylinder.py	(revision dcdf29d606cfaa2f3bb9d578ac394c6f167fa360)
@@ -125,11 +125,10 @@
 """
 category = "shape:cylinder"
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type", "description"],
+parameters = [["sld",         "1e-6/Ang^2", 4, [-inf, inf], "",       "Cylinder scattering length density"],
+              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "",       "Solvent scattering length density"],
+              ["radius",      "Ang",       20, [0, inf],    "volume", "Cylinder radius"],
 
-#             ["name", "units", default, [lower, upper], "type", "description"],
-parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "",
-               "Cylinder scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "",
-               "Solvent scattering length density"],
-              ["radius", "Ang", 20, [0, inf], "volume", "Cylinder radius"],
               # TODO: use an expression for cap radius with fixed bounds.
               # The current form requires cap radius R bigger than cylinder radius r.
@@ -141,9 +140,10 @@
               # in the capped cylinder, and zero for no bar in the barbell model.  In
               # both models, one would be a pill.
-              ["cap_radius", "Ang", 20, [0, inf], "volume", "Cap radius"],
-              ["length", "Ang", 400, [0, inf], "volume", "Cylinder length"],
-              ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"],
-              ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"],
+              ["cap_radius", "Ang",     20, [0, inf],    "volume", "Cap radius"],
+              ["length",     "Ang",    400, [0, inf],    "volume", "Cylinder length"],
+              ["theta",      "degrees", 60, [-inf, inf], "orientation", "In plane angle"],
+              ["phi",        "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"],
              ]
+# pylint: enable=bad-whitespace, line-too-long
 
 source = ["lib/J1.c", "lib/gauss76.c", "capped_cylinder.c"]
@@ -159,4 +159,7 @@
             phi_pd=15, phi_pd_n=1)
 oldname = 'CappedCylinderModel'
-oldpars = dict(sld='sld_capcyl', solvent_sld='sld_solv',
-               length='len_cyl', radius='rad_cyl', cap_radius='rad_cap')
+oldpars = dict(sld='sld_capcyl',
+               solvent_sld='sld_solv',
+               length='len_cyl',
+               radius='rad_cyl',
+               cap_radius='rad_cap')
Index: sasmodels/models/core_shell_bicelle.c
===================================================================
--- sasmodels/models/core_shell_bicelle.c	(revision 8007311db20a23588c62551c8a69f089cf8835e1)
+++ sasmodels/models/core_shell_bicelle.c	(revision 8007311db20a23588c62551c8a69f089cf8835e1)
@@ -0,0 +1,209 @@
+double form_volume(double radius, double rim_thickness, double face_thickness, double length);
+double Iq(double q,
+          double radius,
+          double rim_thickness,
+          double face_thickness,
+          double length,
+          double core_sld,
+          double face_sld,
+          double rim_sld,
+          double solvent_sld);
+
+
+double Iqxy(double qx, double qy,
+          double radius,
+          double rim_thickness,
+          double face_thickness,
+          double length,
+          double core_sld,
+          double face_sld,
+          double rim_sld,
+          double solvent_sld,
+          double theta,
+          double phi);
+
+
+double form_volume(double radius, double rim_thickness, double face_thickness, double length)
+{
+    return M_PI*(radius+rim_thickness)*(radius+rim_thickness)*(length+2*face_thickness);
+}
+
+static double
+bicelle_kernel(double qq,
+              double rad,
+              double radthick,
+              double facthick,
+              double length,
+              double rhoc,
+              double rhoh,
+              double rhor,
+              double rhosolv,
+              double dum)
+{
+	double dr1,dr2,dr3;
+	double besarg1,besarg2;
+	double vol1,vol2,vol3;
+	double sinarg1,sinarg2;
+	double t1,t2,t3;
+	double retval,si1,si2,be1,be2;
+
+	dr1 = rhoc-rhoh;
+	dr2 = rhor-rhosolv;
+	dr3=  rhoh-rhor;
+	vol1 = M_PI*rad*rad*(2.0*length);
+	vol2 = M_PI*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick);
+	vol3= M_PI*(rad)*(rad)*(2.0*length+2.0*facthick);
+	besarg1 = qq*rad*sin(dum);
+	besarg2 = qq*(rad+radthick)*sin(dum);
+	sinarg1 = qq*length*cos(dum);
+	sinarg2 = qq*(length+facthick)*cos(dum);
+
+	if(besarg1 == 0) {
+		be1 = 0.5;
+	} else {
+		be1 = J1(besarg1)/besarg1;
+	}
+	if(besarg2 == 0) {
+		be2 = 0.5;
+	} else {
+		be2 = J1(besarg2)/besarg2;
+	}
+	if(sinarg1 == 0) {
+		si1 = 1.0;
+	} else {
+		si1 = sin(sinarg1)/sinarg1;
+	}
+	if(sinarg2 == 0) {
+		si2 = 1.0;
+	} else {
+		si2 = sin(sinarg2)/sinarg2;
+	}
+	t1 = 2.0*vol1*dr1*si1*be1;
+	t2 = 2.0*vol2*dr2*si2*be2;
+	t3 = 2.0*vol3*dr3*si2*be1;
+
+	retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum);
+	return(retval);
+
+}
+
+static double
+bicelle_integration(double qq,
+                   double rad,
+                   double radthick,
+                   double facthick,
+                   double length,
+                   double rhoc,
+                   double rhoh,
+                   double rhor,
+                   double rhosolv)
+{
+
+
+	double answer,halfheight;
+	double lolim,uplim,summ,yyy,zi;
+	int nord,i;
+
+	// set up the integration end points
+	nord = 76;
+	lolim = 0.0;
+	uplim = M_PI/2;
+	halfheight = length/2.0;
+
+	summ = 0.0;
+	i=0;
+	for(i=0;i<nord;i++) {
+		zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
+		yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick,
+		                     halfheight, rhoc, rhoh, rhor,rhosolv, zi);
+		summ += yyy;
+	}
+
+	// calculate value of integral to return
+	answer = (uplim-lolim)/2.0*summ;
+	return(answer);
+}
+
+static double
+bicelle_kernel_2d(double q, double q_x, double q_y,
+          double radius,
+          double rim_thickness,
+          double face_thickness,
+          double length,
+          double core_sld,
+          double face_sld,
+          double rim_sld,
+          double solvent_sld,
+          double theta,
+          double phi)
+{
+    double cyl_x, cyl_y;
+    double alpha, cos_val;
+    double answer;
+
+    //convert angle degree to radian
+    theta *= M_PI/180.0;
+    phi *= M_PI/180.0;
+
+    // Cylinder orientation
+    cyl_x = cos(theta) * cos(phi);
+    cyl_y = sin(theta);
+
+    // Compute the angle btw vector q and the axis of the cylinder
+    cos_val = cyl_x*q_x + cyl_y*q_y;
+    alpha = acos( cos_val );
+
+    // Get the kernel
+    answer = bicelle_kernel(q, radius, rim_thickness, face_thickness,
+                           length/2.0, core_sld, face_sld, rim_sld,
+                           solvent_sld, alpha) / fabs(sin(alpha));
+
+    answer *= 1.0e-4;
+
+    return answer;
+}
+
+double Iq(double q,
+          double radius,
+          double rim_thickness,
+          double face_thickness,
+          double length,
+          double core_sld,
+          double face_sld,
+          double rim_sld,
+          double solvent_sld)
+{
+    double intensity = bicelle_integration(q, radius, rim_thickness, face_thickness,
+                       length, core_sld, face_sld, rim_sld, solvent_sld);
+    return intensity*1.0e-4;
+}
+
+
+double Iqxy(double qx, double qy,
+          double radius,
+          double rim_thickness,
+          double face_thickness,
+          double length,
+          double core_sld,
+          double face_sld,
+          double rim_sld,
+          double solvent_sld,
+          double theta,
+          double phi)
+{
+    double q;
+    q = sqrt(qx*qx+qy*qy);
+    double intensity = bicelle_kernel_2d(q, qx/q, qy/q,
+                      radius,
+                      rim_thickness,
+                      face_thickness,
+                      length,
+                      core_sld,
+                      face_sld,
+                      rim_sld,
+                      solvent_sld,
+                      theta,
+                      phi);
+
+    return intensity;
+}
Index: sasmodels/models/core_shell_bicelle.py
===================================================================
--- sasmodels/models/core_shell_bicelle.py	(revision 8007311db20a23588c62551c8a69f089cf8835e1)
+++ sasmodels/models/core_shell_bicelle.py	(revision 8007311db20a23588c62551c8a69f089cf8835e1)
@@ -0,0 +1,151 @@
+r"""
+
+Definition
+----------
+This model provides the form factor for a circular cylinder with a core-shell
+scattering length density profile.
+The form factor is normalized by the particle volume.
+
+.. _core-shell-cylinder-geometry:
+
+.. figure:: img/core_shell_bicelle_geometry.png
+
+    (Graphic from DOI: 10.1039/C0NP00002G)
+
+
+The output of the 1D scattering intensity function for randomly oriented
+cylinders is then given by the equation above.
+
+The *theta* and *phi* parameters are not used for the 1D output.
+Our implementation of the scattering kernel and the 1D scattering intensity
+use the c-library from NIST.
+
+.. figure:: img/core_shell_bicelle_1d.jpg
+
+    1D plot using the default values (w/200 data point).
+
+.. figure:: img/core_shell_bicelle_fig1.jpg
+
+    Definition of the angles for the oriented CoreShellBicelleModel.
+
+.. figure:: img/core_shell_bicelle_fig2.jpg
+
+    Examples of the angles for oriented pp against the detector plane.
+
+References
+----------
+
+L A Feigin and D I Svergun,
+*Structure Analysis by Small-Angle X-Ray and Neutron Scattering,*
+Plenum Press, New York, (1987)
+
+"""
+
+from numpy import inf, sin, cos
+
+name = "core_shell_bicelle"
+title = "Circular cylinder with a core-shell scattering length density profile.."
+description = """
+    P(q,alpha)= scale/Vs*f(q)^(2) + bkg,  where: f(q)= 2(core_sld
+    - solvant_sld)* Vc*sin[qLcos(alpha/2)]
+    /[qLcos(alpha/2)]*J1(qRsin(alpha))
+    /[qRsin(alpha)]+2(shell_sld-solvent_sld)
+    *Vs*sin[q(L+T)cos(alpha/2)][[q(L+T)
+    *cos(alpha/2)]*J1(q(R+T)sin(alpha))
+    /q(R+T)sin(alpha)]
+
+    alpha:is the angle between the axis of
+    the cylinder and the q-vector
+    Vs: the volume of the outer shell
+    Vc: the volume of the core
+    L: the length of the core
+    shell_sld: the scattering length density
+    of the shell
+    solvent_sld: the scattering length density
+    of the solvent
+    bkg: the background
+    T: the thickness
+    R+T: is the outer radius
+    L+2T: The total length of the outershell
+    J1: the first order Bessel function
+    theta: axis_theta of the cylinder
+    phi: the axis_phi of the cylinder...
+        """
+category = "shape:cylinder"
+
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type", "description"],
+parameters = [
+    ["radius",         "Ang",       20, [0, inf],    "volume",      "Cylinder core radius"],
+    ["rim_thickness",  "Ang",       10, [0, inf],    "volume",      "Rim shell thickness"],
+    ["face_thickness", "Ang",       10, [0, inf],    "volume",      "Cylinder face thickness"],
+    ["length",         "Ang",      400, [0, inf],    "volume",      "Cylinder length"],
+    ["core_sld",       "1e-6/Ang^2", 1, [-inf, inf], "",            "Cylinder core scattering length density"],
+    ["face_sld",       "1e-6/Ang^2", 4, [-inf, inf], "",            "Cylinder face scattering length density"],
+    ["rim_sld",        "1e-6/Ang^2", 4, [-inf, inf], "",            "Cylinder rim scattering length density"],
+    ["solvent_sld",    "1e-6/Ang^2", 1, [-inf, inf], "",            "Solvent scattering length density"],
+    ["theta",          "degrees",   90, [-inf, inf], "orientation", "In plane angle"],
+    ["phi",            "degrees",    0, [-inf, inf], "orientation", "Out of plane angle"],
+    ]
+
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lib/J1.c", "lib/gauss76.c", "core_shell_bicelle.c"]
+
+demo = dict(scale=1, background=0,
+            radius=20.0,
+            rim_thickness=10.0,
+            face_thickness=10.0,
+            length=400.0,
+            core_sld=1.0,
+            face_sld=4.0,
+            rim_sld=4.0,
+            solvent_sld=1.0,
+            theta=90,
+            phi=0)
+
+oldname = 'CoreShellBicelleModel'
+
+oldpars = dict(rim_thickness='rim_thick',
+               face_thickness='face_thick',
+               theta='axis_theta',
+               phi='axis_phi')
+
+qx, qy = 0.4 * cos(90), 0.5 * sin(0)
+tests = [
+    # Accuracy tests based on content in test/utest_other_models.py
+    [{'radius': 20.0,
+      'rim_thickness': 10.0,
+      'face_thickness': 10.0,
+      'length': 400.0,
+      'core_sld': 1.0,
+      'face_sld': 4.0,
+      'rim_sld': 4.0,
+      'solvent_sld': 1.0,
+      'background': 0.0,
+     }, 0.001, 353.550],
+
+    [{'radius': 20.0,
+      'rim_thickness': 10.0,
+      'face_thickness': 10.0,
+      'length': 400.0,
+      'core_sld': 1.0,
+      'face_sld': 4.0,
+      'rim_sld': 4.0,
+      'solvent_sld': 1.0,
+      'theta': 90.0,
+      'phi': 0.0,
+      'background': 0.00,
+     }, (qx, qy), 24.9167],
+
+    # Additional tests with larger range of parameters
+    [{'radius': 3.0,
+      'rim_thickness': 100.0,
+      'face_thickness': 100.0,
+      'length': 1200.0,
+      'core_sld': 5.0,
+      'face_sld': 41.0,
+      'rim_sld': 42.0,
+      'solvent_sld': 21.0,
+     }, 0.05, 1670.1828],
+    ]
Index: sasmodels/models/hardsphere.py
===================================================================
--- sasmodels/models/hardsphere.py	(revision 7ed702f2d0d4ecbcb0115970cdabf0d05ac08de0)
+++ sasmodels/models/hardsphere.py	(revision 7f47777b1828197366e17d650a92f794a3433fcc)
@@ -106,2 +106,7 @@
 oldpars = dict()
 
+tests = [
+        [ {'scale': 1.0, 'background' : 0.0, 'effect_radius' : 50.0, 'volfraction' : 0.2,
+           'effect_radius_pd' : 0}, [0.001], [0.209128]]
+        ]
+
Index: sasmodels/models/lamellarCaille.py
===================================================================
--- sasmodels/models/lamellarCaille.py	(revision d18f8a8e10d65f0ace6fdc9360abd60dfcd6dcce)
+++ sasmodels/models/lamellarCaille.py	(revision 7f47777b1828197366e17d650a92f794a3433fcc)
@@ -122,2 +122,8 @@
 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]]
+        ]
Index: sasmodels/models/lamellarCailleHG.py
===================================================================
--- sasmodels/models/lamellarCailleHG.py	(revision d4666ca59132213b621b9cd27fda96521894152e)
+++ sasmodels/models/lamellarCailleHG.py	(revision 7f47777b1828197366e17d650a92f794a3433fcc)
@@ -139,6 +139,13 @@
 
 oldname = 'LamellarPSHGModel'
+
 oldpars = dict(
     tail_length='deltaT', head_length='deltaH', Nlayers='n_plates',
     Caille_parameter='caille', sld='sld_tail', head_sld='sld_head',
     solvent_sld='sld_solvent')
+#
+tests = [
+        [ {'scale': 1.0, 'background' : 0.0, 'tail_length' : 10.0, 'head_length' : 2.0,'Nlayers' : 30.0, 'spacing' : 40.,
+            'Caille_parameter' : 0.001, 'sld' : 0.4, 'head_sld' : 2.0, 'solvent_sld' : 6.0,
+            'tail_length_pd' : 0.0, 'head_length_pd' : 0.0, 'spacing_pd' : 0.0 }, [0.001], [6838238.571488]]
+        ]
Index: sasmodels/models/lamellarFFHG.py
===================================================================
--- sasmodels/models/lamellarFFHG.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/lamellarFFHG.py	(revision 7f47777b1828197366e17d650a92f794a3433fcc)
@@ -121,3 +121,9 @@
 oldpars = dict(head_length='h_thickness', tail_length='t_length',
                sld='sld_tail', head_sld='sld_head', solvent_sld='sld_solvent')
+#
+tests = [
+        [ {'scale': 1.0, 'background' : 0.0, 'tail_length' : 15.0, 'head_length' : 10.0,'sld' : 0.4,
+         'head_sld' : 3.0, 'solvent_sld' : 6.0, }, [0.001], [653143.9209]]
+        ]
 
+
Index: sasmodels/models/lamellarPC.py
===================================================================
--- sasmodels/models/lamellarPC.py	(revision d18f8a8e10d65f0ace6fdc9360abd60dfcd6dcce)
+++ sasmodels/models/lamellarPC.py	(revision 7f47777b1828197366e17d650a92f794a3433fcc)
@@ -148,2 +148,8 @@
 oldpars = dict(spacing_polydisp='pd_spacing', sld='sld_layer',
                solvent_sld='sld_solvent')
+#
+tests = [
+        [ {'scale': 1.0, 'background' : 0.0, 'thickness' : 33.,'Nlayers' : 20.0, 'spacing' : 250., 'spacing_polydisp' : 0.2,
+            'sld' : 1.0, 'solvent_sld' : 6.34,
+            'thickness_pd' : 0.0, 'thickness_pd_n' : 40 }, [0.001, 0.215268], [21829.3, 0.00487686]]
+        ]
Index: sasmodels/models/linear_pearls.c
===================================================================
--- sasmodels/models/linear_pearls.c	(revision 577912b439fe119a0663629c075b6b67323269ca)
+++ sasmodels/models/linear_pearls.c	(revision 577912b439fe119a0663629c075b6b67323269ca)
@@ -0,0 +1,123 @@
+double form_volume(double radius, double num_pearls);
+
+double Iq(double q,
+            double radius,
+            double edge_sep,
+            double num_pearls,
+            double pearl_sld,
+            double solvent_sld);
+
+double Iqxy(double qx, double qy,
+            double radius,
+            double edge_sep,
+            double num_pearls,
+            double pearl_sld,
+            double solvent_sld);
+
+double linear_pearls_kernel(double q,
+            double radius,
+            double edge_sep,
+            double num_pearls,
+            double pearl_sld,
+            double solvent_sld);
+
+
+double form_volume(double radius, double num_pearls)
+{
+    // Pearl volume
+    double pearl_vol = 4.0 /3.0 * M_PI * pow(radius, 3.0);
+    // Return total volume
+    return num_pearls * pearl_vol;;
+}
+
+// If used elsewhere - factor out to lib/
+static
+double sinc(double x)
+{
+  if (x==0.0){
+    return 1.0;
+  }
+  return sin(x)/x;
+}
+
+double linear_pearls_kernel(double q,
+            double radius,
+            double edge_sep,
+            double num_pearls,
+            double pearl_sld,
+            double solvent_sld)
+{
+    double n_contrib;
+    //relative sld
+    double contrast_pearl = pearl_sld - solvent_sld;
+    //each volume
+    double pearl_vol = 4.0 /3.0 * M_PI * pow(radius, 3.0);
+    //total volume
+    double tot_vol = num_pearls * pearl_vol;
+    //mass
+    double m_s = contrast_pearl * pearl_vol;
+    //center to center distance between the neighboring pearls
+    double separation = edge_sep + 2.0 * radius;
+
+    double x=q*radius;
+
+    // Try Taylor on x*xos(x)
+	// double out_cos = x - pow(x,3)/2 + pow(x,5)/24 - pow(x,7)/720 + pow(x,9)/40320;
+    // psi -= x*out_cos;
+
+    //sine functions of a pearl
+    double psi = sin(q * radius);
+    psi -= x * cos(x);
+    psi /= pow((q * radius), 3.0);
+
+    // N pearls contribution
+    int n_max = num_pearls - 1;
+    n_contrib = num_pearls;
+    for(int num=1; num<=n_max; num++) {
+        n_contrib += (2.0*(num_pearls-num)*sinc(q*separation*num));
+    }
+    // form factor for num_pearls
+    double form_factor = n_contrib;
+    form_factor *= pow((m_s*psi*3.0), 2.0);
+    form_factor /= (tot_vol * 1.0e4);
+
+    return form_factor;
+}
+
+double Iq(double q,
+            double radius,
+            double edge_sep,
+            double num_pearls,
+            double pearl_sld,
+            double solvent_sld)
+{
+
+	double result = linear_pearls_kernel(q,
+                    radius,
+                    edge_sep,
+                    num_pearls,
+                    pearl_sld,
+                    solvent_sld);
+
+	return result;
+}
+
+double Iqxy(double qx, double qy,
+            double radius,
+            double edge_sep,
+            double num_pearls,
+            double pearl_sld,
+            double solvent_sld)
+{
+	double q;
+	q = sqrt(qx*qx+qy*qy);
+
+	double result = linear_pearls_kernel(q,
+                    radius,
+                    edge_sep,
+                    num_pearls,
+                    pearl_sld,
+                    solvent_sld);
+
+	return result;
+}
Index: sasmodels/models/linear_pearls.py
===================================================================
--- sasmodels/models/linear_pearls.py	(revision 577912b439fe119a0663629c075b6b67323269ca)
+++ sasmodels/models/linear_pearls.py	(revision 577912b439fe119a0663629c075b6b67323269ca)
@@ -0,0 +1,110 @@
+r"""
+This model provides the form factor for $N$ spherical pearls of radius $R$
+linearly joined by short strings (or segment length or edge separation)
+$l$ $(= A - 2R)$. $A$ is the center-to-center pearl separation distance.
+The thickness of each string is assumed to be negligible.
+
+.. figure:: img/linear_pearls_fig1.jpg
+
+
+Definition
+----------
+
+The output of the scattering intensity function for the linear_pearls model
+is given by (Dobrynin, 1996)
+
+.. math::
+
+    P(Q) = \frac{scale}{V}\left[ m_{p}^2
+    \left(N+2\sum_{n-1}^{N-1}(N-n)\frac{sin(qnl)}{qnl}\right)
+    \left( 3\frac{sin(qR)-qRcos(qR)}{(qr)^3}\right)^2\right]
+
+where the mass $m_p$ is $(SLD_{pearl}-SLD_{solvent})*(volume\ of\ N\ pearls)$.
+V is the total volume.
+
+The 2D scattering intensity is the same as P(q) above,
+regardless of the orientation of the q vector.
+
+.. figure:: img/linear_pearls_1d.jpg
+
+    1D plot using the default values (w/500 data point).
+
+References
+----------
+
+A V Dobrynin, M Rubinstein and S P Obukhov, *Macromol.*,
+29 (1996) 2974-2979
+
+"""
+
+from numpy import inf
+
+name = "linear_pearls"
+title = "Linear pearls model of scattering from spherical pearls."
+description = """
+    Calculate form factor for Pearl Necklace Model
+    [Macromol. 1996, 29, 2974-2979]
+    Parameters:
+
+    sld_pearl: the SLD of the pearl spheres
+    sld_solv: the SLD of the solvent
+    num_pearls: number of the pearls
+    radius: the radius of a pearl
+    edge_separation: the length of string segment; surface to surface
+    """
+category = "shape:sphere"
+
+# pylint: disable=bad-whitespace, line-too-long
+#            ["name", "units", default, [lower, upper], "type", "description"],
+parameters = [
+    ["radius",      "Ang",       80.0, [0, inf],     "", "Radius of the pearls"],
+    ["edge_sep",    "Ang",      350.0, [0, inf],     "", "Length of the string segment - surface to surface"],
+    ["num_pearls",  "",           3.0, [0, inf],     "", "Number of the pearls"],
+    ["pearl_sld",   "1e-6/Ang^2", 1.0, [-inf, inf],  "", "SLD of the pearl spheres"],
+    ["solvent_sld", "1e-6/Ang^2", 6.3, [-inf, inf],  "", "SLD of the solvent"],
+    ]
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["linear_pearls.c"]
+
+demo = dict(scale=1.0, background=0.0,
+            radius=80.0,
+            edge_sep=350.0,
+            num_pearls=3,
+            pearl_sld=1.0,
+            solvent_sld=6.3)
+
+oldname = "LinearPearlsModel"
+
+oldpars = dict(edge_sep='edge_separation',
+               pearl_sld='sld_pearl',
+               solvent_sld='sld_solv')
+
+"""
+Tests temporarily disabled, until single-double precision accuracy issue solved.
+
+tests = [
+    # Accuracy tests based on content in test/utest_model_pearlnecklace.py
+    [{'radius':      20.0,
+      'num_pearls':   2.0,
+      'pearl_sld':    1.0,
+      'solvent_sld':  6.3,
+      'edge_sep':   400.0,
+     }, 0.001, 185.135],
+
+    # Additional tests with larger range of parameters
+    [{'radius':     120.0,
+      'num_pearls':   5.0,
+      'pearl_sld':    2.0,
+      'solvent_sld':  2.3,
+      'edge_sep':   100.0,
+     }, 0.01, 45.4984],
+
+    [{'radius':       7.0,
+      'num_pearls':   2.0,
+      'pearl_sld':   10.0,
+      'solvent_sld': 99.3,
+      'edge_sep':    20.0,
+     }, 1.0, 0.632811],
+    ]
+"""
Index: sasmodels/models/stickyhardsphere.py
===================================================================
--- sasmodels/models/stickyhardsphere.py	(revision d18f8a8e10d65f0ace6fdc9360abd60dfcd6dcce)
+++ sasmodels/models/stickyhardsphere.py	(revision 7f47777b1828197366e17d650a92f794a3433fcc)
@@ -182,3 +182,9 @@
 demo = dict(effect_radius=200, volfraction=0.2, perturb=0.05,
             stickiness=0.2, effect_radius_pd=0.1, effect_radius_pd_n=40)
+#
+tests = [
+        [ {'scale': 1.0, 'background' : 0.0, 'effect_radius' : 50.0, 'perturb' : 0.05, 'stickiness' : 0.2, 'volfraction' : 0.1,
+           'effect_radius_pd' : 0}, [0.001], [1.09718]]
+        ]
 
+
