Index: doc/guide/magnetism/magnetism.rst
===================================================================
--- doc/guide/magnetism/magnetism.rst	(revision befe905ba1f75b9bc31808657037e4ea891bdf28)
+++ doc/guide/magnetism/magnetism.rst	(revision df87acf3bc0a5479966307d713d69f552ccaf908)
@@ -89,14 +89,14 @@
 
 ===========   ================================================================
- M0:sld       $D_M M_0$
- mtheta:sld   $\theta_M$
- mphi:sld     $\phi_M$
- up:angle     $\theta_\mathrm{up}$
- up:frac_i    $u_i$ = (spin up)/(spin up + spin down) *before* the sample
- up:frac_f    $u_f$ = (spin up)/(spin up + spin down) *after* the sample
+ sld_M0       $D_M M_0$
+ sld_mtheta   $\theta_M$
+ sld_mphi     $\phi_M$
+ up_frac_i    $u_i$ = (spin up)/(spin up + spin down) *before* the sample
+ up_frac_f    $u_f$ = (spin up)/(spin up + spin down) *after* the sample
+ up_angle     $\theta_\mathrm{up}$
 ===========   ================================================================
 
 .. note::
-    The values of the 'up:frac_i' and 'up:frac_f' must be in the range 0 to 1.
+    The values of the 'up_frac_i' and 'up_frac_f' must be in the range 0 to 1.
 
 *Document History*
Index: doc/guide/plugin.rst
===================================================================
--- doc/guide/plugin.rst	(revision f796469ab131e5fc2424ecdfcd201846be55c8e1)
+++ doc/guide/plugin.rst	(revision 57c609bdedf506354431a1f62e27c73c368bb043)
@@ -423,4 +423,20 @@
 calculations, but instead rely on numerical integration to compute the
 appropriately smeared pattern.
+
+Each .py file also contains a function::
+
+	def random():
+	...
+
+This function provides a model-specific random parameter set which shows model
+features in the USANS to SANS range.  For example, core-shell sphere sets the
+outer radius of the sphere logarithmically in `[20, 20,000]`, which sets the Q
+value for the transition from flat to falling.  It then uses a beta distribution
+to set the percentage of the shape which is shell, giving a preference for very
+thin or very thick shells (but never 0% or 100%).  Using `-sets=10` in sascomp
+should show a reasonable variety of curves over the default sascomp q range.
+The parameter set is returned as a dictionary of `{parameter: value, ...}`.
+Any model parameters not included in the dictionary will default according to
+the code in the `_randomize_one()` function from sasmodels/compare.py.
 
 Python Models
@@ -685,6 +701,6 @@
     erf, erfc, tgamma, lgamma:  **do not use**
         Special functions that should be part of the standard, but are missing
-        or inaccurate on some platforms. Use sas_erf, sas_erfc and sas_gamma
-        instead (see below). Note: lgamma(x) has not yet been tested.
+        or inaccurate on some platforms. Use sas_erf, sas_erfc, sas_gamma
+        and sas_lgamma instead (see below).
 
 Some non-standard constants and functions are also provided:
@@ -753,9 +769,27 @@
         Gamma function sas_gamma\ $(x) = \Gamma(x)$.
 
-        The standard math function, tgamma(x) is unstable for $x < 1$
+        The standard math function, tgamma(x), is unstable for $x < 1$
         on some platforms.
 
         :code:`source = ["lib/sas_gamma.c", ...]`
         (`sas_gamma.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gamma.c>`_)
+
+    sas_gammaln(x):
+        log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$.
+
+        The standard math function, lgamma(x), is incorrect for single
+        precision on some platforms.
+
+        :code:`source = ["lib/sas_gammainc.c", ...]`
+        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_)
+
+    sas_gammainc(a, x), sas_gammaincc(a, x):
+        Incomplete gamma function
+        sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$
+        and complementary incomplete gamma function
+        sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$
+
+        :code:`source = ["lib/sas_gammainc.c", ...]`
+        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_)
 
     sas_erf(x), sas_erfc(x):
@@ -795,4 +829,6 @@
         If $n$ = 0 or 1, it uses sas_J0($x$) or sas_J1($x$), respectively.
 
+        Warning: JN(n,x) can be very inaccurate (0.1%) for x not in [0.1, 100].
+
         The standard math function jn(n, x) is not available on all platforms.
 
@@ -803,7 +839,9 @@
         Sine integral Si\ $(x) = \int_0^x \tfrac{\sin t}{t}\,dt$.
 
+        Warning: Si(x) can be very inaccurate (0.1%) for x in [0.1, 100].
+
         This function uses Taylor series for small and large arguments:
 
-        For large arguments,
+        For large arguments use the following Taylor series,
 
         .. math::
@@ -813,5 +851,5 @@
              - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right)
 
-        For small arguments,
+        For small arguments ,
 
         .. math::
Index: explore/precision.py
===================================================================
--- explore/precision.py	(revision 2a7e20e9c4ba72da5c6088e743a64f3345ee554d)
+++ explore/precision.py	(revision fba9ca0239129ed28284a0b435ba2a9524913f36)
@@ -95,4 +95,7 @@
             neg:    [-100,100]
 
+        For arbitrary range use "start:stop:steps:scale" where scale is
+        one of log, lin, or linear.
+
         *diff* is "relative", "absolute" or "none"
 
@@ -102,14 +105,21 @@
         linear = not xrange.startswith("log")
         if xrange == "zoom":
-            lin_min, lin_max, lin_steps = 1000, 1010, 2000
+            start, stop, steps = 1000, 1010, 2000
         elif xrange == "neg":
-            lin_min, lin_max, lin_steps = -100.1, 100.1, 2000
+            start, stop, steps = -100.1, 100.1, 2000
         elif xrange == "linear":
-            lin_min, lin_max, lin_steps = 1, 1000, 2000
-            lin_min, lin_max, lin_steps = 0.001, 2, 2000
+            start, stop, steps = 1, 1000, 2000
+            start, stop, steps = 0.001, 2, 2000
         elif xrange == "log":
-            log_min, log_max, log_steps = -3, 5, 400
+            start, stop, steps = -3, 5, 400
         elif xrange == "logq":
-            log_min, log_max, log_steps = -4, 1, 400
+            start, stop, steps = -4, 1, 400
+        elif ':' in xrange:
+            parts = xrange.split(':')
+            linear = parts[3] != "log" if len(parts) == 4 else True
+            steps = int(parts[2]) if len(parts) > 2 else 400
+            start = float(parts[0])
+            stop = float(parts[1])
+
         else:
             raise ValueError("unknown range "+xrange)
@@ -121,17 +131,17 @@
             # value to x in the given precision.
             if linear:
-                lin_min = max(lin_min, self.limits[0])
-                lin_max = min(lin_max, self.limits[1])
-                qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='single')
-                #qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='double')
+                start = max(start, self.limits[0])
+                stop = min(stop, self.limits[1])
+                qrf = np.linspace(start, stop, steps, dtype='single')
+                #qrf = np.linspace(start, stop, steps, dtype='double')
                 qr = [mp.mpf(float(v)) for v in qrf]
-                #qr = mp.linspace(lin_min, lin_max, lin_steps)
+                #qr = mp.linspace(start, stop, steps)
             else:
-                log_min = np.log10(max(10**log_min, self.limits[0]))
-                log_max = np.log10(min(10**log_max, self.limits[1]))
-                qrf = np.logspace(log_min, log_max, log_steps, dtype='single')
-                #qrf = np.logspace(log_min, log_max, log_steps, dtype='double')
+                start = np.log10(max(10**start, self.limits[0]))
+                stop = np.log10(min(10**stop, self.limits[1]))
+                qrf = np.logspace(start, stop, steps, dtype='single')
+                #qrf = np.logspace(start, stop, steps, dtype='double')
                 qr = [mp.mpf(float(v)) for v in qrf]
-                #qr = [10**v for v in mp.linspace(log_min, log_max, log_steps)]
+                #qr = [10**v for v in mp.linspace(start, stop, steps)]
 
         target = self.call_mpmath(qr, bits=500)
