Index: explore/realspace.py
===================================================================
--- explore/realspace.py	(revision 362a66f0d054d3213d4e56147f9378d15b5616f0)
+++ explore/realspace.py	(revision 57782799b2cffcbf9e50fd15e7e764c516ec7e9b)
@@ -9,5 +9,5 @@
 import numpy as np
 from numpy import pi, radians, sin, cos, sqrt
-from numpy.random import poisson, uniform, randn, rand
+from numpy.random import poisson, uniform, randn, rand, randint
 from numpy.polynomial.legendre import leggauss
 from scipy.integrate import simps
@@ -78,8 +78,15 @@
 
 
+I3 = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]])
+
 class Shape:
-    rotation = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]])
+    rotation = I3
     center = np.array([0., 0., 0.])[:, None]
     r_max = None
+    lattice_size = np.array((1, 1, 1))
+    lattice_spacing = np.array((1., 1., 1.))
+    lattice_distortion = 0.0
+    lattice_rotation = 0.0
+    lattice_type = ""
 
     def volume(self):
@@ -96,5 +103,6 @@
 
     def rotate(self, theta, phi, psi):
-        self.rotation = rotation(theta, phi, psi) * self.rotation
+        if theta != 0. or phi != 0. or psi != 0.:
+            self.rotation = rotation(theta, phi, psi) * self.rotation
         return self
 
@@ -103,14 +111,79 @@
         return self
 
+    def lattice(self, size=(1, 1, 1), spacing=(2, 2, 2), type="sc",
+                distortion=0.0, rotation=0.0):
+        self.lattice_size = np.asarray(size, 'i')
+        self.lattice_spacing = np.asarray(spacing, 'd')
+        self.lattice_type = type
+        self.lattice_distortion = distortion
+        self.lattice_rotation = rotation
+
     def _adjust(self, points):
-        points = np.asarray(self.rotation * np.matrix(points.T)) + self.center
+        if self.rotation is I3:
+            points = points.T + self.center
+        else:
+            points = np.asarray(self.rotation * np.matrix(points.T)) + self.center
+        if self.lattice_type:
+            points = self._apply_lattice(points)
         return points.T
 
-    def r_bins(self, q, over_sampling=1, r_step=0.):
-        r_max = min(2 * pi / q[0], self.r_max)
+    def r_bins(self, q, over_sampling=10, r_step=0.):
+        if self.lattice_type:
+            r_max = np.sqrt(np.sum(self.lattice_size*self.lattice_spacing*self.dims)**2)/2
+        else:
+            r_max = self.r_max
+        #r_max = min(2 * pi / q[0], r_max)
         if r_step == 0.:
             r_step = 2 * pi / q[-1] / over_sampling
         #r_step = 0.01
         return np.arange(r_step, r_max, r_step)
+
+    def _apply_lattice(self, points):
+        """Spread points to different lattice positions"""
+        size = self.lattice_size
+        spacing = self.lattice_spacing
+        shuffle = self.lattice_distortion
+        rotate = self.lattice_rotation
+        lattice = self.lattice_type
+
+        if rotate != 0:
+            # To vectorize the rotations we will need to unwrap the matrix multiply
+            raise NotImplementedError("don't handle rotations yet")
+
+        # Determine the number of lattice points in the lattice
+        shapes_per_cell = 2 if lattice == "bcc" else 4 if lattice == "fcc" else 1
+        number_of_lattice_points = np.prod(size) * shapes_per_cell
+
+        # For each point in the original shape, figure out which lattice point
+        # to translate it to.  This is both cell index (i*ny*nz + j*nz  + k) as
+        # well as the point in the cell (corner, body center or face center).
+        nsamples = points.shape[1]
+        lattice_point = randint(number_of_lattice_points, size=nsamples)
+
+        # Translate the cell index into the i,j,k coordinates of the senter
+        cell_index = lattice_point // shapes_per_cell
+        center = np.vstack((cell_index//(size[1]*size[2]),
+                            (cell_index%(size[1]*size[2]))//size[2],
+                            cell_index%size[2]))
+        center = np.asarray(center, dtype='d')
+        if lattice == "bcc":
+            center[:, lattice_point % shapes_per_cell == 1] += [[0.5], [0.5], [0.5]]
+        elif lattice == "fcc":
+            center[:, lattice_point % shapes_per_cell == 1] += [[0.0], [0.5], [0.5]]
+            center[:, lattice_point % shapes_per_cell == 2] += [[0.5], [0.0], [0.5]]
+            center[:, lattice_point % shapes_per_cell == 3] += [[0.5], [0.5], [0.0]]
+
+        # Each lattice point has its own displacement from the ideal position.
+        # Not checking that shapes do not overlap if displacement is too large.
+        offset = shuffle*(randn(3, number_of_lattice_points) if shuffle < 0.3
+                          else rand(3, number_of_lattice_points))
+        center += offset[:, cell_index]
+
+        # Each lattice point has its own rotation.  Rotate the point prior to
+        # applying any displacement.
+        # rotation = rotate*(randn(size=(shapes, 3)) if shuffle < 30 else rand(size=(nsamples, 3)))
+        # for k in shapes: points[k] = rotation[k]*points[k]
+        points += center*(np.array([spacing])*np.array(self.dims)).T
+        return points
 
 class Composite(Shape):
