Index: sasmodels/models/gauss_lorentz_gel.py
===================================================================
--- sasmodels/models/gauss_lorentz_gel.py	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/gauss_lorentz_gel.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -3,5 +3,5 @@
 but typically a physical rather than chemical network.
 It is modeled as a sum of a low-q exponential decay (which happens to
-give a functional form similar to Guinier scattering, so interpret with 
+give a functional form similar to Guinier scattering, so interpret with
 care) plus a Lorentzian at higher-q values. See also the gel_fit model.
 
@@ -88,4 +88,21 @@
 
 
+def random():
+    import numpy as np
+    gauss_scale = 10**np.random.uniform(1, 3)
+    lorentz_scale = 10**np.random.uniform(1, 3)
+    cor_length_static = 10**np.random.uniform(0, 3)
+    cor_length_dynamic = 10**np.random.uniform(0, 3)
+    pars = dict(
+        #background=0,
+        scale=1,
+        gauss_scale=gauss_scale,
+        lorentz_scale=lorentz_scale,
+        cor_length_static=cor_length_static,
+        cor_length_dynamic=cor_length_dynamic,
+    )
+    return pars
+
+
 demo = dict(scale=1, background=0.1,
             gauss_scale=100.0,
Index: sasmodels/models/gaussian_peak.py
===================================================================
--- sasmodels/models/gaussian_peak.py	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/gaussian_peak.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -51,3 +51,14 @@
 # VR defaults to 1.0
 
-demo = dict(scale=1, background=0, peak_pos=0.05, sigma=0.005)
+def random():
+    import numpy as np
+    peak_pos = 10**np.random.uniform(-3, -1)
+    sigma = 10**np.random.uniform(-1.3, -0.3)*peak_pos
+    scale = 10**np.random.uniform(0, 4)
+    pars = dict(
+        #background=1e-8,
+        scale=scale,
+        peak_pos=peak_pos,
+        sigam=sigma,
+    )
+    return pars
Index: sasmodels/models/gel_fit.c
===================================================================
--- sasmodels/models/gel_fit.c	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/gel_fit.c	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -10,10 +10,10 @@
     double lorentzian_term = square(q*cor_length);
     lorentzian_term = 1.0 + ((fractal_dim + 1.0)/3.0)*lorentzian_term;
-    lorentzian_term = pow(lorentzian_term, (fractal_dim/2.0) );
+    lorentzian_term = pow(lorentzian_term, fractal_dim/2.0 );
 
     // Exponential Term
     ////////////////////////double d(x[i]*x[i]*rg*rg);
     double exp_term = square(q*rg);
-    exp_term = exp(-1.0*(exp_term/3.0));
+    exp_term = exp(-exp_term/3.0);
 
     // Scattering Law
Index: sasmodels/models/gel_fit.py
===================================================================
--- sasmodels/models/gel_fit.py	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/gel_fit.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -71,4 +71,22 @@
 source = ["gel_fit.c"]
 
+def random():
+    import numpy as np
+    guinier_scale = 10**np.random.uniform(1, 3)
+    lorentz_scale = 10**np.random.uniform(1, 3)
+    rg = 10**np.random.uniform(1, 5)
+    fractal_dim = np.random.uniform(0, 6)
+    cor_length = 10**np.random.uniform(0, 3)
+    pars = dict(
+        #background=0,
+        scale=1,
+        guinier_scale=guinier_scale,
+        lorentz_scale=lorentz_scale,
+        rg=rg,
+        fractal_dim=fractal_dim,
+        cor_length=cor_length
+    )
+    return pars
+
 demo = dict(background=0.01,
             guinier_scale=1.7,
Index: sasmodels/models/guinier.py
===================================================================
--- sasmodels/models/guinier.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
+++ sasmodels/models/guinier.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -51,6 +51,8 @@
 def random():
     import numpy as np
-    scale = 10**np.random.uniform(1, 5)
-    # Note: compare.py has rg cutoff for guinier, so use that
+    scale = 10**np.random.uniform(1, 4)
+    # Note: compare.py has Rg cutoff of 1e-30 at q=1 for guinier, so use that
+    # log_10 Ae^(-(q Rg)^2/3) = log_10(A) - (q Rg)^2/ (3 ln 10) > -30
+    #   => log_10(A) > Rg^2/(3 ln 10) - 30
     q_max = 1.0
     rg_max = np.sqrt(90*np.log(10) + 3*np.log(scale))/q_max
Index: sasmodels/models/guinier_porod.py
===================================================================
--- sasmodels/models/guinier_porod.py	(revision cdcebf1033e85ca2bb11488372c82a9050ec942a)
+++ sasmodels/models/guinier_porod.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -4,5 +4,5 @@
 and dimensionality of scattering objects, including asymmetric objects
 such as rods or platelets, and shapes intermediate between spheres
-and rods or between rods and platelets, and overcomes some of the 
+and rods or between rods and platelets, and overcomes some of the
 deficiencies of the (Beaucage) Unified_Power_Rg model (see Hammouda, 2010).
 
@@ -77,5 +77,5 @@
                         scale = Guinier Scale
                         s = Dimension Variable
-                        Rg = Radius of Gyration [A] 
+                        Rg = Radius of Gyration [A]
                         porod_exp = Porod Exponent
                         background  = Background [1/cm]"""
@@ -114,4 +114,17 @@
 Iq.vectorized = True # Iq accepts an array of q values
 
+def random():
+    import numpy as np
+    rg = 10**np.random.uniform(1, 5)
+    s = np.random.uniform(0, 3)
+    porod_exp = s + np.random.uniform(0, 3)
+    pars = dict(
+        #scale=1, background=0,
+        rg=rg,
+        s=s,
+        porod_exp=porod_exp,
+    )
+    return pars
+
 demo = dict(scale=1.5, background=0.5, rg=60, s=1.0, porod_exp=3.0)
 
Index: sasmodels/models/line.py
===================================================================
--- sasmodels/models/line.py	(revision 404ebbd7d77c5752ee99a8b39bf8172dd0c0e551)
+++ sasmodels/models/line.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -70,9 +70,11 @@
 def random():
     import numpy as np
-    slope = 1e5*np.random.uniform(-1, 1)
+    scale = 10**np.random.uniform(0, 3)
+    slope = np.random.uniform(-1, 1)*1e2
     offset = 1e-5 + (0 if slope > 0 else -slope)
-    intercept = 10**np.random.uniform(0, 5) + offset
+    intercept = 10**np.random.uniform(0, 1) + offset
     pars = dict(
-        scale=1, background=0,
+        #background=0,
+        scale=scale,
         slope=slope,
         intercept=intercept,
Index: sasmodels/models/spinodal.py
===================================================================
--- sasmodels/models/spinodal.py	(revision bba9361c45987d76f3673d7d0a4eb16d28a8b701)
+++ sasmodels/models/spinodal.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -3,14 +3,14 @@
 ----------
 
-This model calculates the SAS signal of a phase separating solution under spinodal decomposition. 
+This model calculates the SAS signal of a phase separating solution under spinodal decomposition.
 The scattering intensity $I(q)$ is calculated as
 
-.. math:: 
+.. math::
     I(q) = I_{max}\frac{(1+\gamma/2)x^2}{\gamma/2+x^{2+\gamma}}+B
 
-where $x=q/q_0$ and $B$ is a flat background. The characteristic structure length 
-scales with the correlation peak at $q_0$. The exponent $\gamma$ is equal to 
+where $x=q/q_0$ and $B$ is a flat background. The characteristic structure length
+scales with the correlation peak at $q_0$. The exponent $\gamma$ is equal to
 $d+1$ with d the dimensionality of the off-critical concentration mixtures.
-A transition to $\gamma=2d$ is seen near the percolation threshold into the 
+A transition to $\gamma=2d$ is seen near the percolation threshold into the
 critical concentration regime.
 
@@ -18,5 +18,5 @@
 ----------
 
-H. Furukawa. Dynamics-scaling theory for phase-separating unmixing mixtures: 
+H. Furukawa. Dynamics-scaling theory for phase-separating unmixing mixtures:
 Growth rates of droplets and scaling properties of autocorrelation functions. Physica A 123,497 (1984).
 
@@ -25,6 +25,6 @@
 
 * **Author:** Dirk Honecker **Date:** Oct 7, 2016
-* **Last Modified by:** 
-* **Last Reviewed by:** 
+* **Last Modified by:**
+* **Last Reviewed by:**
 """
 
@@ -46,6 +46,5 @@
 # pylint: disable=bad-whitespace, line-too-long
 #             ["name", "units", default, [lower, upper], "type", "description"],
-parameters = [["scale",    "",      1.0, [-inf, inf], "", "Scale factor"],
-              ["gamma",      "",    3.0, [-inf, inf], "", "Exponent"],
+parameters = [["gamma",      "",    3.0, [-inf, inf], "", "Exponent"],
               ["q_0",  "1/Ang",     0.1, [-inf, inf], "", "Correlation peak position"]
              ]
@@ -53,21 +52,28 @@
 
 def Iq(q,
-       scale=1.0,
        gamma=3.0,
        q_0=0.1):
     """
     :param q:              Input q-value
-    :param scale:          Scale factor
     :param gamma:          Exponent
     :param q_0:            Correlation peak position
     :return:               Calculated intensity
     """
-    
+
     with errstate(divide='ignore'):
         x = q/q_0
-        inten = scale * ((1 + gamma / 2) * x ** 2) / (gamma / 2 + x ** (2 + gamma)) 
+        inten = ((1 + gamma / 2) * x ** 2) / (gamma / 2 + x ** (2 + gamma))
     return inten
 Iq.vectorized = True  # Iq accepts an array of q values
 
+def random():
+    import numpy as np
+    pars = dict(
+        scale=10**np.random.uniform(1, 3),
+        gamma=np.random.uniform(0, 6),
+        q_0=10**np.random.uniform(-3, -1),
+    )
+    return pars
+
 demo = dict(scale=1, background=0,
             gamma=1, q_0=0.1)
Index: sasmodels/models/star_polymer.py
===================================================================
--- sasmodels/models/star_polymer.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/models/star_polymer.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -59,7 +59,13 @@
 source = ["star_polymer.c"]
 
-demo = dict(scale=1, background=0,
-            rg_squared=100.0,
-            arms=3.0)
+def random():
+    import numpy as np
+    pars = dict(
+        #background=0,
+        scale=10**np.random.uniform(1, 4),
+        rg_squared=10**np.random.uniform(1, 8),
+        arms=np.random.uniform(1, 6),
+    )
+    return pars
 
 tests = [[{'rg_squared': 2.0,
Index: sasmodels/models/surface_fractal.py
===================================================================
--- sasmodels/models/surface_fractal.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/surface_fractal.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -74,6 +74,17 @@
 source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "surface_fractal.c"]
 
-demo = dict(scale=1, background=1e-5,
-            radius=10, fractal_dim_surf=2.0, cutoff_length=500)
+def random():
+    import numpy as np
+    radius = 10**np.random.uniform(1, 4)
+    fractal_dim_surf = np.random.uniform(1, 3-1e-6)
+    cutoff_length = 1e6  # Sets the low q limit; keep it big for sim
+    pars = dict(
+        #background=0,
+        scale=1,
+        radius=radius,
+        fractal_dim_surf=fractal_dim_surf,
+        cutoff_length=cutoff_length,
+    )
+    return pars
 
 tests = [
Index: sasmodels/models/teubner_strey.py
===================================================================
--- sasmodels/models/teubner_strey.py	(revision 8393c74ed553108e237a0c93b1e6c65ff328e00d)
+++ sasmodels/models/teubner_strey.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -102,6 +102,21 @@
 Iq.vectorized = True  # Iq accepts an array of q values
 
+def random():
+    import numpy as np
+    d = 10**np.random.uniform(1, 4)
+    xi = 10**np.random.uniform(-0.3, 2)*d
+    pars = dict(
+        #background=0,
+        scale=100,
+        volfraction_a=10**np.random.uniform(-3, 0),
+        sld_a=np.random.uniform(-0.5, 12),
+        sld_b=np.random.uniform(-0.5, 12),
+        d=d,
+        xi=xi,
+    )
+    return pars
+
 demo = dict(scale=1, background=0, volfraction_a=0.5,
-                     sld_a=0.3, sld_b=6.3,
-                     d=100.0, xi=30.0)
+            sld_a=0.3, sld_b=6.3,
+            d=100.0, xi=30.0)
 tests = [[{}, 0.06, 41.5918888453]]
Index: sasmodels/models/two_lorentzian.py
===================================================================
--- sasmodels/models/two_lorentzian.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
+++ sasmodels/models/two_lorentzian.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -93,4 +93,22 @@
 Iq.vectorized = True  # Iq accepts an array of q values
 
+def random():
+    import numpy as np
+    scale = 10**np.random.uniform(0, 4, 2)
+    length = 10**np.random.uniform(1, 4, 2)
+    expon = np.random.uniform(1, 6, 2)
+
+    pars = dict(
+        #background=0,
+        scale=1, # scale provided in model
+        lorentz_scale_1=scale[0],
+        lorentz_length_1=length[0],
+        lorentz_exp_1=expon[0],
+        lorentz_scale_2=scale[1],
+        lorentz_length_2=length[1],
+        lorentz_exp_2=expon[1],
+    )
+    return pars
+
 
 demo = dict(scale=1, background=0.1,
Index: sasmodels/models/two_power_law.py
===================================================================
--- sasmodels/models/two_power_law.py	(revision bb4b50990533a15ceda4443988db938e13a28a2b)
+++ sasmodels/models/two_power_law.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -98,4 +98,19 @@
 Iq.vectorized = True  # Iq accepts an array of q values
 
+def random():
+    import numpy as np
+    coefficient_1 = 1
+    crossover = 10**np.random.uniform(-3, -1)
+    power_1 = np.random.uniform(1, 6)
+    power_2 = np.random.uniform(1, 6)
+    pars = dict(
+        scale=1, #background=0,
+        coefficient_1=coefficient_1,
+        crossover=crossover,
+        power_1=power_1,
+        power_2=power_2,
+    )
+    return pars
+
 demo = dict(scale=1, background=0.0,
             coefficent_1=1.0,
Index: sasmodels/models/unified_power_Rg.py
===================================================================
--- sasmodels/models/unified_power_Rg.py	(revision 66ca2a6cf404d09b601f9022e85222304a98dc8e)
+++ sasmodels/models/unified_power_Rg.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -3,7 +3,7 @@
 ----------
 
-This model employs the empirical multiple level unified Exponential/Power-law 
-fit method developed by Beaucage. Four functions are included so that 1, 2, 3, 
-or 4 levels can be used. In addition a 0 level has been added which simply 
+This model employs the empirical multiple level unified Exponential/Power-law
+fit method developed by Beaucage. Four functions are included so that 1, 2, 3,
+or 4 levels can be used. In addition a 0 level has been added which simply
 calculates
 
@@ -16,8 +16,8 @@
 (Debye equation), ellipsoidal particles, etc.
 
-The model works best for mass fractal systems characterized by Porod exponents 
-between 5/3 and 3. It should not be used for surface fractal systems. Hammouda 
-(2010) has pointed out a deficiency in the way this model handles the 
-transitioning between the Guinier and Porod regimes and which can create 
+The model works best for mass fractal systems characterized by Porod exponents
+between 5/3 and 3. It should not be used for surface fractal systems. Hammouda
+(2010) has pointed out a deficiency in the way this model handles the
+transitioning between the Guinier and Porod regimes and which can create
 artefacts that appear as kinks in the fitted model function.
 
@@ -88,5 +88,5 @@
 # pylint: disable=bad-whitespace, line-too-long
 parameters = [
-    ["level",     "",     1,      [0, 6], "", "Level number"],
+    ["level",     "",     1,      [1, 6], "", "Level number"],
     ["rg[level]", "Ang",  15.8,   [0, inf], "", "Radius of gyration"],
     ["power[level]", "",  4,      [-inf, inf], "", "Power"],
@@ -118,4 +118,25 @@
 Iq.vectorized = True
 
+def random():
+    import numpy as np
+    from scipy.special import gamma
+    level = np.minimum(np.random.poisson(0.5) + 1, 6)
+    n = level
+    power = np.random.uniform(1.6, 3, n)
+    rg = 10**np.random.uniform(1, 5, n)
+    G = np.random.uniform(0.1, 10, n)**2 * 10**np.random.uniform(0.3, 3, n)
+    B = G * power / rg**power * gamma(power/2)
+    scale = 10**np.random.uniform(1, 4)
+    pars = dict(
+        #background=0,
+        scale=scale,
+        level=level,
+    )
+    pars.update(("power%d"%(k+1), v) for k, v in enumerate(power))
+    pars.update(("rg%d"%(k+1), v) for k, v in enumerate(rg))
+    pars.update(("B%d"%(k+1), v) for k, v in enumerate(B))
+    pars.update(("G%d"%(k+1), v) for k, v in enumerate(G))
+    return pars
+
 demo = dict(
     level=2,
Index: sasmodels/models/vesicle.py
===================================================================
--- sasmodels/models/vesicle.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/vesicle.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
@@ -120,4 +120,18 @@
     return whole, whole - core
 
+def random():
+    import numpy as np
+    total_radius = 10**np.random.uniform(1.3, 5)
+    radius = total_radius * np.random.uniform(0, 1)
+    thickness = total_radius - radius
+    volfraction = 10**np.random.uniform(-3, -1)
+    pars = dict(
+        #background=0,
+        scale=1,  # volfraction is part of the model, so scale=1
+        radius=radius,
+        thickness=thickness,
+        volfraction=volfraction,
+    )
+    return pars
 
 # parameters for demo
