Index: sasmodels/models/dab.py
===================================================================
--- sasmodels/models/dab.py	(revision 4962519cc1910422a467e4f643dba223623e0e09)
+++ sasmodels/models/dab.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -61,7 +61,17 @@
     double numerator   = cube(cor_length);
     double denominator = square(1 + square(q*cor_length));
-    
+
     return numerator / denominator ;
     """
+
+def random():
+    import numpy as np
+    pars = dict(
+        scale=10**np.random.uniform(1, 4),
+        cor_length=10**np.random.uniform(0.3, 3),
+        #background = 0,
+    )
+    pars['scale'] /= pars['cor_length']**3
+    return pars
 
 # ER defaults to 1.0
Index: sasmodels/models/ellipsoid.py
===================================================================
--- sasmodels/models/ellipsoid.py	(revision 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
+++ sasmodels/models/ellipsoid.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -2,5 +2,5 @@
 # Note: model title and parameter table are inserted automatically
 r"""
-The form factor is normalized by the particle volume 
+The form factor is normalized by the particle volume
 
 Definition
@@ -119,4 +119,5 @@
 * **Last Modified by:** Paul Kienzle **Date:** March 22, 2017
 """
+from __future__ import division
 
 from numpy import inf, sin, cos, pi
@@ -177,8 +178,21 @@
     delta = 0.75 * b1 * b2
 
-    ddd = np.zeros_like(radius_polar)
+    #ddd = np.zeros_like(radius_polar)
     ddd[valid] = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial ** 2
     return 0.5 * ddd ** (1.0 / 3.0)
 
+def random():
+    import numpy as np
+    V = 10**np.random.uniform(4, 12)
+    radius_polar = 10**np.random.uniform(1.3, 4)
+    radius_equatorial = np.sqrt(V/radius_polar) # ignore 4/3 pi
+    Vf = 10**np.random.uniform(-4, -2)
+    pars = dict(
+        #background=0, sld=0, sld_solvent=1,
+        scale=1e9*Vf/V,
+        radius_polar=radius_polar,
+        radius_equatorial=radius_equatorial,
+    )
+    return pars
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/elliptical_cylinder.py
===================================================================
--- sasmodels/models/elliptical_cylinder.py	(revision 9802ab3945a154b871264715dcb463d4ec644afa)
+++ sasmodels/models/elliptical_cylinder.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -33,7 +33,7 @@
 
     a = qr'\sin(\alpha)
-    
+
     b = q\frac{L}{2}\cos(\alpha)
-    
+
     r'=\frac{r_{minor}}{\sqrt{2}}\sqrt{(1+\nu^{2}) + (1-\nu^{2})cos(\psi)}
 
@@ -57,5 +57,5 @@
 define the axis of the cylinder using two angles $\theta$, $\phi$ and $\Psi$
 (see :ref:`cylinder orientation <cylinder-angle-definition>`). The angle
-$\Psi$ is the rotational angle around its own long_c axis. 
+$\Psi$ is the rotational angle around its own long_c axis.
 
 All angle parameters are valid and given only for 2D calculation; ie, an
@@ -72,11 +72,11 @@
     detector plane, with $\Psi$ = 0.
 
-The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data. 
+The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.
 On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will
-appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, the $b$ and $a$ axes of the 
-cylinder cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.) 
-The third orientation distribution, in $\psi$, is about the $c$ axis of the particle. Some experimentation may be required to 
-understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation 
-distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.) 
+appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, the $b$ and $a$ axes of the
+cylinder cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.)
+The third orientation distribution, in $\psi$, is about the $c$ axis of the particle. Some experimentation may be required to
+understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation
+distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.)
 
 NB: The 2nd virial coefficient of the cylinder is calculated based on the
@@ -110,5 +110,5 @@
 
 * **Author:**
-* **Last Modified by:** 
+* **Last Modified by:**
 * **Last Reviewed by:**  Richard Heenan - corrected equation in docs **Date:** December 21, 2016
 
@@ -156,4 +156,22 @@
                            + (length + radius) * (length + pi * radius))
     return 0.5 * (ddd) ** (1. / 3.)
+
+def random():
+    import numpy as np
+    # V = pi * radius_major * radius_minor * length;
+    V = 10**np.random.uniform(3, 9)
+    length = 10**np.random.uniform(1, 3)
+    axis_ratio = 10**(4*np.random.uniform(-2, 2)
+    radius_minor = np.sqrt(V/length/axis_ratio)
+    Vf = 10**np.random.uniform(-4, -2)
+    pars = dict(
+        #background=0, sld=0, sld_solvent=1,
+        scale=1e9*Vf/V,
+        length=length,
+        radius_minor=radius_minor,
+        axis_ratio=axis_ratio,
+    )
+    return pars
+
 q = 0.1
 # april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct!
Index: sasmodels/models/fcc_paracrystal.py
===================================================================
--- sasmodels/models/fcc_paracrystal.py	(revision 69e1afc57270f67b2a0d7157a6d263fc2db98cc0)
+++ sasmodels/models/fcc_paracrystal.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -78,5 +78,5 @@
 .. figure:: img/parallelepiped_angle_definition.png
 
-    Orientation of the crystal with respect to the scattering plane, when 
+    Orientation of the crystal with respect to the scattering plane, when
     $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis).
 
@@ -119,4 +119,24 @@
 source = ["lib/sas_3j1x_x.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc_paracrystal.c"]
 
+def random():
+    import numpy as np
+    # Define lattice spacing as a multiple of the particle radius
+    # using the formulat a = 4 r/sqrt(3).  Systems which are ordered
+    # are probably mostly filled, so use a distribution which goes from
+    # zero to one, but leaving 90% of them within 80% of the
+    # maximum bcc packing.  Lattice distortion values are empirically
+    # useful between 0.01 and 0.7.  Use an exponential distribution
+    # in this range 'cuz its easy.
+    dnn_fraction = np.random.beta(a=10, b=1)
+    pars = dict(
+        #sld=1, sld_solvent=0, scale=1, background=1e-32,
+        radius=10**np.random.uniform(1.3, 4),
+        d_factor=10**np.random.uniform(-2, -0.7),  # sigma_d in 0.01-0.7
+    )
+    pars['dnn'] = pars['radius']*4/np.sqrt(3)/dnn_fraction
+    #pars['scale'] = 1/(0.68*dnn_fraction**3)  # bcc packing fraction is 0.68
+    pars['scale'] = 1
+    return pars
+
 # parameters for demo
 demo = dict(scale=1, background=0,
Index: sasmodels/models/fractal.py
===================================================================
--- sasmodels/models/fractal.py	(revision df89d77375bb6b013466ad838acfb111b7a0ab90)
+++ sasmodels/models/fractal.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -27,6 +27,6 @@
 
 where $\xi$ is the correlation length representing the cluster size and $D_f$
-is the fractal dimension, representing the self similarity of the structure. 
-Note that S(q) here goes negative if $D_f$ is too large, and the Gamma function 
+is the fractal dimension, representing the self similarity of the structure.
+Note that S(q) here goes negative if $D_f$ is too large, and the Gamma function
 diverges at $D_f=0$ and $D_f=1$.
 
@@ -55,4 +55,5 @@
 
 """