@@ -669,6 +742,6 @@
     Iq = 100 * np.ones_like(qx)
     data = Data2D(x=qx, y=qy, z=Iq, dx=None, dy=None, dz=np.sqrt(Iq))
-    data.x_bins = qx[0,:]
-    data.y_bins = qy[:,0]
+    data.x_bins = qx[0, :]
+    data.y_bins = qy[:, 0]
     data.filename = "fake data"
 
@@ -695,5 +768,7 @@
     return shape, fn, fn_xy
 
-def build_sphere(radius=125, rho=2):
+DEFAULT_SPHERE_RADIUS = 125
+DEFAULT_SPHERE_CONTRAST = 2
+def build_sphere(radius=DEFAULT_SPHERE_RADIUS, rho=DEFAULT_SPHERE_CONTRAST):
     shape = TriaxialEllipsoid(radius, radius, radius, rho)
     fn = lambda q: sphere_Iq(q, radius)*rho**2
@@ -751,8 +826,8 @@
     return shape, fn, fn_xy
 
-def build_cubic_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2,
-                  shuffle=0, rotate=0):
+def build_sc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2,
+                        shuffle=0, rotate=0):
     a, b, c = shape.dims
-    shapes = [copy(shape)
+    corners= [copy(shape)
               .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a,
                      (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b,
@@ -762,7 +837,66 @@
               for iy in range(ny)
               for iz in range(nz)]
-    lattice = Composite(shapes)
+    lattice = Composite(corners)
     return lattice
 
+def build_bcc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2,
+                      shuffle=0, rotate=0):
+    a, b, c = shape.dims
+    corners = [copy(shape)
+               .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a,
+                      (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b,
+                      (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c)
+               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate))
+               for ix in range(nx)
+               for iy in range(ny)
+               for iz in range(nz)]
+    centers = [copy(shape)
+               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a,
+                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b,
+                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c)
+               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate))
+               for ix in range(nx)
+               for iy in range(ny)
+               for iz in range(nz)]
+    lattice = Composite(corners + centers)
+    return lattice
+
+def build_fcc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2,
+                      shuffle=0, rotate=0):
+    a, b, c = shape.dims
+    corners = [copy(shape)
+               .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a,
+                      (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b,
+                      (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c)
+               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate))
+               for ix in range(nx)
+               for iy in range(ny)
+               for iz in range(nz)]
+    faces_a = [copy(shape)
+               .shift((ix+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a,
+                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b,
+                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c)
+               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate))
+               for ix in range(nx)
+               for iy in range(ny)
+               for iz in range(nz)]
+    faces_b = [copy(shape)
+               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a,
+                      (iy+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b,
+                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c)
+               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate))
+               for ix in range(nx)
+               for iy in range(ny)
+               for iz in range(nz)]
+    faces_c = [copy(shape)
+               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a,
+                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b,
+                      (iz+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c)
+               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate))
+               for ix in range(nx)
+               for iy in range(ny)
+               for iz in range(nz)]
+    lattice = Composite(corners + faces_a + faces_b + faces_c)
+    return lattice
 
 SHAPE_FUNCTIONS = OrderedDict([
@@ -775,4 +909,10 @@
 ])
 SHAPES = list(SHAPE_FUNCTIONS.keys())