@@ -176,5 +186,5 @@
     """
     if diff == "relative":
-        err = np.array([abs((t-a)/t) for t, a in zip(target, actual)], 'd')
+        err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd')
         #err = np.clip(err, 0, 1)
         pylab.loglog(x, err, '-', label=label)
@@ -197,4 +207,21 @@
     return model_info
 
+# Hack to allow second parameter A in two parameter functions
+A = 1
+def parse_extra_pars():
+    global A
+
+    A_str = str(A)
+    pop = []
+    for k, v in enumerate(sys.argv[1:]):
+        if v.startswith("A="):
+            A_str = v[2:]
+            pop.append(k+1)
+    if pop:
+        sys.argv = [v for k, v in enumerate(sys.argv) if k not in pop]
+        A = float(A_str)
+
+parse_extra_pars()
+
 
 # =============== FUNCTION DEFINITIONS ================
@@ -297,4 +324,23 @@
     ocl_function=make_ocl("return sas_gamma(q);", "sas_gamma", ["lib/sas_gamma.c"]),
     limits=(-3.1, 10),
+)
+add_function(
+    name="gammaln(x)",
+    mp_function=mp.loggamma,
+    np_function=scipy.special.gammaln,
+    ocl_function=make_ocl("return sas_gammaln(q);", "sas_gammaln", ["lib/sas_gammainc.c"]),
+    #ocl_function=make_ocl("return lgamma(q);", "sas_gammaln"),
+)
+add_function(
+    name="gammainc(x)",
+    mp_function=lambda x, a=A: mp.gammainc(a, a=0, b=x)/mp.gamma(a),
+    np_function=lambda x, a=A: scipy.special.gammainc(a, x),
+    ocl_function=make_ocl("return sas_gammainc(%.15g,q);"%A, "sas_gammainc", ["lib/sas_gammainc.c"]),
+)
+add_function(
+    name="gammaincc(x)",
+    mp_function=lambda x, a=A: mp.gammainc(a, a=x, b=mp.inf)/mp.gamma(a),
+    np_function=lambda x, a=A: scipy.special.gammaincc(a, x),
+    ocl_function=make_ocl("return sas_gammaincc(%.15g,q);"%A, "sas_gammaincc", ["lib/sas_gammainc.c"]),
 )
 add_function(
@@ -463,6 +509,6 @@
 lanczos_gamma = """\
     const double coeff[] = {
-            76.18009172947146,     -86.50532032941677,
-            24.01409824083091,     -1.231739572450155,
+            76.18009172947146, -86.50532032941677,
+            24.01409824083091, -1.231739572450155,
             0.1208650973866179e-2,-0.5395239384953e-5
             };
@@ -475,5 +521,5 @@
 """
 add_function(
-    name="log gamma(x)",
+    name="loggamma(x)",
     mp_function=mp.loggamma,
     np_function=scipy.special.gammaln,
@@ -599,5 +645,5 @@
 
 ALL_FUNCTIONS = set(FUNCTIONS.keys())
-ALL_FUNCTIONS.discard("loggamma")  # OCL version not ready yet
+ALL_FUNCTIONS.discard("loggamma")  # use cephes-based gammaln instead
 ALL_FUNCTIONS.discard("3j1/x:taylor")
 ALL_FUNCTIONS.discard("3j1/x:trig")
@@ -615,10 +661,15 @@
     -r indicates that the relative error should be plotted (default),
     -x<range> indicates the steps in x, where <range> is one of the following
-      log indicates log stepping in [10^-3, 10^5] (default)
-      logq indicates log stepping in [10^-4, 10^1]
-      linear indicates linear stepping in [1, 1000]
-      zoom indicates linear stepping in [1000, 1010]
-      neg indicates linear stepping in [-100.1, 100.1]
-and name is "all" or one of:
+        log indicates log stepping in [10^-3, 10^5] (default)
+        logq indicates log stepping in [10^-4, 10^1]
+        linear indicates linear stepping in [1, 1000]
+        zoom indicates linear stepping in [1000, 1010]
+        neg indicates linear stepping in [-100.1, 100.1]
+        start:stop:n[:stepping] indicates an n-step plot in [start, stop]
+            or [10^start, 10^stop] if stepping is "log" (default n=400)
+Some functions (notably gammainc/gammaincc) have an additional parameter A
+which can be set from the command line as A=value.  Default is A=1.
+
+Name is one of:
     """+names)
     sys.exit(1)