+from __future__ import division
 
 from numpy import inf
@@ -98,4 +99,24 @@
 source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "lib/fractal_sq.c", "fractal.c"]
 
+def random():
+    import numpy as np
+    radius = 10**np.random.uniform(0.7, 4)
+    #radius = 5
+    cor_length = 10**np.random.uniform(0.7, 2)*radius
+    #cor_length = 20*radius
+    volfraction = 10**np.random.uniform(-3, -1)
+    #volfraction = 0.05
+    fractal_dim = 2*np.random.beta(3, 4) + 1
+    #fractal_dim = 2
+    pars = dict(
+        #background=0, sld_block=1, sld_solvent=0,
+        scale=1e4/radius**(fractal_dim/2),
+        volfraction=volfraction,
+        radius=radius,
+        cor_length=cor_length,
+        fractal_dim=fractal_dim,
+    )
+    return pars
+
 demo = dict(volfraction=0.05,
             radius=5.0,
Index: sasmodels/models/fractal_core_shell.py
===================================================================
--- sasmodels/models/fractal_core_shell.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/fractal_core_shell.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -24,5 +24,5 @@
     \frac{\sin(qr_s)-qr_s\cos(qr_s)}{(qr_s)^3}\right]^2
 
-    S(q) &= 1 + \frac{D_f\ \Gamma\!(D_f-1)}{[1+1/(q\xi)^2]^{(D_f-1)/2}} 
+    S(q) &= 1 + \frac{D_f\ \Gamma\!(D_f-1)}{[1+1/(q\xi)^2]^{(D_f-1)/2}}
     \frac{\sin[(D_f-1)\tan^{-1}(q\xi)]}{(qr_s)^{D_f}}
 