+LATTICE_FUNCTIONS = OrderedDict([
+    ("sc", build_sc_lattice),
+    ("bcc", build_bcc_lattice),
+    ("fcc", build_fcc_lattice),
+])
+LATTICE_TYPES = list(LATTICE_FUNCTIONS.keys())
 
 def check_shape(title, shape, fn=None, show_points=False,
@@ -783,5 +923,7 @@
     r = shape.r_bins(q, r_step=r_step)
     sampling_density = samples / shape.volume
+    print("sampling points")
     rho, points = shape.sample(sampling_density)
+    print("calculating Pr")
     t0 = time.time()
     Pr = calc_Pr(r, rho-rho_solvent, points)
@@ -792,5 +934,5 @@
     import pylab
     if show_points:
-         plot_points(rho, points); pylab.figure()
+        plot_points(rho, points); pylab.figure()
     plot_calc(r, Pr, q, Iq, theory=theory, title=title)
     pylab.gcf().canvas.set_window_title(title)
@@ -806,4 +948,5 @@
     Qx, Qy = np.meshgrid(qx, qy)
     sampling_density = samples / shape.volume
+    print("sampling points")
     t0 = time.time()
     rho, points = shape.sample(sampling_density)
@@ -844,5 +987,8 @@
                         help='lattice size')
     parser.add_argument('-z', '--spacing', type=str, default='2,2,2',
-                        help='lattice spacing')
+                        help='lattice spacing (relative to shape)')
+    parser.add_argument('-t', '--type', choices=LATTICE_TYPES,
+                        default=LATTICE_TYPES[0],
+                        help='lattice type')
     parser.add_argument('-r', '--rotate', type=float, default=0.,
                         help="rotation relative to lattice, gaussian < 30 degrees, uniform otherwise")
@@ -858,8 +1004,36 @@
     nx, ny, nz = [int(v) for v in opts.lattice.split(',')]
     dx, dy, dz = [float(v) for v in opts.spacing.split(',')]
-    shuffle, rotate = opts.shuffle, opts.rotate
+    distortion, rotation = opts.shuffle, opts.rotate
     shape, fn, fn_xy = SHAPE_FUNCTIONS[opts.shape](**pars)
-    if nx > 1 or ny > 1 or nz > 1:
-        shape = build_cubic_lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate)
+    view = tuple(float(v) for v in opts.view.split(','))
+    # If comparing a sphere in a cubic lattice, compare against the
+    # corresponding paracrystalline model.
+    if opts.shape == "sphere" and dx == dy == dz and nx*ny*nz > 1:
+        radius = pars.get('radius', DEFAULT_SPHERE_RADIUS)
+        model_name = opts.type + "_paracrystal"
+        model_pars = {
+            "scale": 1.,
+            "background": 0.,
+            "lattice_spacing": 2*radius*dx,
+            "lattice_distortion": distortion,
+            "radius": radius,
+            "sld": pars.get('rho', DEFAULT_SPHERE_CONTRAST),
+            "sld_solvent": 0.,
+            "theta": view[0],
+            "phi": view[1],
+            "psi": view[2],
+        }
+        fn, fn_xy = wrap_sasmodel(model_name, **model_pars)
+    if nx*ny*nz > 1:
+        if rotation != 0:
+            print("building %s lattice"%opts.type)
+            build_lattice = LATTICE_FUNCTIONS[opts.type]
+            shape = build_lattice(shape, nx, ny, nz, dx, dy, dz,
+                                  distortion, rotation)
+        else:
+            shape.lattice(size=(nx, ny, nz), spacing=(dx, dy, dz),
+                          type=opts.type,
+                          rotation=rotation, distortion=distortion)
+
     title = "%s(%s)" % (opts.shape, " ".join(opts.pars))
     if opts.dim == 1:
@@ -867,5 +1041,4 @@
                     mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples)
     else:
-        view = tuple(float(v) for v in opts.view.split(','))
         check_shape_2d(title, shape, fn_xy, view=view, show_points=opts.plot,
                        mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples)
Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision 610ef23645ef28d7c1e8c87aab1e03c8cd123a93)
+++ sasmodels/convert.py	(revision 610ef23645ef28d7c1e8c87aab1e03c8cd123a93)
@@ -166,5 +166,12 @@
         oldpars = _hand_convert_3_1_2_to_4_1(name, oldpars)
     if version < (4, 2, 0):
+        oldpars = _hand_convert_4_1_to_4_2(name, oldpars)
         oldpars = _rename_magnetic_pars(oldpars)