Index: sasmodels/__init__.py
===================================================================
--- sasmodels/__init__.py	(revision e65c3ba5944b729f651cf8747f7b64bed2cbb826)
+++ sasmodels/__init__.py	(revision a1ec9081782bea4f4d869eec0518e8eeebd50f18)
@@ -14,5 +14,5 @@
 defining new models.
 """
-__version__ = "0.97"
+__version__ = "0.98"
 
 def data_files():
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision bd7630de7da8bd67aa24b8fc036b2dcc1bc72a16)
+++ sasmodels/compare.py	(revision 610ef23645ef28d7c1e8c87aab1e03c8cd123a93)
@@ -368,5 +368,5 @@
 
     # Limit magnetic SLDs to a smaller range, from zero to iron=5/A^2
-    if par.name.startswith('M0:'):
+    if par.name.endswith('_M0'):
         return np.random.uniform(0, 5)
 
@@ -538,5 +538,5 @@
     magnetic_pars = []
     for p in parameters.user_parameters(pars, is2d):
-        if any(p.id.startswith(x) for x in ('M0:', 'mtheta:', 'mphi:')):
+        if any(p.id.endswith(x) for x in ('_M0', '_mtheta', '_mphi')):
             continue
         if p.id.startswith('up:'):
@@ -550,7 +550,7 @@
             pdtype=pars.get(p.id+"_pd_type", 'gaussian'),
             relative_pd=p.relative_pd,
-            M0=pars.get('M0:'+p.id, 0.),
-            mphi=pars.get('mphi:'+p.id, 0.),
-            mtheta=pars.get('mtheta:'+p.id, 0.),
+            M0=pars.get(p.id+'_M0', 0.),
+            mphi=pars.get(p.id+'_mphi', 0.),
+            mtheta=pars.get(p.id+'_mtheta', 0.),
         )
         lines.append(_format_par(p.name, **fields))
@@ -618,5 +618,5 @@
     if suppress:
         for p in pars:
-            if p.startswith("M0:"):
+            if p.endswith("_M0"):
                 pars[p] = 0
     else:
@@ -624,5 +624,5 @@
         first_mag = None
         for p in pars:
-            if p.startswith("M0:"):
+            if p.endswith("_M0"):
                 any_mag |= (pars[p] != 0)
                 if first_mag is None:
Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision 78f830855309bd88b3236276481ada35aeafea06)
+++ sasmodels/convert.py	(revision 610ef23645ef28d7c1e8c87aab1e03c8cd123a93)
@@ -167,4 +167,5 @@
     if version < (4, 2, 0):
         oldpars = _hand_convert_4_1_to_4_2(name, oldpars)
+        oldpars = _rename_magnetic_pars(oldpars)
     return oldpars
 
@@ -174,4 +175,20 @@
         oldpars['lattice_distortion'] = oldpars.pop('d_factor')
     return oldpars
+
+def _rename_magnetic_pars(pars):
+    """
+    Change from M0:par to par_M0, etc.
+    """
+    keys = list(pars.items())
+    for k in keys:
+        if k.startswith('M0:'):
+            pars[k[3:]+'_M0'] = pars.pop(k)
+        elif k.startswith('mtheta:'):
+            pars[k[7:]+'_mtheta'] = pars.pop(k)
+        elif k.startswith('mphi:'):
+            pars[k[5:]+'_mphi'] = pars.pop(k)
+        elif k.startswith('up:'):
+            pars['up_'+k[3:]] = pars.pop(k)
+    return pars
 
 def _hand_convert_3_1_2_to_4_1(name, oldpars):
Index: sasmodels/custom/__init__.py
===================================================================
--- sasmodels/custom/__init__.py	(revision 0f48f1e9570fb194a19f0d8faf8c23c879a61d8b)
+++ sasmodels/custom/__init__.py	(revision d321747456846675ebcfb426f57c90d80040210d)
@@ -12,5 +12,5 @@
 import sys
 import os
-from os.path import basename, splitext
+from os.path import basename, splitext, join as joinpath, exists, dirname
 
 try:
@@ -18,4 +18,5 @@
     from importlib.util import spec_from_file_location, module_from_spec  # type: ignore
     def load_module_from_path(fullname, path):
+        # type: (str, str) -> "module"
         """load module from *path* as *fullname*"""
         spec = spec_from_file_location(fullname, os.path.expanduser(path))
@@ -27,4 +28,5 @@
     import imp
     def load_module_from_path(fullname, path):
+        # type: (str, str) -> "module"
         """load module from *path* as *fullname*"""
         # Clear out old definitions, if any
@@ -37,10 +39,83 @@
         return module
 
+_MODULE_CACHE = {} # type: Dict[str, Tuple("module", int)]
+_MODULE_DEPENDS = {} # type: Dict[str, List[str]]
+_MODULE_DEPENDS_STACK = [] # type: List[str]
 def load_custom_kernel_module(path):
+    # type: str -> "module"
     """load SAS kernel from *path* as *sasmodels.custom.modelname*"""
     # Pull off the last .ext if it exists; there may be others
     name = basename(splitext(path)[0])
-    # Placing the model in the 'sasmodels.custom' name space.
-    kernel_module = load_module_from_path('sasmodels.custom.'+name,
-                                          os.path.expanduser(path))
-    return kernel_module
+    path = os.path.expanduser(path)
+
+    # Reload module if necessary.
+    if need_reload(path):
+        # Assume the module file is the only dependency
+        _MODULE_DEPENDS[path] = set([path])
+
+        # Load the module while pushing it onto the dependency stack.  If
+        # this triggers any submodules, then they will add their dependencies
+        # to this module as the "working_on" parent.  Pop the stack when the
+        # module is loaded.
+        _MODULE_DEPENDS_STACK.append(path)
+        module = load_module_from_path('sasmodels.custom.'+name, path)
+        _MODULE_DEPENDS_STACK.pop()
+
+        # Include external C code in the dependencies.  We are looking
+        # for module.source and assuming that it is a list of C source files
+        # relative to the module itself.  Any files that do not exist,
+        # such as those in the standard libraries, will be ignored.
+        # TODO: look in builtin module path for standard c sources
+        # TODO: share code with generate.model_sources
+        c_sources = getattr(module, 'source', None)
+        if isinstance(c_sources, (list, tuple)):
+            _MODULE_DEPENDS[path].update(_find_sources(path, c_sources))
+
+        # Cache the module, and tag it with the newest timestamp
+        timestamp = max(os.path.getmtime(f) for f in _MODULE_DEPENDS[path])
+        _MODULE_CACHE[path] = module, timestamp
+
+        #print("loading", os.path.basename(path), _MODULE_CACHE[path][1],
+        #    [os.path.basename(p) for p in _MODULE_DEPENDS[path]])
+
+    # Add path and all its dependence to the parent module, if there is one.
+    if _MODULE_DEPENDS_STACK:
+        working_on = _MODULE_DEPENDS_STACK[-1]
+        _MODULE_DEPENDS[working_on].update(_MODULE_DEPENDS[path])
+
+    return _MODULE_CACHE[path][0]
+
+def need_reload(path):
+    # type: str -> bool
+    """
+    Return True if any path dependencies have a timestamp newer than the time
+    when the path was most recently loaded.
+    """
+    # TODO: fails if a dependency has a modification time in the future
+    # If the newest dependency has a time stamp in the future, then this
+    # will be recorded as the cached time.  When a second dependency
+    # is updated to the current time stamp, it will still be considered
+    # older than the current build and the reload will not be triggered.
+    # Could instead treat all future times as 0 here and in the code above
+    # which records the newest timestamp.  This will force a reload when
+    # the future time is reached, but other than that should perform
+    # correctly.  Probably not worth the extra code...
+    _, cache_time = _MODULE_CACHE.get(path, (None, -1))
+    depends = _MODULE_DEPENDS.get(path, [path])
+    #print("reload", any(cache_time < os.path.getmtime(p) for p in depends))
+    #for f in depends: print(">>>  ", f, os.path.getmtime(f))
+    return any(cache_time < os.path.getmtime(p) for p in depends)
+
+def _find_sources(path, source_list):
+    # type: (str, List[str]) -> List[str]
+    """
+    Return a list of the sources relative to base file; ignore any that
+    are not found.
+    """
+    root = dirname(path)
+    found = []
+    for source_name in source_list:
+        source_path = joinpath(root, source_name)
+        if exists(source_path):
+            found.append(source_path)
+    return found
Index: sasmodels/kernelpy.py
===================================================================
--- sasmodels/kernelpy.py	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
+++ sasmodels/kernelpy.py	(revision 91bd550ddd44f57ea241f0f2e3f7bb1d1ab4f5ee)
@@ -37,5 +37,5 @@
         self.info = model_info
         self.dtype = np.dtype('d')
-        logger.info("load python model " + self.info.name)
+        logger.info("make python model " + self.info.name)
 
     def make_kernel(self, q_vectors):
Index: sasmodels/model_test.py
===================================================================
--- sasmodels/model_test.py	(revision 012cd342ba6d27483147e3401f01fb7b93ea589c)
+++ sasmodels/model_test.py	(revision 12eec1eafef77dae90813c12e016e43de8ffe1c4)
@@ -47,4 +47,5 @@
 import sys
 import unittest
+import traceback
 
 try:
@@ -74,5 +75,4 @@
 # pylint: enable=unused-import
 
-
 def make_suite(loaders, models):
     # type: (List[str], List[str]) -> unittest.TestSuite
@@ -86,5 +86,4 @@
     *models* is the list of models to test, or *["all"]* to test all models.
     """
-    ModelTestCase = _hide_model_case_from_nose()
     suite = unittest.TestSuite()
 
@@ -95,70 +94,74 @@
         skip = []
     for model_name in models:
-        if model_name in skip:
-            continue
-        model_info = load_model_info(model_name)
-
-        #print('------')
-        #print('found tests in', model_name)
-        #print('------')
-
-        # if ispy then use the dll loader to call pykernel
-        # don't try to call cl kernel since it will not be
-        # available in some environmentes.
-        is_py = callable(model_info.Iq)
-
-        # Some OpenCL drivers seem to be flaky, and are not producing the
-        # expected result.  Since we don't have known test values yet for
-        # all of our models, we are instead going to compare the results
-        # for the 'smoke test' (that is, evaluation at q=0.1 for the default
-        # parameters just to see that the model runs to completion) between
-        # the OpenCL and the DLL.  To do this, we define a 'stash' which is
-        # shared between OpenCL and DLL tests.  This is just a list.  If the
-        # list is empty (which it will be when DLL runs, if the DLL runs
-        # first), then the results are appended to the list.  If the list
-        # is not empty (which it will be when OpenCL runs second), the results
-        # are compared to the results stored in the first element of the list.
-        # This is a horrible stateful hack which only makes sense because the
-        # test suite is thrown away after being run once.
-        stash = []
-
-        if is_py:  # kernel implemented in python
-            test_name = "%s-python"%model_name
-            test_method_name = "test_%s_python" % model_info.id
+        if model_name not in skip:
+            model_info = load_model_info(model_name)
+            _add_model_to_suite(loaders, suite, model_info)
+
+    return suite
+
+def _add_model_to_suite(loaders, suite, model_info):
+    ModelTestCase = _hide_model_case_from_nose()
+
+    #print('------')
+    #print('found tests in', model_name)
+    #print('------')
+
+    # if ispy then use the dll loader to call pykernel
+    # don't try to call cl kernel since it will not be
+    # available in some environmentes.
+    is_py = callable(model_info.Iq)
+
+    # Some OpenCL drivers seem to be flaky, and are not producing the
+    # expected result.  Since we don't have known test values yet for
+    # all of our models, we are instead going to compare the results
+    # for the 'smoke test' (that is, evaluation at q=0.1 for the default
+    # parameters just to see that the model runs to completion) between
+    # the OpenCL and the DLL.  To do this, we define a 'stash' which is
+    # shared between OpenCL and DLL tests.  This is just a list.  If the
+    # list is empty (which it will be when DLL runs, if the DLL runs
+    # first), then the results are appended to the list.  If the list
+    # is not empty (which it will be when OpenCL runs second), the results
+    # are compared to the results stored in the first element of the list.
+    # This is a horrible stateful hack which only makes sense because the
+    # test suite is thrown away after being run once.
+    stash = []
+
+    if is_py:  # kernel implemented in python
+        test_name = "%s-python"%model_info.name
+        test_method_name = "test_%s_python" % model_info.id
+        test = ModelTestCase(test_name, model_info,
+                                test_method_name,
+                                platform="dll",  # so that
+                                dtype="double",
+                                stash=stash)
+        suite.addTest(test)
+    else:   # kernel implemented in C
+
+        # test using dll if desired
+        if 'dll' in loaders or not use_opencl():
+            test_name = "%s-dll"%model_info.name
+            test_method_name = "test_%s_dll" % model_info.id
             test = ModelTestCase(test_name, model_info,
-                                 test_method_name,
-                                 platform="dll",  # so that
-                                 dtype="double",
-                                 stash=stash)
+                                    test_method_name,
+                                    platform="dll",
+                                    dtype="double",
+                                    stash=stash)
             suite.addTest(test)