@@ -33,5 +33,5 @@
 of the whole particle respectively, $D_f$ is the fractal dimension, and |xi| the
 correlation length.
- 
+
 Polydispersity of radius and thickness are also provided for.
 
@@ -98,4 +98,28 @@
           "lib/fractal_sq.c", "fractal_core_shell.c"]
 
+def random():
+    import numpy as np
+    total_radius = 10**np.random.uniform(0.7, 4)
+    core_portion = np.random.uniform(0, 1)
+    radius = total_radius * core_portion
+    thickness = total_radius - radius
+    #radius = 5
+    cor_length = 10**np.random.uniform(0.7, 2)*radius
+    #cor_length = 20*radius
+    volfraction = 10**np.random.uniform(-3, -1)
+    #volfraction = 0.05
+    fractal_dim = 2*np.random.beta(3, 4) + 1
+    #fractal_dim = 2
+    pars = dict(
+        #background=0, sld_block=1, sld_solvent=0,
+        scale=1e3/total_radius**(fractal_dim/2),
+        volfraction=volfraction,
+        radius=radius,
+        cor_length=cor_length,
+        fractal_dim=fractal_dim,
+    )
+    return pars
+
+
 demo = dict(scale=0.05,
             background=0,
Index: sasmodels/models/guinier.py
===================================================================
--- sasmodels/models/guinier.py	(revision 3330bb44faf647c84e605c5a75bed3f385a25a52)
+++ sasmodels/models/guinier.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -33,5 +33,5 @@
 description = """
  I(q) = scale.exp ( - rg^2 q^2 / 3.0 )
- 
+
     List of default parameters:
     scale = scale
@@ -49,4 +49,18 @@
 """
 
+def random():
+    import numpy as np
+    scale = 10**np.random.uniform(1, 5)
+    # Note: compare.py has rg cutoff for guinier, so use that
+    q_max = 1.0
+    rg_max = np.sqrt(90*np.log(10) + 3*np.log(scale))/q_max
+    rg = 10**np.random.uniform(0, np.log10(rg_max))
+    pars = dict(
+        #background=0,
+        scale=scale,
+        rg=rg,
+    )
+    return pars
+
 # parameters for demo
 demo = dict(scale=1.0, rg=60.0)
Index: sasmodels/models/lamellar.py
===================================================================
--- sasmodels/models/lamellar.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/models/lamellar.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -89,4 +89,14 @@
     """
 
+def random():
+    import numpy as np
+    thickness = 10**np.random.uniform(1, 4)
+    scale = thickness * 10**np.random.uniform(-7, -4)
+    pars = dict(
+        scale=scale,
+        thickness=thickness,
+    )
+    return pars
+
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/lamellar_hg.py
===================================================================
--- sasmodels/models/lamellar_hg.py	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/lamellar_hg.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -98,4 +98,17 @@
     """
 
+def random():
+    import numpy as np
+    thickness = 10**np.random.uniform(1, 4)
+    length_head = thickness * np.random.uniform(0, 1)
+    length_tail = thickness - length_head
+    scale = thickness * 10**np.random.uniform(-7, -4)
+    pars = dict(
+        scale=scale,
+        length_head=length_head,
+        length_tail=length_tail,
+    )
+    return pars
+
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/lamellar_hg_stack_caille.py
===================================================================
--- sasmodels/models/lamellar_hg_stack_caille.py	(revision a57b31dc1634e8d710e6f72d02be91edfc5b23eb)
+++ sasmodels/models/lamellar_hg_stack_caille.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -123,4 +123,24 @@
 # VR defaults to 1.0
 
+def random():
+    import numpy as np
+    total_thickness = 10**np.random.uniform(2, 4.7)
+    Nlayers = np.random.randint(2, 200)
+    d_spacing = total_thickness / Nlayers
+    thickness = d_spacing * np.random.uniform(0, 1)
+    length_head = thickness * np.random.uniform(0, 1)
+    length_tail = thickness - length_head
+    Caille_parameter = np.random.uniform(0, 0.8)
+    scale = thickness * 10**np.random.uniform(-7, -4)
+    pars = dict(
+        scale=1,
+        length_head=length_head,
+        length_tail=length_tail,
+        Nlayers=Nlayers,
+        d_spacing=d_spacing,
+        Caille_parameter=Caille_parameter,
+    )
+    return pars
+
 demo = dict(
     scale=1, background=0,
Index: sasmodels/models/lamellar_stack_caille.py
===================================================================
--- sasmodels/models/lamellar_stack_caille.py	(revision a57b31dc1634e8d710e6f72d02be91edfc5b23eb)
+++ sasmodels/models/lamellar_stack_caille.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -98,4 +98,21 @@
 source = ["lamellar_stack_caille.c"]
 
+def random():
+    import numpy as np
+    total_thickness = 10**np.random.uniform(2, 4.7)
+    Nlayers = np.random.randint(2, 200)
+    d_spacing = total_thickness / Nlayers
+    thickness = d_spacing * np.random.uniform(0, 1)
+    Caille_parameter = np.random.uniform(0, 0.8)
+    scale = thickness * 10**np.random.uniform(-7, -4)
+    pars = dict(
+        scale=1,
+        thickness=thickness,
+        Nlayers=Nlayers,
+        d_spacing=d_spacing,
+        Caille_parameter=Caille_parameter,
+    )
+    return pars
+
 # No volume normalization despite having a volume parameter
 # This should perhaps be volume normalized?
Index: sasmodels/models/lamellar_stack_paracrystal.py
===================================================================
--- sasmodels/models/lamellar_stack_paracrystal.py	(revision a0168e84cf86a1efabf5abe8430ab9266ddf1879)
+++ sasmodels/models/lamellar_stack_paracrystal.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -130,5 +130,23 @@
 form_volume = """
     return 1.0;
-    """
+"""
+
+def random():
+    import numpy as np
+    total_thickness = 10**np.random.uniform(2, 4.7)
+    Nlayers = np.random.randint(2, 200)
+    d_spacing = total_thickness / Nlayers
+    thickness = d_spacing * np.random.uniform(0, 1)
+    # Let polydispersity peak around 15%; 95% < 0.4; max=100%
+    sigma_d = np.random.beta(1.5, 7)
+    scale = thickness * 10**np.random.uniform(-7, -4)
+    pars = dict(
+        scale=1,
+        thickness=thickness,
+        Nlayers=Nlayers,
+        d_spacing=d_spacing,
+        sigma_d=sigma_d,
+    )
+    return pars
 
 # ER defaults to 0.0
Index: sasmodels/models/line.py
===================================================================
--- sasmodels/models/line.py	(revision 3330bb44faf647c84e605c5a75bed3f385a25a52)
+++ sasmodels/models/line.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -68,4 +68,15 @@
 Iqxy.vectorized = True  # Iqxy accepts an array of qx qy values
 
+def random():
+    import numpy as np
+    slope = 1e5*np.random.uniform(-1, 1)
+    offset = 1e-5 + (0 if slope > 0 else -slope)
+    intercept = 10**np.random.uniform(0, 5) + offset
+    pars = dict(
+        scale=1, background=0,
+        slope=slope,
+        intercept=intercept,
+    )
+    return pars
 
 tests = [
Index: sasmodels/models/lorentz.py
===================================================================
--- sasmodels/models/lorentz.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
+++ sasmodels/models/lorentz.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -30,11 +30,11 @@
 description = """
 Model that evaluates a Lorentz (Ornstein-Zernicke) model.
-        
+
 I(q) = scale/( 1 + (q*L)^2 ) + bkd
-        
-The model has three parameters: 
-    length     =  screening Length
-    scale  =  scale factor
-    background    =  incoherent background
+
+The model has three parameters:
+    length = screening Length
+    scale = scale factor
+    background = incoherent background
 """
 category = "shape-independent"
@@ -48,4 +48,13 @@
 """
 
+def random():
+    import numpy as np
+    pars = dict(
+        #background=0,
+        scale=10**np.random.uniform(1, 4),
+        cor_length=10**np.random.uniform(0, 3),
+    )
+    return pars
+
 # parameters for demo
 demo = dict(scale=1.0, background=0.0, cor_length=50.0)
Index: sasmodels/models/mass_fractal.py
===================================================================
--- sasmodels/models/mass_fractal.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/mass_fractal.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -88,4 +88,19 @@
 source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "mass_fractal.c"]
 
+def random():
+    import numpy as np
+    radius = 10**np.random.uniform(0.7, 4)
+    cutoff_length = 10**np.random.uniform(0.7, 2)*radius
+    fractal_dim_mass = 2*np.random.beta(3, 4) + 1
+    Vf = 10**np.random.uniform(-4, -1)
+    pars = dict(
+        #background=0,
+        scale=1, #1e5*Vf/radius**(fractal_dim_mass),
+        radius=radius,
+        cutoff_length=cutoff_length,
+        fractal_dim_mass=fractal_dim_mass,
+    )
+    return pars
+
 demo = dict(scale=1, background=0,
             radius=10.0,
Index: sasmodels/models/mass_surface_fractal.py
===================================================================
--- sasmodels/models/mass_surface_fractal.py	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
+++ sasmodels/models/mass_surface_fractal.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -89,4 +89,24 @@
 source = ["mass_surface_fractal.c"]
 
+def random():
+    import numpy as np
+    fractal_dim = np.random.uniform(1, 6)
+    surface_portion = np.random.uniform(0, 1)
+    fractal_dim_surf = fractal_dim*surface_portion
+    fractal_dim_mass = fractal_dim - fractal_dim_surf
+    rg_cluster = 10**np.random.uniform(1, 5)
+    rg_primary = rg_cluster*10**np.random.uniform(-4, -1)
+    scale = 10**np.random.uniform(2, 5)
+    pars = dict(
+        #background=0,
+        scale=scale,
+        fractal_dim_mass=fractal_dim_mass,
+        fractal_dim_surf=fractal_dim_surf,
+        rg_cluster=rg_cluster,
+        rg_primary=rg_primary,
+    )
+    return pars
+
+
 demo = dict(scale=1, background=0,
             fractal_dim_mass=1.8,
Index: sasmodels/models/mono_gauss_coil.py
===================================================================
--- sasmodels/models/mono_gauss_coil.py	(revision 3330bb44faf647c84e605c5a75bed3f385a25a52)
+++ sasmodels/models/mono_gauss_coil.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -62,5 +62,5 @@
 
 description = """
-    Evaluates the scattering from 
+    Evaluates the scattering from
     monodisperse polymer chains.
     """
@@ -86,4 +86,15 @@
 Iq.vectorized = True # Iq accepts an array of q values
 
+def random():
+    import numpy as np
+    rg = 10**np.random.uniform(0, 4)
+    #rg = 1e3
+    pars = dict(
+        #scale=1, background=0,
+        i_zero=1e7, # i_zero is a simple scale
+        rg=rg,
+    )
+    return pars
+
 demo = dict(scale=1.0, i_zero=70.0, rg=75.0, background=0.0)
 
Index: sasmodels/models/peak_lorentz.py
===================================================================
--- sasmodels/models/peak_lorentz.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
+++ sasmodels/models/peak_lorentz.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -59,4 +59,16 @@
 Iq.vectorized = True  # Iq accepts an array of q values
 
+def random():
+    import numpy as np
+    peak_pos = 10**np.random.uniform(-3, -1)
+    peak_hwhm = peak_pos * 10**np.random.uniform(-3, 0)
+    pars = dict(
+        #background=0,
+        scale=10**np.random.uniform(2, 6),
+        peak_pos=peak_pos,
+        peak_hwhm=peak_hwhm,
+    )
+    return pars
+
 demo = dict(scale=100, background=1.0,
             peak_pos=0.05, peak_hwhm=0.005)
Index: sasmodels/models/poly_gauss_coil.py
===================================================================
--- sasmodels/models/poly_gauss_coil.py	(revision 3330bb44faf647c84e605c5a75bed3f385a25a52)
+++ sasmodels/models/poly_gauss_coil.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -66,5 +66,5 @@
 
 description = """
-    Evaluates the scattering from 
+    Evaluates the scattering from
     polydisperse polymer chains.
     """
@@ -96,9 +96,9 @@
         p = [
             #(-1 - 20*u - 155*u**2 - 580*u**3 - 1044*u**4 - 720*u**5) / 2520.,
-            #( 1 + 14*u + 71*u**2 + 154*u**3 + 120*u**4) / 360.,
+            #(+1 + 14*u + 71*u**2 + 154*u**3 + 120*u**4) / 360.,
             #(-1 - 9*u - 26*u**2 - 24*u**3) / 60.,
-            ( 1 + 5*u + 6*u**2) / 12.,
+            (+1 + 5*u + 6*u**2) / 12.,
             (-1 - 2*u) / 3.,
-            ( 1 ),
+            (+1),
             ]
         result = 2.0 * (power(1.0 + u*z, -1.0/u) + z - 1.0) / (1.0 + u)
@@ -108,4 +108,17 @@
     return i_zero * result
 Iq.vectorized = True  # Iq accepts an array of q values
+
+def random():
+    import numpy as np
+    rg = 10**np.random.uniform(0, 4)
+    #rg = 1e3
+    polydispersity = 10**np.random.uniform(0, 3)
+    pars = dict(
+        #scale=1, background=0,
+        i_zero=1e7, # i_zero is a simple scale
+        rg=rg,
+        polydispersity=polydispersity,
+    )
+    return pars
 
 demo = dict(scale=1.0,
Index: sasmodels/models/polymer_excl_volume.py
===================================================================
--- sasmodels/models/polymer_excl_volume.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/models/polymer_excl_volume.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -108,6 +108,6 @@
 #   ["name", "units", default, [lower, upper], "type", "description"],
 parameters = [
-    ["rg",        "Ang", 60.0, [0, inf],    "", "Radius of Gyration"],
-    ["porod_exp", "",     3.0, [-inf, inf], "", "Porod exponent"],
+    ["rg",        "Ang", 60.0, [0, inf], "", "Radius of Gyration"],
+    ["porod_exp", "",     3.0, [0, inf], "", "Porod exponent"],
 ]
 # pylint: enable=bad-whitespace, line-too-long
@@ -133,4 +133,17 @@
 Iq.vectorized = True  # Iq accepts an array of q values
 
+def random():
+    import numpy as np
+    rg = 10**np.random.uniform(0, 4)
+    porod_exp = np.random.uniform(1e-3, 6)
+    scale = 10**np.random.uniform(1, 5)
+    pars = dict(
+        #background=0,
+        scale=scale,
+        rg=rg,
+        porod_exp=porod_exp,
+    )
+    return pars
+
 tests = [
     # Accuracy tests based on content in test/polyexclvol_default_igor.txt
Index: sasmodels/models/polymer_micelle.py
===================================================================
--- sasmodels/models/polymer_micelle.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/polymer_micelle.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -16,5 +16,5 @@
 which then defines a micelle core of $radius\_core$, which is a separate parameter
 even though it could be directly determined.
-The Gaussian random coil tails, of gyration radius $rg$, are imagined uniformly 
+The Gaussian random coil tails, of gyration radius $rg$, are imagined uniformly
 distributed around the spherical core, centred at a distance $radius\_core + d\_penetration.rg$
 from the micelle centre, where $d\_penetration$ is of order unity.
@@ -27,12 +27,12 @@
 .. math::
     P(q) = N^2\beta^2_s\Phi(qR)^2+N\beta^2_cP_c(q)+2N^2\beta_s\beta_cS_{sc}s_c(q)+N(N-1)\beta_c^2S_{cc}(q)
-    
+
     \beta_s = v\_core(sld\_core - sld\_solvent)
-    
+
     \beta_c = v\_corona(sld\_corona - sld\_solvent)
 
-where $N = n\_aggreg$, and for the spherical core of radius $R$ 
+where $N = n\_aggreg$, and for the spherical core of radius $R$
 
-.. math::   
+.. math::
    \Phi(qR)= \frac{\sin(qr) - qr\cos(qr)}{(qr)^3}
 
@@ -49,9 +49,9 @@
 
 .. math::
-   
+
    S_{sc}(q)=\Phi(qR)\psi(Z)\frac{sin(q(R+d.R_g))}{q(R+d.R_g)}
-   
+
    S_{cc}(q)=\psi(Z)^2\left[\frac{sin(q(R+d.R_g))}{q(R+d.R_g)} \right ]^2
-   
+
    \psi(Z)=\frac{[1-exp^{-Z}]}{Z}
 
@@ -60,5 +60,5 @@
 
 $P(q)$ above is multiplied by $ndensity$, and a units conversion of 10^{-13}, so $scale$
-is likely 1.0 if the scattering data is in absolute units. This model has not yet been 
+is likely 1.0 if the scattering data is in absolute units. This model has not yet been
 independently validated.
 
@@ -71,5 +71,5 @@
 """
 
-from numpy import inf
+from numpy import inf, pi
 
 name = "polymer_micelle"
@@ -80,5 +80,5 @@
 to block copolymer micelles. To work well the Gaussian chains must be much
 smaller than the core, which is often not the case.  Please study the
-reference to Pedersen and full documentation carefully. 
+reference to Pedersen and full documentation carefully.
     """
 
@@ -106,16 +106,31 @@
 source = ["lib/sas_3j1x_x.c", "polymer_micelle.c"]
 
-demo = dict(scale=1, background=0,
-            ndensity=8.94,
-            v_core=62624.0,
-            v_corona=61940.0,
-            sld_solvent=6.4,
-            sld_core=0.34,
-            sld_corona=0.8,
-            radius_core=45.0,
-            rg=20.0,
-            d_penetration=1.0,
-            n_aggreg=6.0)
-
+def random():
+    import numpy as np
+    radius_core = 10**np.random.uniform(1, 3)
+    rg = radius_core * 10**np.random.uniform(-2, -0.3)
+    d_penetration = np.random.randn()*0.05 + 1
+    n_aggreg = np.random.randint(3, 30)
+    # volume of head groups is the core volume over the number of groups,
+    # with a correction for packing fraction of the head groups.
+    v_core = 4*pi/3*radius_core**3/n_aggreg * 0.68
+    # Rg^2 for gaussian coil is a^2n/6 => a^2 = 6 Rg^2/n
+    # a=2r => r = Rg sqrt(3/2n)
+    # v = 4/3 pi r^3 n => v = 4/3 pi Rg^3 (3/2n)^(3/2) n = pi Rg^3 sqrt(6/n)
+    tail_segments = np.random.randint(6, 30)
+    v_corona = pi * rg**3 * np.sqrt(6/tail_segments)
+    V = 4*pi/3*(radius_core + rg)**3
+    pars = dict(
+        background=0,
+        scale=1e7/V,
+        ndensity=8.94,
+        v_core=v_core,
+        v_corona=v_corona,
+        radius_core=radius_core,
+        rg=rg,
+        d_penetration=d_penetration,
+        n_aggreg=n_aggreg,
+    )
+    return pars
 
 tests = [
Index: sasmodels/models/power_law.py
===================================================================
--- sasmodels/models/power_law.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/models/power_law.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -50,4 +50,13 @@
 Iq.vectorized = True  # Iq accepts an array of q values
 
+def random():
+    import numpy as np
+    power = np.random.uniform(1, 6)
+    pars = dict(
+        scale=0.1**power*10**np.random.uniform(-4, 2),
+        power=power,
+    )
+    return pars
+
 demo = dict(scale=1.0, power=4.0, background=0.0)
 
Index: sasmodels/models/sc_paracrystal.py
===================================================================
--- sasmodels/models/sc_paracrystal.py	(revision 69e1afc57270f67b2a0d7157a6d263fc2db98cc0)
+++ sasmodels/models/sc_paracrystal.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -138,4 +138,24 @@
 source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c", "lib/gauss150.c", "sc_paracrystal.c"]
 
+def random():
+    import numpy as np
+    # Define lattice spacing as a multiple of the particle radius
+    # using the formulat a = 4 r/sqrt(3).  Systems which are ordered
+    # are probably mostly filled, so use a distribution which goes from
+    # zero to one, but leaving 90% of them within 80% of the
+    # maximum bcc packing.  Lattice distortion values are empirically
+    # useful between 0.01 and 0.7.  Use an exponential distribution
+    # in this range 'cuz its easy.
+    dnn_fraction = np.random.beta(a=10, b=1)
+    pars = dict(
+        #sld=1, sld_solvent=0, scale=1, background=1e-32,
+        radius=10**np.random.uniform(1.3, 4),
+        d_factor=10**np.random.uniform(-2, -0.7),  # sigma_d in 0.01-0.7
+    )
+    pars['dnn'] = pars['radius']*4/np.sqrt(3)/dnn_fraction
+    #pars['scale'] = 1/(0.68*dnn_fraction**3)  # bcc packing fraction is 0.68
+    pars['scale'] = 1
+    return pars
+
 demo = dict(scale=1, background=0,
             dnn=220.0,
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/sphere.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
@@ -84,4 +84,16 @@
     return radius
 
+def random():
+    import numpy as np
+    Vf = 10**np.random.uniform(-4, -2)
+    radius = 10**np.random.uniform(1.3, 4)
+    V = radius**3
+    pars = dict(
+        #background=0, sld=1, sld_solvent=0,
+        scale=1e10*Vf/V,
+        radius=radius,
+    )
+    return pars
+
 # VR defaults to 1.0
 