+    return oldpars
+
+def _hand_convert_4_1_to_4_2(name, oldpars):
+    if name in ('bcc_paracrystal', 'fcc_paracrystal', 'sc_paracrystal'):
+        oldpars['lattice_spacing'] = oldpars.pop('dnn')
+        oldpars['lattice_distortion'] = oldpars.pop('d_factor')
     return oldpars
 
Index: sasmodels/models/bcc_paracrystal.c
===================================================================
--- sasmodels/models/bcc_paracrystal.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
+++ sasmodels/models/bcc_paracrystal.c	(revision 642046e89e9bce7059d153e0f8506ea74d7c1420)
@@ -1,4 +1,4 @@
 static double
-bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
+bcc_Zq(double qa, double qb, double qc, double lattice_spacing, double lattice_distortion)
 {
     // Equations from Matsuoka 26-27-28, multiplied by |q|
@@ -17,10 +17,10 @@
     //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1)
     //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1
-    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    const double arg = -0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3);
     const double exp_arg = exp(arg);
     const double Zq = -cube(expm1(2.0*arg))
-        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0)
-          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0)
-          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0));
+        / ( ((exp_arg - 2.0*cos(lattice_spacing*a1))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(lattice_spacing*a2))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(lattice_spacing*a3))*exp_arg + 1.0));
 
 #elif 0
@@ -36,17 +36,17 @@
     //            = tanh(-a) / [1 - cos(d a_k)/cosh(-a)]
     //
-    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    const double arg = 0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3);
     const double sinh_qd = sinh(arg);
     const double cosh_qd = cosh(arg);
-    const double Zq = sinh_qd/(cosh_qd - cos(dnn*a1))
-                    * sinh_qd/(cosh_qd - cos(dnn*a2))
-                    * sinh_qd/(cosh_qd - cos(dnn*a3));
+    const double Zq = sinh_qd/(cosh_qd - cos(lattice_spacing*a1))
+                    * sinh_qd/(cosh_qd - cos(lattice_spacing*a2))
+                    * sinh_qd/(cosh_qd - cos(lattice_spacing*a3));
 #else
-    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    const double arg = 0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3);
     const double tanh_qd = tanh(arg);
     const double cosh_qd = cosh(arg);
-    const double Zq = tanh_qd/(1.0 - cos(dnn*a1)/cosh_qd)
-                    * tanh_qd/(1.0 - cos(dnn*a2)/cosh_qd)
-                    * tanh_qd/(1.0 - cos(dnn*a3)/cosh_qd);
+    const double Zq = tanh_qd/(1.0 - cos(lattice_spacing*a1)/cosh_qd)
+                    * tanh_qd/(1.0 - cos(lattice_spacing*a2)/cosh_qd)
+                    * tanh_qd/(1.0 - cos(lattice_spacing*a3)/cosh_qd);
 #endif
 
@@ -57,7 +57,7 @@
 // occupied volume fraction calculated from lattice symmetry and sphere radius
 static double
-bcc_volume_fraction(double radius, double dnn)
+bcc_volume_fraction(double radius, double lattice_spacing)
 {
-    return 2.0*sphere_volume(sqrt(0.75)*radius/dnn);
+    return 2.0*sphere_volume(radius/lattice_spacing);
 }
 
@@ -69,6 +69,6 @@
 
 
-static double Iq(double q, double dnn,
-    double d_factor, double radius,
+static double Iq(double q, double lattice_spacing,
+    double lattice_distortion, double radius,
     double sld, double solvent_sld)
 {
@@ -94,5 +94,5 @@
             const double qa = qab*cos_phi;
             const double qb = qab*sin_phi;
-            const double form = bcc_Zq(qa, qb, qc, dnn, d_factor);
+            const double form = bcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion);
             inner_sum += GAUSS_W[j] * form;
         }
@@ -103,15 +103,15 @@
     const double Zq = outer_sum/(4.0*M_PI);
     const double Pq = sphere_form(q, radius, sld, solvent_sld);
-    return bcc_volume_fraction(radius, dnn) * Pq * Zq;
+    return bcc_volume_fraction(radius, lattice_spacing) * Pq * Zq;
 }
 
 
 static double Iqabc(double qa, double qb, double qc,
-    double dnn, double d_factor, double radius,
+    double lattice_spacing, double lattice_distortion, double radius,
     double sld, double solvent_sld)
 {
     const double q = sqrt(qa*qa + qb*qb + qc*qc);
-    const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor);
+    const double Zq = bcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion);
     const double Pq = sphere_form(q, radius, sld, solvent_sld);