-        else:   # kernel implemented in C
-
-            # test using dll if desired
-            if 'dll' in loaders or not use_opencl():
-                test_name = "%s-dll"%model_name
-                test_method_name = "test_%s_dll" % model_info.id
-                test = ModelTestCase(test_name, model_info,
-                                     test_method_name,
-                                     platform="dll",
-                                     dtype="double",
-                                     stash=stash)
-                suite.addTest(test)
-
-            # test using opencl if desired and available
-            if 'opencl' in loaders and use_opencl():
-                test_name = "%s-opencl"%model_name
-                test_method_name = "test_%s_opencl" % model_info.id
-                # Using dtype=None so that the models that are only
-                # correct for double precision are not tested using
-                # single precision.  The choice is determined by the
-                # presence of *single=False* in the model file.
-                test = ModelTestCase(test_name, model_info,
-                                     test_method_name,
-                                     platform="ocl", dtype=None,
-                                     stash=stash)
-                #print("defining", test_name)
-                suite.addTest(test)
-
-    return suite
+
+        # test using opencl if desired and available
+        if 'opencl' in loaders and use_opencl():
+            test_name = "%s-opencl"%model_info.name
+            test_method_name = "test_%s_opencl" % model_info.id
+            # Using dtype=None so that the models that are only
+            # correct for double precision are not tested using
+            # single precision.  The choice is determined by the
+            # presence of *single=False* in the model file.
+            test = ModelTestCase(test_name, model_info,
+                                    test_method_name,
+                                    platform="ocl", dtype=None,
+                                    stash=stash)
+            #print("defining", test_name)
+            suite.addTest(test)
+
 
 def _hide_model_case_from_nose():
@@ -348,11 +351,23 @@
     return abs(target-actual)/shift < 1.5*10**-digits
 
-def run_one(model):
-    # type: (str) -> str
-    """
-    Run the tests for a single model, printing the results to stdout.
-
-    *model* can by a python file, which is handy for checking user defined
-    plugin models.
+# CRUFT: old interface; should be deprecated and removed
+def run_one(model_name):
+    # msg = "use check_model(model_info) rather than run_one(model_name)"
+    # warnings.warn(msg, category=DeprecationWarning, stacklevel=2)
+    try:
+        model_info = load_model_info(model_name)
+    except Exception:
+        output = traceback.format_exc()
+        return output
+
+    success, output = check_model(model_info)
+    return output
+
+def check_model(model_info):
+    # type: (ModelInfo) -> str
+    """
+    Run the tests for a single model, capturing the output.
+
+    Returns success status and the output string.
     """
     # Note that running main() directly did not work from within the
@@ -369,11 +384,6 @@
     # Build a test suite containing just the model
     loaders = ['opencl'] if use_opencl() else ['dll']
-    models = [model]
-    try:
-        suite = make_suite(loaders, models)
-    except Exception:
-        import traceback
-        stream.writeln(traceback.format_exc())
-        return
+    suite = unittest.TestSuite()
+    _add_model_to_suite(loaders, suite, model_info)
 
     # Warn if there are no user defined tests.
@@ -390,5 +400,5 @@
     for test in suite:
         if not test.info.tests:
-            stream.writeln("Note: %s has no user defined tests."%model)
+            stream.writeln("Note: %s has no user defined tests."%model_info.name)
         break
     else:
@@ -406,5 +416,5 @@
     output = stream.getvalue()
     stream.close()
-    return output
+    return result.wasSuccessful(), output
 
 
Index: sasmodels/modelinfo.py
===================================================================
--- sasmodels/modelinfo.py	(revision 7b9e4ddca7ad3551d8fe3696d2afe75b0018599e)
+++ sasmodels/modelinfo.py	(revision bd547d022d7afce0c8a39ebfc0c2551bc9339186)
@@ -466,5 +466,5 @@
         self.is_asymmetric = any(p.name == 'psi' for p in self.kernel_parameters)
         self.magnetism_index = [k for k, p in enumerate(self.call_parameters)
-                                if p.id.startswith('M0:')]
+                                if p.id.endswith('_M0')]
 
         self.pd_1d = set(p.name for p in self.call_parameters
@@ -586,9 +586,9 @@
         if self.nmagnetic > 0:
             full_list.extend([
-                Parameter('up:frac_i', '', 0., [0., 1.],
+                Parameter('up_frac_i', '', 0., [0., 1.],
                           'magnetic', 'fraction of spin up incident'),
-                Parameter('up:frac_f', '', 0., [0., 1.],
+                Parameter('up_frac_f', '', 0., [0., 1.],
                           'magnetic', 'fraction of spin up final'),
-                Parameter('up:angle', 'degrees', 0., [0., 360.],
+                Parameter('up_angle', 'degrees', 0., [0., 360.],
                           'magnetic', 'spin up angle'),
             ])
@@ -596,9 +596,9 @@
             for p in slds:
                 full_list.extend([
-                    Parameter('M0:'+p.id, '1e-6/Ang^2', 0., [-np.inf, np.inf],
+                    Parameter(p.id+'_M0', '1e-6/Ang^2', 0., [-np.inf, np.inf],
                               'magnetic', 'magnetic amplitude for '+p.description),
-                    Parameter('mtheta:'+p.id, 'degrees', 0., [-90., 90.],
+                    Parameter(p.id+'_mtheta', 'degrees', 0., [-90., 90.],
                               'magnetic', 'magnetic latitude for '+p.description),
-                    Parameter('mphi:'+p.id, 'degrees', 0., [-180., 180.],
+                    Parameter(p.id+'_mphi', 'degrees', 0., [-180., 180.],
                               'magnetic', 'magnetic longitude for '+p.description),
                 ])
@@ -640,6 +640,6 @@
 
         Parameters marked as sld will automatically have a set of associated
-        magnetic parameters (m0:p, mtheta:p, mphi:p), as well as polarization
-        information (up:theta, up:frac_i, up:frac_f).
+        magnetic parameters (p_M0, p_mtheta, p_mphi), as well as polarization
+        information (up_theta, up_frac_i, up_frac_f).
         """
         # control parameters go first
@@ -668,7 +668,7 @@
             result.append(expanded_pars[name])
             if is2d:
-                for tag in 'M0:', 'mtheta:', 'mphi:':
-                    if tag+name in expanded_pars:
-                        result.append(expanded_pars[tag+name])
+                for tag in '_M0', '_mtheta', '_mphi':
+                    if name+tag in expanded_pars:
+                        result.append(expanded_pars[name+tag])
 
         # Gather the user parameters in order
@@ -703,9 +703,9 @@
                 append_group(p.id)
 
-        if is2d and 'up:angle' in expanded_pars:
+        if is2d and 'up_angle' in expanded_pars:
             result.extend([
-                expanded_pars['up:frac_i'],
-                expanded_pars['up:frac_f'],
-                expanded_pars['up:angle'],
+                expanded_pars['up_frac_i'],
+                expanded_pars['up_frac_f'],
+                expanded_pars['up_angle'],
             ])
 
@@ -793,4 +793,6 @@
     info.structure_factor = getattr(kernel_module, 'structure_factor', False)
     info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y'])
+    # Note: custom.load_custom_kernel_module assumes the C sources are defined
+    # by this attribute.
     info.source = getattr(kernel_module, 'source', [])
     info.c_code = getattr(kernel_module, 'c_code', None)
@@ -1014,3 +1016,8 @@
                          for k in range(control+1, p.length+1)
                          if p.length > 1)
+            for p in self.parameters.kernel_parameters:
+                if p.length > 1 and p.type == "sld":
+                    for k in range(control+1, p.length+1):
+                        base = p.id+str(k)
+                        hidden.update((base+"_M0", base+"_mtheta", base+"_mphi"))
         return hidden
Index: sasmodels/models/bcc_paracrystal.py
===================================================================
--- sasmodels/models/bcc_paracrystal.py	(revision e7e9231f417c4c697a93a5690bbe7f81fad8f175)
+++ sasmodels/models/bcc_paracrystal.py	(revision da7b26b2f3ecdd4d6eea063846bb148e3def2477)
@@ -1,3 +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 
+             SasView Developers. September 2018.*
+
 Definition
 ----------
@@ -13,5 +18,4 @@
 
     I(q) = \frac{\text{scale}}{V_p} V_\text{lattice} P(q) Z(q)
-
 
 where *scale* is the volume fraction of spheres, $V_p$ is the volume of the
@@ -97,5 +101,5 @@
 
 Authorship and Verification
-----------------------------
+---------------------------
 
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
Index: sasmodels/models/be_polyelectrolyte.py
===================================================================
--- sasmodels/models/be_polyelectrolyte.py	(revision ef07e950965796e2a3a8156fcade98ce5e403001)
+++ sasmodels/models/be_polyelectrolyte.py	(revision ca77fc11c7ce03a565ec1aaf8755f560748546d0)
@@ -1,3 +1,5 @@
 r"""
+.. note:: Please read the Validation section below.
+
 Definition
 ----------
@@ -11,14 +13,14 @@
 
     I(q) = K\frac{q^2+k^2}{4\pi L_b\alpha ^2}
-    \frac{1}{1+r_{0}^2(q^2+k^2)(q^2-12hC_a/b^2)} + background
+    \frac{1}{1+r_{0}^4(q^2+k^2)(q^2-12hC_a/b^2)} + background
 
     k^2 = 4\pi L_b(2C_s + \alpha C_a)
 
-    r_{0}^2 = \frac{1}{\alpha \sqrt{C_a} \left( b/\sqrt{48\pi L_b}\right)}
+    r_{0}^2 = \frac{b}{\alpha \sqrt{C_a 48\pi L_b}}
 
 where
 
 $K$ is the contrast factor for the polymer which is defined differently than in
-other models and is given in barns where $1 barn = 10^{-24} cm^2$.  $K$ is
+other models and is given in barns where 1 $barn = 10^{-24}$ $cm^2$.  $K$ is
 defined as:
 
@@ -29,20 +31,29 @@
     a = b_p - (v_p/v_s) b_s
 
-where $b_p$ and $b_s$ are sum of the scattering lengths of the atoms
-constituting the monomer of the polymer and the sum of the scattering lengths
-of the atoms constituting the solvent molecules respectively, and $v_p$ and
-$v_s$ are the partial molar volume of the polymer and the solvent respectively
-
-$L_b$ is the Bjerrum length(|Ang|) - **Note:** This parameter needs to be
-kept constant for a given solvent and temperature!
-
-$h$ is the virial parameter (|Ang^3|/mol) - **Note:** See [#Borue]_ for the
-correct interpretation of this parameter.  It incorporates second and third
-virial coefficients and can be Negative.
-
-$b$ is the monomer length(|Ang|), $C_s$ is the concentration of monovalent
-salt(mol/L), $\alpha$ is the ionization degree (ionization degree : ratio of
-charged monomers  to total number of monomers), $C_a$ is the polymer molar
-concentration(mol/L), and $background$ is the incoherent background.
+where:
+
+- $b_p$ and $b_s$ are **sum of the scattering lengths of the atoms**
+  constituting the polymer monomer and the solvent molecules, respectively.
+
+- $v_p$ and $v_s$ are the partial molar volume of the polymer and the 
+  solvent, respectively.
+
+- $L_b$ is the Bjerrum length (|Ang|) - **Note:** This parameter needs to be
+  kept constant for a given solvent and temperature!
+
+- $h$ is the virial parameter (|Ang^3|) - **Note:** See [#Borue]_ for the
+  correct interpretation of this parameter.  It incorporates second and third
+  virial coefficients and can be *negative*.
+
+- $b$ is the monomer length (|Ang|).
+
+- $C_s$ is the concentration of monovalent salt(1/|Ang^3| - internally converted from mol/L).
+
+- $\alpha$ is the degree of ionization (the ratio of charged monomers to the total 
+  number of monomers)
+
+- $C_a$ is the polymer molar concentration (1/|Ang^3| - internally converted from mol/L)
+
+- $background$ is the incoherent background.
 
 For 2D data the scattering intensity is calculated in the same way as 1D,
@@ -52,4 +63,22 @@
 
     q = \sqrt{q_x^2 + q_y^2}
+
+Validation
+----------
+
+As of the last revision, this code is believed to be correct.  However it
+needs further validation and should be used with caution at this time.  The
+history of this code goes back to a 1998 implementation. It was recently noted
+that in that implementation, while both the polymer concentration and salt
+concentration were converted from experimental units of mol/L to more
+dimensionally useful units of 1/|Ang^3|, only the converted version of the
+polymer concentration was actually being used in the calculation while the
+unconverted salt concentration (still in apparent units of mol/L) was being 
+used.  This was carried through to Sasmodels as used for SasView 4.1 (though 
+the line of code converting the salt concentration to the new units was removed 
+somewhere along the line). Simple dimensional analysis of the calculation shows 
+that the converted salt concentration should be used, which the original code 
+suggests was the intention, so this has now been corrected (for SasView 4.2). 
+Once better validation has been performed this note will be removed.
 
 References
@@ -66,6 +95,6 @@
 
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
-* **Last Modified by:** Paul Kienzle **Date:** July 24, 2016
-* **Last Reviewed by:** Paul Butler and Richard Heenan **Date:** October 07, 2016
+* **Last Modified by:** Paul Butler **Date:** September 25, 2018
+* **Last Reviewed by:** Paul Butler **Date:** September 25, 2018
 """
 
@@ -92,5 +121,5 @@
     ["contrast_factor",       "barns",   10.0,  [-inf, inf], "", "Contrast factor of the polymer"],
     ["bjerrum_length",        "Ang",      7.1,  [0, inf],    "", "Bjerrum length"],
-    ["virial_param",          "Ang^3/mol", 12.0,  [-inf, inf], "", "Virial parameter"],
+    ["virial_param",          "Ang^3", 12.0,  [-inf, inf], "", "Virial parameter"],
     ["monomer_length",        "Ang",     10.0,  [0, inf],    "", "Monomer length"],
     ["salt_concentration",    "mol/L",    0.0,  [-inf, inf], "", "Concentration of monovalent salt"],
@@ -102,29 +131,27 @@
 
 def Iq(q,
-       contrast_factor=10.0,
-       bjerrum_length=7.1,
-       virial_param=12.0,
-       monomer_length=10.0,
-       salt_concentration=0.0,
-       ionization_degree=0.05,
-       polymer_concentration=0.7):
+       contrast_factor,
+       bjerrum_length,
+       virial_param,
+       monomer_length,
+       salt_concentration,
+       ionization_degree,
+       polymer_concentration):
     """
-    :param q:                     Input q-value
-    :param contrast_factor:       Contrast factor of the polymer
-    :param bjerrum_length:        Bjerrum length
-    :param virial_param:          Virial parameter
-    :param monomer_length:        Monomer length
-    :param salt_concentration:    Concentration of monovalent salt
-    :param ionization_degree:     Degree of ionization
-    :param polymer_concentration: Polymer molar concentration
-    :return:                      1-D intensity
+    :params: see parameter table
+    :return: 1-D form factor for polyelectrolytes in low salt
+    
+    parameter names, units, default values, and behavior (volume, sld etc) are
+    defined in the parameter table.  The concentrations are converted from
+    experimental mol/L to dimensionaly useful 1/A3 in first two lines
     """
 
-    concentration = polymer_concentration * 6.022136e-4
-
-    k_square = 4.0 * pi * bjerrum_length * (2*salt_concentration +
-                                            ionization_degree * concentration)
-
-    r0_square = 1.0/ionization_degree/sqrt(concentration) * \
+    concentration_pol = polymer_concentration * 6.022136e-4
+    concentration_salt = salt_concentration * 6.022136e-4
+
+    k_square = 4.0 * pi * bjerrum_length * (2*concentration_salt +
+                                            ionization_degree * concentration_pol)
+
+    r0_square = 1.0/ionization_degree/sqrt(concentration_pol) * \
                 (monomer_length/sqrt((48.0*pi*bjerrum_length)))
 
@@ -133,5 +160,5 @@
 
     term2 = 1.0 + r0_square**2 * (q**2 + k_square) * \
-        (q**2 - (12.0 * virial_param * concentration/(monomer_length**2)))
+        (q**2 - (12.0 * virial_param * concentration_pol/(monomer_length**2)))
 
     return term1/term2
@@ -174,4 +201,9 @@
 
     # Accuracy tests based on content in test/utest_other_models.py
+    # Note that these should some day be validated beyond this self validation
+    # (circular reasoning). -- i.e. the "good value," at least for those with
+    # non zero salt concentrations, were obtained by running the current
+    # model in SasView and copying the appropriate result here.
+    #    PDB -- Sep 26, 2018
     [{'contrast_factor':       10.0,
       'bjerrum_length':         7.1,
@@ -184,33 +216,32 @@
      }, 0.001, 0.0948379],
 