-    return bcc_volume_fraction(radius, dnn) * Pq * Zq;
+    return bcc_volume_fraction(radius, lattice_spacing) * Pq * Zq;
 }
Index: sasmodels/models/bcc_paracrystal.py
===================================================================
--- sasmodels/models/bcc_paracrystal.py	(revision da7b26b2f3ecdd4d6eea063846bb148e3def2477)
+++ sasmodels/models/bcc_paracrystal.py	(revision da7b26b2f3ecdd4d6eea063846bb148e3def2477)
@@ -34,5 +34,5 @@
 .. math::
 
-    V_\text{lattice} = \frac{16\pi}{3} \frac{R^3}{\left(D\sqrt{2}\right)^3}
+    V_\text{lattice} = \frac{8\pi}{3} \frac{R^3}{\left(2D/\sqrt{3}\right)^3}
 
 
@@ -104,6 +104,6 @@
 
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
-* **Last Modified by:** Paul Butler **Date:** September 29, 2016
-* **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016
+* **Last Modified by:** Paul Butler **Date:** September 16, 2018
+* **Last Reviewed by:** Paul Butler **Date:** September 16, 2018
 """
 
@@ -127,6 +127,6 @@
 # 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"],
+parameters = [["lattice_spacing",         "Ang",       220,    [-inf, inf], "",            "Lattice spacing"],
+              ["lattice_distortion",    "",            0.06, [-inf, inf], "",            "Paracrystal distortion factor"],
               ["radius",      "Ang",        40,    [0, inf],    "volume",      "Particle radius"],
               ["sld",         "1e-6/Ang^2",  4,    [-inf, inf], "sld",         "Particle scattering length density"],
@@ -149,11 +149,11 @@
     # in this range 'cuz its easy.
     radius = 10**np.random.uniform(1.3, 4)
-    d_factor = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7
-    dnn_fraction = np.random.beta(a=10, b=1)
-    dnn = radius*4/np.sqrt(3)/dnn_fraction
+    lattice_distortion = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7
+    lattice_spacing_fraction = np.random.beta(a=10, b=1)
+    lattice_spacing = radius*4/np.sqrt(3)/lattice_spacing_fraction
     pars = dict(
         #sld=1, sld_solvent=0, scale=1, background=1e-32,
-        dnn=dnn,
-        d_factor=d_factor,
+        lattice_spacing=lattice_spacing,
+        lattice_distortion=lattice_distortion,
         radius=radius,
     )
Index: sasmodels/models/fcc_paracrystal.c
===================================================================
--- sasmodels/models/fcc_paracrystal.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/fcc_paracrystal.c	(revision 642046e89e9bce7059d153e0f8506ea74d7c1420)
@@ -1,4 +1,4 @@
 static double
-fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
+fcc_Zq(double qa, double qb, double qc, double lattice_spacing, double lattice_distortion)
 {
     // Equations from Matsuoka 17-18-19, multiplied by |q|
@@ -16,10 +16,10 @@
     //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1)
     //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1
-    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    const double arg = -0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3);
     const double exp_arg = exp(arg);
     const double Zq = -cube(expm1(2.0*arg))
-        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0)
-          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0)
-          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0));
+        / ( ((exp_arg - 2.0*cos(lattice_spacing*a1))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(lattice_spacing*a2))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(lattice_spacing*a3))*exp_arg + 1.0));
 
     return Zq;
@@ -29,7 +29,7 @@
 // occupied volume fraction calculated from lattice symmetry and sphere radius
 static double
-fcc_volume_fraction(double radius, double dnn)
+fcc_volume_fraction(double radius, double lattice_spacing)
 {
-    return 4.0*sphere_volume(M_SQRT1_2*radius/dnn);
+    return 4.0*sphere_volume(radius/lattice_spacing);
 }
 
@@ -41,6 +41,6 @@
 
 
-static double Iq(double q, double dnn,
-  double d_factor, double radius,
+static double Iq(double q, double lattice_spacing,
+  double lattice_distortion, double radius,
   double sld, double solvent_sld)
 {
@@ -66,5 +66,5 @@
             const double qa = qab*cos_phi;
             const double qb = qab*sin_phi;
-            const double form = fcc_Zq(qa, qb, qc, dnn, d_factor);
+            const double form = fcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion);
             inner_sum += GAUSS_W[j] * form;
         }
@@ -76,14 +76,14 @@
     const double Pq = sphere_form(q, radius, sld, solvent_sld);
 
-    return fcc_volume_fraction(radius, dnn) * Pq * Zq;
+    return fcc_volume_fraction(radius, lattice_spacing) * Pq * Zq;
 }
 
 static double Iqabc(double qa, double qb, double qc,
-    double dnn, double d_factor, double radius,
+    double lattice_spacing, double lattice_distortion, double radius,
     double sld, double solvent_sld)
 {
     const double q = sqrt(qa*qa + qb*qb + qc*qc);
     const double Pq = sphere_form(q, radius, sld, solvent_sld);
-    const double Zq = fcc_Zq(qa, qb, qc, dnn, d_factor);
-    return fcc_volume_fraction(radius, dnn) * Pq * Zq;
+    const double Zq = fcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion);
+    return fcc_volume_fraction(radius, lattice_spacing) * Pq * Zq;
 }
Index: sasmodels/models/fcc_paracrystal.py
===================================================================
--- sasmodels/models/fcc_paracrystal.py	(revision da7b26b2f3ecdd4d6eea063846bb148e3def2477)
+++ sasmodels/models/fcc_paracrystal.py	(revision da7b26b2f3ecdd4d6eea063846bb148e3def2477)
@@ -3,7 +3,7 @@
 #note - calculation requires double precision
 r"""