-    # Additional tests with larger range of parameters
     [{'contrast_factor':       10.0,
       'bjerrum_length':       100.0,
       'virial_param':           3.0,
-      'monomer_length':         1.0,
-      'salt_concentration':    10.0,
-      'ionization_degree':      2.0,
-      'polymer_concentration': 10.0,
+      'monomer_length':         5.0,
+      'salt_concentration':     1.0,
+      'ionization_degree':      0.1,
+      'polymer_concentration':  1.0,
       'background':             0.0,
-     }, 0.1, -3.75693800588],
+     }, 0.1, 0.253469484],
 
     [{'contrast_factor':       10.0,
       'bjerrum_length':       100.0,
       'virial_param':           3.0,
-      'monomer_length':         1.0,
-      'salt_concentration':    10.0,
-      'ionization_degree':      2.0,
-      'polymer_concentration': 10.0,
-      'background':           100.0
-     }, 5.0, 100.029142149],
+      'monomer_length':         5.0,
+      'salt_concentration':     1.0,
+      'ionization_degree':      0.1,
+      'polymer_concentration':  1.0,
+      'background':             1.0,
+     }, 0.05, 1.738358122],
 
     [{'contrast_factor':     100.0,
       'bjerrum_length':       10.0,
-      'virial_param':        180.0,
-      'monomer_length':        1.0,
+      'virial_param':         12.0,
+      'monomer_length':       10.0,
       'salt_concentration':    0.1,
       'ionization_degree':     0.5,
       'polymer_concentration': 0.1,
-      'background':             0.0,
-     }, 200., 1.80664667511e-06],
+      'background':           0.01,
+     }, 0.5, 0.012881893],
     ]
Index: sasmodels/models/fcc_paracrystal.py
===================================================================
--- sasmodels/models/fcc_paracrystal.py	(revision e7e9231f417c4c697a93a5690bbe7f81fad8f175)
+++ sasmodels/models/fcc_paracrystal.py	(revision da7b26b2f3ecdd4d6eea063846bb148e3def2477)
@@ -3,4 +3,9 @@
 #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
+             SasView Developers. September 2018.*
+
 Definition
 ----------
@@ -10,5 +15,5 @@
 negligible, and the size of the paracrystal is infinitely large.
 Paracrystalline distortion is assumed to be isotropic and characterized by
-a Gaussian distribution. 
+a Gaussian distribution.
 
 The scattering intensity $I(q)$ is calculated as
Index: sasmodels/models/lib/sas_gammainc.c
===================================================================
--- sasmodels/models/lib/sas_gammainc.c	(revision fba9ca0239129ed28284a0b435ba2a9524913f36)
+++ sasmodels/models/lib/sas_gammainc.c	(revision fba9ca0239129ed28284a0b435ba2a9524913f36)
@@ -0,0 +1,458 @@
+#if FLOAT_SIZE>4  // double precision
+// based on cephes/double/igam.c
+double gammaln(double x);
+double gammainc(double a, double x);
+double gammaincc(double a, double x);
+
+double cephes_igamc(double a, double x);
+double cephes_igam(double a, double x);
+double cephes_lgam(double x);
+double cephes_lgam2(double x);
+
+double gammainc(double a, double x)
+{
+    if ((x <= 0) || ( a <= 0)) return 0.0;
+    if ((x > 1.0) && (x > a)) return 1.0 - cephes_igamc(a, x);
+    return cephes_igam(a, x);
+}
+
+double gammaincc(double a, double x)
+{
+    if ((x <= 0) || (a <= 0)) return 1.0;
+    if ((x < 1.0) || (x < a)) return 1.0 - cephes_igam(a, x);
+    return cephes_igamc(a, x);
+}
+
+double gammaln(double x)
+{
+    if (isnan(x)) return(x);
+    if (!isfinite(x)) return(INFINITY);
+    return cephes_lgam(x);
+}
+
+double cephes_igamc(double a, double x)
+{
+    const double MACHEP = 1.11022302462515654042E-16; // IEEE 2**-53
+    const double MAXLOG = 7.09782712893383996843E2; // IEEE log(2**1024) denormalized
+    const double BIG = 4.503599627370496e15;
+    const double BIGINV = 2.22044604925031308085e-16;
+    double ans, ax, c, yc, r, t, y, z;
+    double pk, pkm1, pkm2, qk, qkm1, qkm2;
+
+    /* Compute  x**a * exp(-x) / gamma(a)  */
+    ax = a * log(x) - x - cephes_lgam(a);
+    if (ax < -MAXLOG) return 0.0;  // underflow
+    ax = exp(ax);
+
+    /* continued fraction */
+    y = 1.0 - a;
+    z = x + y + 1.0;
+    c = 0.0;
+    pkm2 = 1.0;
+    qkm2 = x;
+    pkm1 = x + 1.0;
+    qkm1 = z * x;
+    ans = pkm1/qkm1;
+
+    do {
+        c += 1.0;
+        y += 1.0;
+        z += 2.0;
+        yc = y * c;
+        pk = pkm1 * z  -  pkm2 * yc;
+        qk = qkm1 * z  -  qkm2 * yc;
+        if (qk != 0) {
+            r = pk/qk;
+            t = fabs( (ans - r)/r );
+            ans = r;
+        } else {
+            t = 1.0;
+        }
+        pkm2 = pkm1;
+        pkm1 = pk;
+        qkm2 = qkm1;
+        qkm1 = qk;
+        if (fabs(pk) > BIG) {
+            pkm2 *= BIGINV;
+            pkm1 *= BIGINV;
+            qkm2 *= BIGINV;
+            qkm1 *= BIGINV;
+        }
+    } while( t > MACHEP );
+
+    return( ans * ax );
+}
+
+
+double cephes_igam(double a, double x)
+{
+    const double MACHEP = 1.11022302462515654042E-16; // IEEE 2**-53
+    const double MAXLOG = 7.09782712893383996843E2; // IEEE log(2**1024) denormalized
+    double ans, ax, c, r;
+
+    /* Compute  x**a * exp(-x) / gamma(a)  */
+    ax = a * log(x) - x - cephes_lgam(a);
+    if (ax < -MAXLOG) return 0.0; // underflow
+    ax = exp(ax);
+
+    /* power series */
+    r = a;
+    c = 1.0;
+    ans = 1.0;
+
+    do {
+        r += 1.0;
+        c *= x/r;
+        ans += c;
+    } while (c/ans > MACHEP);
+
+    return ans * ax/a;
+}
+
+/* Logarithm of gamma function */
+
+
+double cephes_lgam(double x)
+{
+    const double LOGPI = 1.14472988584940017414; // log(pi)
+    int sgngam = 1;
+    double p, q, w, z;
+    int i;
+
+    if (x < -34.0) {
+        q = -x;
+        w = cephes_lgam2(q);
+        p = floor(q);
+        if (p == q) return INFINITY;
+        i = p;
+        sgngam = ((i&1) == 0) ? -1 : 1;
+        z = q - p;
+        if (z > 0.5) {
+            p += 1.0;
+            z = p - q;
+        }
+        z = q * sin(M_PI * z);
+        if (z == 0.0) return INFINITY;
+        z = LOGPI - log(z) - w;
+        return z;
+    } else {
+        return cephes_lgam2(x);
+    }
+}
+
+double cephes_lgam2(double x)
+{
+    const double LS2PI = 0.91893853320467274178; // log(sqrt(2*pi))
+    const double MAXLGM = 2.556348e305;
+    int sgngam = 1;
+    double p, q, u, z;
+    double A, B, C;
+
+    if (x < 13.0) {
+        z = 1.0;
+        p = 0.0;
+        u = x;
+        while (u >= 3.0) {
+            p -= 1.0;
+            u = x + p;
+            z *= u;
+        }
+        while (u < 2.0) {
+            if (u == 0.0) return INFINITY;
+            z /= u;
+            p += 1.0;
+            u = x + p;
+        }
+        if (z < 0.0) {
+            sgngam = -1;
+            z = -z;
+        } else {
+            sgngam = 1;
+        }
+        if (u == 2.0) return log(z);
+        p -= 2.0;
+        x = x + p;
+        B = ((((((
+            -1.37825152569120859100E3)*x
+            -3.88016315134637840924E4)*x
+            -3.31612992738871184744E5)*x
+            -1.16237097492762307383E6)*x
+            -1.72173700820839662146E6)*x
+            -8.53555664245765465627E5);
+        C = ((((((
+            /* 1.00000000000000000000E0)* */x
+            -3.51815701436523470549E2)*x
+            -1.70642106651881159223E4)*x
+            -2.20528590553854454839E5)*x
+            -1.13933444367982507207E6)*x
+            -2.53252307177582951285E6)*x
+            -2.01889141433532773231E6);
+        p = x * B / C;
+        return log(z) + p;
+    }
+
+    if (x > MAXLGM) return sgngam * INFINITY;
+
+    q = (x - 0.5) * log(x) - x + LS2PI;
+    if (x > 1.0e8) return q;
+
+    p = 1.0/(x*x);
+    if (x >= 1000.0) {
+        q += ((7.9365079365079365079365e-4 * p
+            - 2.7777777777777777777778e-3) * p
+            + 0.0833333333333333333333) / x;
+    } else {
+        A = (((((
+            +8.11614167470508450300E-4)*p
+            -5.95061904284301438324E-4)*p
+            +7.93650340457716943945E-4)*p
+            -2.77777777730099687205E-3)*p
+            +8.33333333333331927722E-2);
+        q += A / x;
+    }
+    return q;
+}
+
+#else // single precision
+
+// based on cephes/float/igam.c
+float gammalnf(float x);
+float gammaincf(float a, float x);
+float gammainccf(float a, float x);
+
+float cephes_igamcf(float a, float x);
+float cephes_igamf(float a, float x);
+float cephes_lgamf(float x);
+float cephes_lgam2f(float x);
+
+// Note: original uses logf, fabsf, floorf and sinf
+
+float gammaincf(float a, float x)
+{
+    if ((x <= 0) || ( a <= 0)) return 0.0;
+    if ((x > 1.0) && (x > a)) return 1.0 - cephes_igamcf(a, x);
+    return cephes_igamf(a, x);
+}
+
+float gammainccf(float a, float x)
+{
+    if ((x <= 0) || (a <= 0)) return 1.0;
+    if ((x < 1.0) || (x < a)) return 1.0 - cephes_igamf(a, x);
+    return cephes_igamcf(a, x);
+}
+
+float gammalnf(float x)
+{
+    if (isnan(x)) return(x);
+    if (!isfinite(x)) return(INFINITY);
+    return cephes_lgamf(x);
+}
+
+float cephes_igamcf(float a, float x)
+{
+    const float MAXLOGF = 88.72283905206835;
+    const float MACHEPF = 5.9604644775390625E-8;
+    const float BIG = 16777216.;
+    float ans, c, yc, ax, y, z;
+    float pk, pkm1, pkm2, qk, qkm1, qkm2;
+    float r, t;
+
+    ax = a * log(x) - x - cephes_lgamf(a);
+    if (ax < -MAXLOGF) return 0.0; // underflow
+    ax = expf(ax);
+
+    /* continued fraction */
+    y = 1.0 - a;
+    z = x + y + 1.0;
+    c = 0.0;
+    pkm2 = 1.0;
+    qkm2 = x;
+    pkm1 = x + 1.0;
+    qkm1 = z * x;
+    ans = pkm1/qkm1;
+
+    do {
+        c += 1.0;
+        y += 1.0;
+        z += 2.0;
+        yc = y * c;
+        pk = pkm1 * z  -  pkm2 * yc;
+        qk = qkm1 * z  -  qkm2 * yc;
+        if (qk != 0) {
+            r = pk/qk;
+            t = fabs((ans - r)/r);
+            ans = r;
+        } else {
+            t = 1.0;
+        }
+        pkm2 = pkm1;
+        pkm1 = pk;
+        qkm2 = qkm1;
+        qkm1 = qk;
+        if (fabs(pk) > BIG) {
+            pkm2 *= MACHEPF;
+            pkm1 *= MACHEPF;
+            qkm2 *= MACHEPF;
+            qkm1 *= MACHEPF;
+        }
+    } while (t > MACHEPF);
+
+    return ans * ax;
+}
+
+float cephes_igamf(float a, float x)
+{
+    const float MAXLOGF = 88.72283905206835;
+    const float MACHEPF = 5.9604644775390625E-8;
+    float ans, ax, c, r;
+
+    /* Compute  x**a * exp(-x) / gamma(a)  */
+    ax = a * log(x) - x - cephes_lgamf(a);
+    if (ax < -MAXLOGF) return 0.0; // underflow
+    ax = expf(ax);
+
+    /* power series */
+    r = a;
+    c = 1.0;
+    ans = 1.0;
+
+    do {
+        r += 1.0;
+        c *= x/r;
+        ans += c;
+    } while (c/ans > MACHEPF);
+
+    return ans * ax/a;
+}
+
+float cephes_lgamf(float x)
+{
+    const float PIINV =  0.318309886183790671538;
+    float p, q, w, z;
+    int i;
+    int sgngamf;
+
+    if (x < 0.0) {
+        q = -x;
+        w = cephes_lgam2f(q);
+        p = floor(q);
+        if (p == q) return INFINITY;
+        i = p;
+        sgngamf = ((i&1) == 0) ? -1 : 1;
+        z = q - p;
+        if (z > 0.5) {
+            p += 1.0;
+            z = p - q;
+        }
+        z = q * sin(M_PI * z);
+        if (z == 0.0) return sgngamf * INFINITY;
+        z = -log(PIINV*z) - w;
+        return z;
+    } else {
+        return cephes_lgam2f(x);
+    }
+}
+
+float cephes_lgam2f(float x)
+{
+    const float LS2PI = 0.91893853320467274178; // log(sqrt(2*pi))
+    const float MAXLGM = 2.035093e36;
+    float p, q, z;
+    float nx, tx;
+    int direction;
+    int sgngamf = 1;
+
+    if (x < 6.5) {
+        direction = 0;
+        z = 1.0;
+        tx = x;
+        nx = 0.0;
+        if (x >= 1.5) {
+            while (tx > 2.5) {
+                nx -= 1.0;
+                tx = x + nx;
+                z *=tx;
+            }
+            x += nx - 2.0;
+        } else if (x >= 1.25) {
+            z *= x;
+            x -= 1.0; /* x + 1 - 2 */
+            direction = 1;
+        } else if (x >= 0.75) {
+            x -= 1.0;
+            /* log gamma(x+1), -.25 < x < .25 */
+            // p = x * polevlf( x, C, 7 );
+            p = ((((((((
+                +1.369488127325832E-001)*x
+                -1.590086327657347E-001)*x
+                +1.692415923504637E-001)*x
+                -2.067882815621965E-001)*x
+                +2.705806208275915E-001)*x
+                -4.006931650563372E-001)*x
+                +8.224670749082976E-001)*x
+                -5.772156501719101E-001)*x;
+            q = 0.0;
+            return p + q;
+        } else {
+            while (tx < 1.5) {
+                if (tx == 0.0) return sgngamf * INFINITY;
+                z *=tx;
+                nx += 1.0;
+                tx = x + nx;
+            }
+            direction = 1;
+            x += nx - 2.0;
+        }
+
+        /* log gamma(x+2), -.5 < x < .5 */
+        // p = x * polevlf( x, B, 7 );
+        p = ((((((((
+            +6.055172732649237E-004)*x
+            -1.311620815545743E-003)*x
+            +2.863437556468661E-003)*x
+            -7.366775108654962E-003)*x
+            +2.058355474821512E-002)*x
+            -6.735323259371034E-002)*x
+            +3.224669577325661E-001)*x
+            +4.227843421859038E-001)*x;
+
+        if (z < 0.0) {
+            sgngamf = -1;
+            z = -z;
+        } else {
+            sgngamf = 1;
+        }
+        q = log(z);
+        if (direction) q = -q;
+        return p + q;
+    } else if (x > MAXLGM) {
+        return sgngamf * INFINITY;
+    } else {
+        /* Note, though an asymptotic formula could be used for x >= 3,
+        * there is cancellation error in the following if x < 6.5.
+        */
+        q = LS2PI - x;
+        q += (x - 0.5) * log(x);
+
+        if (x <= 1.0e4) {
+            z = 1.0/x;
+            p = z * z;
+            q += ((6.789774945028216E-004 * p
+                - 2.769887652139868E-003 ) * p
+                +  8.333316229807355E-002 ) * z;
+        }
+        return q;
+    }
+}
+
+#endif // !single precision
+
+#if FLOAT_SIZE>4
+#define sas_gammaln gammaln
+#define sas_gammainc gammainc
+#define sas_gammaincc gammaincc
+#else
+#define sas_gammaln gammalnf
+#define sas_gammainc gammaincf
+#define sas_gammaincc gammainccf
+#endif
Index: sasmodels/models/sc_paracrystal.py
===================================================================
--- sasmodels/models/sc_paracrystal.py	(revision 0b906ea4a50ccd003f35113867d84a93d788de83)
+++ sasmodels/models/sc_paracrystal.py	(revision da7b26b2f3ecdd4d6eea063846bb148e3def2477)
@@ -1,3 +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
+             SasView Developers. September 2018.*
+
 Definition
 ----------