-.. warning:: This model and this model description are under review following 
-             concerns raised by SasView users. If you need to use this model, 
-             please email help@sasview.org for the latest situation. *The 
+.. warning:: This model and this model description are under review following
+             concerns raised by SasView users. If you need to use this model,
+             please email help@sasview.org for the latest situation. *The
              SasView Developers. September 2018.*
 
@@ -100,9 +100,9 @@
 
 Authorship and Verification
----------------------------
+----------------------------
 
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
-* **Last Modified by:** Paul Butler **Date:** September 29, 2016
-* **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016
+* **Last Modified by:** Paul Butler **Date:** September 16, 2018
+* **Last Reviewed by:** Paul Butler **Date:** September 16, 2018
 """
 
@@ -123,6 +123,6 @@
 # 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"],
+parameters = [["lattice_spacing", "Ang", 220, [-inf, inf], "", "Lattice spacing"],
+              ["lattice_distortion", "", 0.06, [-inf, inf], "", "Paracrystal distortion factor"],
               ["radius", "Ang", 40, [0, inf], "volume", "Particle radius"],
               ["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Particle scattering length density"],
@@ -139,11 +139,11 @@
     # copied from bcc_paracrystal
     radius = 10**np.random.uniform(1.3, 4)
-    d_factor = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7
-    dnn_fraction = np.random.beta(a=10, b=1)
-    dnn = radius*4/np.sqrt(2)/dnn_fraction
+    lattice_distortion = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7
+    lattice_spacing_fraction = np.random.beta(a=10, b=1)
+    lattice_spacing = radius*4/np.sqrt(2)/lattice_spacing_fraction
     pars = dict(
         #sld=1, sld_solvent=0, scale=1, background=1e-32,
-        dnn=dnn,
-        d_factor=d_factor,
+        latice_spacing=lattice_spacing,
+        lattice_distortion=d_factor,
         radius=radius,
     )
Index: sasmodels/models/sc_paracrystal.c
===================================================================
--- sasmodels/models/sc_paracrystal.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/sc_paracrystal.c	(revision 6530963bb96348a84504e7d230b909ede701366c)
@@ -1,4 +1,4 @@
 static double
-sc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
+sc_Zq(double qa, double qb, double qc, double lattice_spacing, double lattice_distortion)
 {
     // Equations from Matsuoka 9-10-11, multiplied by |q|
@@ -16,10 +16,10 @@
     //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1)
     //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1
-    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    const double arg = -0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3);
     const double exp_arg = exp(arg);
     const double Zq = -cube(expm1(2.0*arg))
-        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0)
-          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0)
-          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0));
+        / ( ((exp_arg - 2.0*cos(lattice_spacing*a1))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(lattice_spacing*a2))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(lattice_spacing*a3))*exp_arg + 1.0));
 
     return Zq;
@@ -28,7 +28,7 @@
 // occupied volume fraction calculated from lattice symmetry and sphere radius
 static double
-sc_volume_fraction(double radius, double dnn)
+sc_volume_fraction(double radius, double lattice_spacing)
 {
-    return sphere_volume(radius/dnn);
+    return sphere_volume(radius/lattice_spacing);
 }
 
@@ -41,6 +41,6 @@
 
 static double
-Iq(double q, double dnn,
-    double d_factor, double radius,
+Iq(double q, double lattice_spacing,
+    double lattice_distortion, double radius,
     double sld, double solvent_sld)
 {
@@ -67,5 +67,5 @@
             const double qa = qab*cos_phi;
             const double qb = qab*sin_phi;
-            const double form = sc_Zq(qa, qb, qc, dnn, d_factor);
+            const double form = sc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion);
             inner_sum += GAUSS_W[j] * form;
         }
@@ -77,15 +77,15 @@
     const double Pq = sphere_form(q, radius, sld, solvent_sld);
 
-    return sc_volume_fraction(radius, dnn) * Pq * Zq;
+    return sc_volume_fraction(radius, lattice_spacing) * Pq * Zq;
 }
 
 static double
 Iqabc(double qa, double qb, double qc,
-    double dnn, double d_factor, double radius,
+    double lattice_spacing, double lattice_distortion, double radius,
     double sld, double solvent_sld)
 {
     const double q = sqrt(qa*qa + qb*qb + qc*qc);
     const double Pq = sphere_form(q, radius, sld, solvent_sld);
-    const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor);
-    return sc_volume_fraction(radius, dnn) * Pq * Zq;
+    const double Zq = sc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion);
+    return sc_volume_fraction(radius, lattice_spacing) * Pq * Zq;
 }
Index: sasmodels/models/sc_paracrystal.py
===================================================================
--- sasmodels/models/sc_paracrystal.py	(revision da7b26b2f3ecdd4d6eea063846bb148e3def2477)
+++ sasmodels/models/sc_paracrystal.py	(revision da7b26b2f3ecdd4d6eea063846bb148e3def2477)
@@ -1,8 +1,8 @@
 r"""