@@ -9,5 +14,5 @@
 by a Gaussian distribution.
 
-he scattering intensity $I(q)$ is calculated as
+The scattering intensity $I(q)$ is calculated as
 
 .. math::
@@ -21,5 +26,6 @@
 
 Equation (16) of the 1987 reference\ [#CIT1987]_ is used to calculate $Z(q)$,
-using equations (13)-(15) from the 1987 paper\ [#CIT1987]_ for Z1, Z2, and Z3.
+using equations (13)-(15) from the 1987 paper\ [#CIT1987]_ for $Z1$, $Z2$, and
+$Z3$.
 
 The lattice correction (the occupied volume of the lattice) for a simple cubic
Index: sasmodels/models/spinodal.py
===================================================================
--- sasmodels/models/spinodal.py	(revision 475ff585df3808f7c560003331e4cf2e1ea87fc3)
+++ sasmodels/models/spinodal.py	(revision 93fe8a1613531fbb5a9270c8cee68cad9b4642cb)
@@ -12,5 +12,18 @@
 where $x=q/q_0$, $q_0$ is the peak position, $I_{max}$ is the intensity 
 at $q_0$ (parameterised as the $scale$ parameter), and $B$ is a flat 
-background. The spinodal wavelength is given by $2\pi/q_0$. 
+background. The spinodal wavelength, $\Lambda$, is given by $2\pi/q_0$. 
+
+The definition of $I_{max}$ in the literature varies. Hashimoto *et al* (1991) 
+define it as 
+
+.. math::
+    I_{max} = \Lambda^3\Delta\rho^2
+    
+whereas Meier & Strobl (1987) give 
+
+.. math::
+    I_{max} = V_z\Delta\rho^2
+    
+where $V_z$ is the volume per monomer unit.
 
 The exponent $\gamma$ is equal to $d+1$ for off-critical concentration 
@@ -28,6 +41,14 @@
 
 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).
+Growth rates of droplets and scaling properties of autocorrelation functions. 
+Physica A 123, 497 (1984).
+
+H. Meier & G. Strobl. Small-Angle X-ray Scattering Study of Spinodal 
+Decomposition in Polystyrene/Poly(styrene-co-bromostyrene) Blends. 
+Macromolecules 20, 649-654 (1987).
+
+T. Hashimoto, M. Takenaka & H. Jinnai. Scattering Studies of Self-Assembling 
+Processes of Polymer Blends in Spinodal Decomposition. 
+J. Appl. Cryst. 24, 457-466 (1991).
 
 Revision History
@@ -35,5 +56,5 @@
 
 * **Author:**  Dirk Honecker **Date:** Oct 7, 2016
-* **Revised:** Steve King    **Date:** Sep 7, 2018
+* **Revised:** Steve King    **Date:** Oct 25, 2018
 """
 
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision d533590d9857d7aef962616cfdaa6ecfc2e6751d)
+++ sasmodels/sasview_model.py	(revision bd547d022d7afce0c8a39ebfc0c2551bc9339186)
@@ -62,6 +62,10 @@
 #: set of defined models (standard and custom)
 MODELS = {}  # type: Dict[str, SasviewModelType]
+# TODO: remove unused MODEL_BY_PATH cache once sasview no longer references it
 #: custom model {path: model} mapping so we can check timestamps
 MODEL_BY_PATH = {}  # type: Dict[str, SasviewModelType]
+#: Track modules that we have loaded so we can determine whether the model
+#: has changed since we last reloaded.
+_CACHED_MODULE = {}  # type: Dict[str, "module"]
 
 def find_model(modelname):
@@ -106,13 +110,19 @@
     Load a custom model given the model path.
     """
-    model = MODEL_BY_PATH.get(path, None)
-    if model is not None and model.timestamp == getmtime(path):
-        #logger.info("Model already loaded %s", path)
-        return model
-
     #logger.info("Loading model %s", path)
+
+    # Load the kernel module.  This may already be cached by the loader, so
+    # only requires checking the timestamps of the dependents.
     kernel_module = custom.load_custom_kernel_module(path)
-    if hasattr(kernel_module, 'Model'):
-        model = kernel_module.Model
+
+    # Check if the module has changed since we last looked.
+    reloaded = kernel_module != _CACHED_MODULE.get(path, None)
+    _CACHED_MODULE[path] = kernel_module
+
+    # Turn the module into a model.  We need to do this in even if the
+    # model has already been loaded so that we can determine the model
+    # name and retrieve it from the MODELS cache.
+    model = getattr(kernel_module, 'Model', None)
+    if model is not None:
         # Old style models do not set the name in the class attributes, so
         # set it here; this name will be overridden when the object is created
@@ -127,5 +137,4 @@
         model_info = modelinfo.make_model_info(kernel_module)
         model = make_model_from_info(model_info)
-    model.timestamp = getmtime(path)
 
     # If a model name already exists and we are loading a different model,
@@ -143,7 +152,9 @@
                     _previous_name, model.name, model.filename)
 
-    MODELS[model.name] = model
-    MODEL_BY_PATH[path] = model
-    return model
+    # Only update the model if the module has changed
+    if reloaded or model.name not in MODELS:
+        MODELS[model.name] = model
+
+    return MODELS[model.name]
 
 
@@ -372,4 +383,10 @@
             hidden.add('background')
             self._model_info.parameters.defaults['background'] = 0.
+
+        # Update the parameter lists to exclude any hidden parameters
+        self.magnetic_params = tuple(pname for pname in self.magnetic_params
+                                     if pname not in hidden)
+        self.orientation_params = tuple(pname for pname in self.orientation_params
+                                        if pname not in hidden)
 
         self._persistency_dict = {}
@@ -786,4 +803,14 @@
             return value, [value], [1.0]
 
+    @classmethod
+    def runTests(cls):
+        """
+        Run any tests built into the model and captures the test output.
+
+        Returns success flag and output
+        """
+        from .model_test import check_model
+        return check_model(cls._model_info)
+
 def test_cylinder():
     # type: () -> float
@@ -873,5 +900,5 @@
     Model = _make_standard_model('sphere')
     model = Model()
-    model.setParam('M0:sld', 8)
+    model.setParam('sld_M0', 8)
     q = np.linspace(-0.35, 0.35, 500)
     qx, qy = np.meshgrid(q, q)
Index: sasmodels/special.py
===================================================================
--- sasmodels/special.py	(revision df69efa2cdb726d21ce081492d47b44552fc0d19)
+++ sasmodels/special.py	(revision fba9ca0239129ed28284a0b435ba2a9524913f36)
@@ -113,4 +113,16 @@
         The standard math function, tgamma(x) is unstable for $x < 1$
         on some platforms.
+
+    sas_gammaln(x):
+        log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$.
+
+        The standard math function, lgamma(x), is incorrect for single
+        precision on some platforms.
+
+    sas_gammainc(a, x), sas_gammaincc(a, x):
+        Incomplete gamma function
+        sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$
+        and complementary incomplete gamma function
+        sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$
 
     sas_erf(x), sas_erfc(x):
@@ -207,4 +219,7 @@
 from numpy import pi, nan, inf
 from scipy.special import gamma as sas_gamma
+from scipy.special import gammaln as sas_gammaln
+from scipy.special import gammainc as sas_gammainc
+from scipy.special import gammaincc as sas_gammaincc
 from scipy.special import erf as sas_erf
 from scipy.special import erfc as sas_erfc
Index: setup.py
===================================================================
--- setup.py	(revision 1f991d6e4007acafaaf7edbc17c549a8b777d374)
+++ setup.py	(revision 783e76f9206beec611a37505457ee1b4561f150f)
@@ -29,4 +29,9 @@
                 return version[1:-1]
     raise RuntimeError("Could not read version from %s/__init__.py"%package)
+
+install_requires = ['numpy', 'scipy']
+
+if sys.platform=='win32' or sys.platform=='cygwin':
+    install_requires.append('tinycc')
 
 setup(
@@ -61,10 +66,9 @@
         'sasmodels': ['*.c', '*.cl'],
     },
-    install_requires=[
-    ],
+    install_requires=install_requires,
     extras_require={
+        'full': ['docutils', 'bumps', 'matplotlib'],
+        'server': ['bumps'],
         'OpenCL': ["pyopencl"],
-        'Bumps': ["bumps"],
-        'TinyCC': ["tinycc"],
     },
     build_requires=['setuptools'],