-.. warning:: This model and this model description are under review following 
-             concerns raised by SasView users. If you need to use this model, 
-             please email help@sasview.org for the latest situation. *The 
+.. warning:: This model and this model description are under review following
+             concerns raised by SasView users. If you need to use this model,
+             please email help@sasview.org for the latest situation. *The
              SasView Developers. September 2018.*
-             
+
 Definition
 ----------
@@ -104,9 +104,9 @@
 
 Authorship and Verification
----------------------------
+----------------------------
 
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
-* **Last Modified by:** Paul Butler **Date:** September 29, 2016
-* **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016
+* **Last Modified by:** Paul Butler **Date:** September 16, 2018
+* **Last Reviewed by:** Paul Butler **Date:** September 16, 2018
 """
 
@@ -129,6 +129,6 @@
         scale: volume fraction of spheres
         bkg:background, R: radius of sphere
-        dnn: Nearest neighbor distance
-        d_factor: Paracrystal distortion factor
+        lattice_spacing: Nearest neighbor distance
+        lattice_distortion: Paracrystal distortion factor
         radius: radius of the spheres
         sldSph: SLD of the sphere
@@ -139,6 +139,6 @@
 # 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"],
+parameters = [["lattice_spacing",         "Ang",       220.0, [0.0, inf],  "",  "Lattice spacing"],
+              ["lattice_distortion",    "",           0.06, [-inf, inf], "",   "Paracrystal distortion factor"],
               ["radius",      "Ang",        40.0, [0.0, inf],  "volume",      "Radius of sphere"],
               ["sld",  "1e-6/Ang^2",         3.0, [0.0, inf],  "sld",         "Sphere scattering length density"],
@@ -155,11 +155,11 @@
     # copied from bcc_paracrystal
     radius = 10**np.random.uniform(1.3, 4)
-    d_factor = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7
-    dnn_fraction = np.random.beta(a=10, b=1)
-    dnn = radius*4/np.sqrt(4)/dnn_fraction
+    lattice_distortion = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7
+    lattice_spacing_fraction = np.random.beta(a=10, b=1)
+    lattice_spacing = radius*4/np.sqrt(4)/lattice_spacing_fraction
     pars = dict(
         #sld=1, sld_solvent=0, scale=1, background=1e-32,
-        dnn=dnn,
-        d_factor=d_factor,
+        lattice_spacing=lattice_spacing,
+        lattice_distortion=lattice_distortion,
         radius=radius,
     )
