Index: .travis.yml
===================================================================
--- .travis.yml	(revision 335271e7aa66a7275b4a987d7b763a60ff47493b)
+++ .travis.yml	(revision 5c36bf10d5206ab16fbf213971eae1414d26e8dc)
@@ -6,5 +6,5 @@
     env:
     - PY=2.7
-    - DEPLOY=True
+    #- DEPLOY=True
   - os: linux
     env:
Index: doc/guide/magnetism/magnetism.rst
===================================================================
--- doc/guide/magnetism/magnetism.rst	(revision 4f5afc935845263676247d906df9688597dfda82)
+++ doc/guide/magnetism/magnetism.rst	(revision befe905ba1f75b9bc31808657037e4ea891bdf28)
@@ -39,8 +39,8 @@
 
 .. math::
-    -- &= ((1-u_i)(1-u_f))^{1/4} \\
-    -+ &= ((1-u_i)(u_f))^{1/4} \\
-    +- &= ((u_i)(1-u_f))^{1/4} \\
-    ++ &= ((u_i)(u_f))^{1/4}
+    -- &= (1-u_i)(1-u_f) \\
+    -+ &= (1-u_i)(u_f) \\
+    +- &= (u_i)(1-u_f) \\
+    ++ &= (u_i)(u_f)
 
 Ideally the experiment would measure the pure spin states independently and
@@ -104,2 +104,3 @@
 | 2015-05-02 Steve King
 | 2017-11-15 Paul Kienzle
+| 2018-06-02 Adam Washington
Index: doc/guide/pd/polydispersity.rst
===================================================================
--- doc/guide/pd/polydispersity.rst	(revision 29afc50fdd56895a3c14d9c1e09b002d087c956a)
+++ doc/guide/pd/polydispersity.rst	(revision d712a0f4699e15d422e39b923be41331fe3922bf)
@@ -20,6 +20,11 @@
   P(q) = \text{scale} \langle F^* F \rangle / V + \text{background}
 
-where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an
-average over the size distribution.
+where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an 
+average over the size distribution $f(x; \bar x, \sigma)$, giving
+
+.. math::
+
+  P(q) = \frac{\text{scale}}{V} \int_\mathbb{R} 
+  f(x; \bar x, \sigma) F^2(q, x)\, dx + \text{background}
 
 Each distribution is characterized by a center value $\bar x$ or
@@ -41,5 +46,5 @@
 with larger values of $N_\sigma$ required for heavier tailed distributions.
 The scattering in general falls rapidly with $qr$ so the usual assumption
-that $G(r - 3\sigma_r)$ is tiny and therefore $f(r - 3\sigma_r)G(r - 3\sigma_r)$
+that $f(r - 3\sigma_r)$ is tiny and therefore $f(r - 3\sigma_r)f(r - 3\sigma_r)$
 will not contribute much to the average may not hold when particles are large.
 This, too, will require increasing $N_\sigma$.
@@ -63,4 +68,13 @@
 
 Additional distributions are under consideration.
+
+.. note:: In 2009 IUPAC decided to introduce the new term 'dispersity' to replace 
+           the term 'polydispersity' (see `Pure Appl. Chem., (2009), 81(2), 
+           351-353 <http://media.iupac.org/publications/pac/2009/pdf/8102x0351.pdf>`_ 
+           in order to make the terminology describing distributions of properties 
+           unambiguous. Throughout the SasView documentation we continue to use the 
+           term polydispersity because one of the consequences of the IUPAC change is 
+           that orientational polydispersity would not meet their new criteria (which 
+           requires dispersity to be dimensionless).
 
 Suggested Applications
Index: doc/guide/plugin.rst
===================================================================
--- doc/guide/plugin.rst	(revision 7e6bc45eddf06819ea47d06f521d0c27c582996e)
+++ doc/guide/plugin.rst	(revision f796469ab131e5fc2424ecdfcd201846be55c8e1)
@@ -822,5 +822,5 @@
 
         :code:`source = ["lib/Si.c", ...]`
-        (`Si.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/Si.c>`_)
+        (`Si.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_Si.c>`_)
 
     sas_3j1x_x(x):
Index: example/model_ellipsoid_hayter_msa.py
===================================================================
--- example/model_ellipsoid_hayter_msa.py	(revision 8a5f0216970a47f3e897a508054840132bd87909)
+++ example/model_ellipsoid_hayter_msa.py	(revision 9a999931f7b7ce4cd7a1ccb4cb1b1c86845e1065)
@@ -16,5 +16,5 @@
 
 # DEFINE THE MODEL
-kernel = load_model('ellipsoid*hayter_msa')
+kernel = load_model('ellipsoid@hayter_msa')
 
 pars = dict(scale=6.4, background=0.06, sld=0.33, sld_solvent=2.15, radius_polar=14.0,
Index: example/slurm_batch.py
===================================================================
--- example/slurm_batch.py	(revision 1d9998c981957f42dd41ba7220fbc68bdfba7539)
+++ example/slurm_batch.py	(revision 1d9998c981957f42dd41ba7220fbc68bdfba7539)
@@ -0,0 +1,241 @@
+#!/usr/bin/env python
+"""
+Submit a batch fit job to the slurm cluster.
+
+Given a model.py file defining a Bumps problem defined on a single data
+file, with the data file specified as a command line argument, run the
+bumps fit as a batch over a set of different datafiles independently.
+An example model is given in model_ellipsoid_hayter_msa.py, which fits 
+the data in 09319*.dat.
+
+To run the fit, use::
+
+    slurm_batch.py [--slurm_opts] model.py *.dat --store=T1 [--bumps_opt ...]
+
+For example::
+
+    slurm_batch.py model_ellipsoid_hayter_msa.py 09319*.dat --store=T1
+
+This creates the T1 subdirectory to hold the fit results and 
+prints the real command that is submitted, as well as the job id.
+
+The store directory T1 contains a copy of the model file and 
+all the data files.  The fit results for each file will be 
+in T1/##/*.  The file T1/files.dat contains the list 
+of "subdirectory filename" pairs indicating which ## directory 
+contains the resuls for which file.  Check for errors using::
+
+    cat T1/slurm*_1.out
+
+The following slurm options are used::
+
+    --array=1-#files     batch size comes from the file list
+    --gres=gpu:1         request a gpu for each fit
+    --job-name=model.py  use model file name for job name
+    --output=...         log into T1/slurm-job_##.out
+    --chdir=...          run fit from store directory
+    --time=2             time as number of hours (can override)
+
+To receive an email on job completion or failure, add the following
+slurm options before the model file::
+
+    --mail-type=END,FAIL --mail-user=user@mail.domain
+
+Bumps options are described at bumps.readthedocs.org, with the
+following set automatically::
+
+    --batch              run in batch mode, without output to .mon
+    --view=log           SAS fits want log plots
+    --time=2-0.1         slurm time minus 6 minutes for cleanup
+
+The --store and --resume options indicate the parent directory for
+the output.  These are modified to store the results in a separate
+subdirectory for each file.  Keep in mind that the fit is run from
+the store directory, so any files or modules referenced from the
+model file will need to use a full path to the original location.
+
+After submitting the job a job id will be printed to the console.
+You can check the status of the job using the usual slurm commands
+such as::
+
+    squeue
+
+or cancel the job using::
+
+    scancel jobid
+
+The slurm_batch program runs directly from the source tree for sasmodels, 
+and requires sasview, bumps and periodictable as sister directories
+accessible on the worker nodes.  You can link it into your bin directory
+using::
+
+     mkdir ~/bin
+     ln -s path/to/slurm_batch.py ~/bin
+
+or if you are a cluster administrator, into /usr/local/bin.
+"""
+
+# If called from command line, this submits a job to the slurm queue, with _this_ file
+# as the batch script. Before calling it on the worker node, slurm sets the
+# SLURM_ARRAY_TASK_ID to the current task so we can tell that we are running
+# as a worker and which file we should be working on.
+
+## SBATCH options as comments do not seem to work.  Maybe they neeed to be before 
+## the doc string?  For now they are hardcoded in the sbatch call in submit_job.
+
+import sys
+import os
+import tempfile
+import shutil
+
+DEFAULT_TIME_LIMIT = 2
+
+def split_args():
+    slurm_opts = []
+    bumps_opts = []
+    model_file = None
+    store = None
+    resume = None
+    data_files = []
+    time_limit = DEFAULT_TIME_LIMIT
+
+    # start with '-' arguments as slurm opts, then after
+    # the model file any '-' arguments are bumps opts.
+    opts = slurm_opts
+    for v in sys.argv[1:]:
+        if v.startswith('--store='):
+            store = os.path.realpath(os.path.abspath(v[8:]))
+        elif v.startswith('--resume='):
+            resume = os.path.realpath(os.path.abspath(v[9:]))
+        elif v.startswith('--time='):
+            time_limit = float(v[7:])
+        elif v[0] == '-':
+            opts.append(v)
+        elif model_file is None:
+            model_file = v
+            opts = bumps_opts
+        else:
+            data_files.append(v)
+
+
+    s = time_limit*3600
+    slurm_opts.append("--time=%d:%02d:%02d"%(s//3600, (s%3600)//60, s%60))
+    bumps_opts.append('--time=%f'%(time_limit - 0.1))  # 6 min to stop cleanly
+
+    return {
+        'slurm': slurm_opts, 
+        'model_file': model_file, 
+        'data_files': data_files, 
+        'store': store, 
+        'resume': resume,
+        'bumps': bumps_opts,
+    }
+
+def dirn(path, n):
+    path = os.path.realpath(os.path.abspath(path))
+    for _ in range(n):
+        path = os.path.dirname(path)
+    return path
+
+def submit_job():
+    # sbatch --array=1-5 ./slurm_batch.py model_ellipsoid_hayter_msa.py 09*.dat --store=T1 --fit=dream
+    opts = split_args()
+    store = opts['store']
+    model_file = opts['model_file']
+    data_files = opts['data_files']
+    bumps_opts = opts['bumps']
+    slurm_opts = opts['slurm']
+
+    # make sure the store directory exists and save the order of the files, as well
+    # as the model and the data files
+    if store is not None:
+        if not os.path.exists(store):
+            os.makedirs(store)
+
+        # save file order
+        with open(os.path.join(store, 'files.dat'), 'w') as fid:
+            for k, f in enumerate(data_files):
+                fid.write("%02d %s\n"%(k+1, f))
+
+        # Copy the model and data files to the root store directory
+        # Since bumps changes into the model directory prior to loading
+        # the datafiles, strip all leading paths from data and model and
+        # set the working directory for the job to the store directory.
+        model_copy = os.path.basename(model_file)
+        shutil.copy(model_file, os.path.join(store, model_copy))
+        data_copy = []
+        for f in data_files:
+            f_copy = os.path.basename(f)
+            shutil.copy(f, os.path.join(store, f_copy))
+            data_copy.append(f_copy) 
+
+        model_file = model_copy
+        data_files = data_copy
+
+
+    # build and run the command
+    SRC = dirn(__file__, 3) # __file__ is $SRC/sasmodels/example/slurm_batch.py
+    parts = [
+        "sbatch",
+        "--array=1-%d"%len(data_files),
+        "--gres=gpu:1",
+        "--job-name="+model_file,
+        ## since we are setting the current working directory, we don't need
+        ## to fiddle the slurm output files
+        "--output=%s/slurm-%%A_%%a.out"%store,
+        "--chdir=%s"%store,
+        ]
+    parts.extend(slurm_opts)
+    parts.append(__file__)
+    # Remember the source root so we can reconstruct the correct python path
+    # This is done after the model file so that it doesn't get interpreted
+    # as a slurm option.
+    parts.append("--source_root=%s"%SRC)
+    parts.append(model_file)
+    parts.extend(data_files)
+    parts.extend(bumps_opts)
+    #if store is not None:
+    #    parts.append("--store=" + store)
+    command = " ".join(parts)
+
+    print(command)
+    os.system(command)
+
+def run_task(task_id):
+    opts = split_args()
+
+    # Set environment put compiled sasmodels in user-specific temporary cache
+    # We need this because users don't have a home directory on the individual
+    # cluster nodes.
+    assert opts['slurm'][0].startswith('--source_root=')
+    SRC = opts['slurm'][0][14:]
+    PACKAGES = ("periodictable", "sasview/src", "bumps", "sasmodels")
+    os.environ['PYTHONPATH'] = ":".join(SRC+"/"+v for v in PACKAGES)
+    TMP = tempfile.gettempdir()
+    cache_path = os.path.join(TMP, os.environ['USER'], '.cache')
+    os.environ['SAS_DLL_PATH'] = cache_path
+    os.environ['XDG_CACHE_HOME'] = cache_path
+
+    #task_store = "%s/%02d"%(opts['store'], task_id)
+    task_store = "%02d"%task_id
+    parts = [
+       "python", os.path.join(SRC, "bumps", "run.py"), "--batch",
+       "--view=log",
+       opts['model_file'],
+       opts['data_files'][task_id-1],
+       ]
+    parts.extend(opts['bumps'])
+    parts.append('--store='+task_store)
+    if opts['resume'] is not None:
+        parts.append('--resume='+os.path.join(opts['resume'], task_store))
+    command = " ".join(parts)
+    print(os.getcwd() + "$ " + command)
+    os.system(command)
+
+
+task_id = int(os.environ.get('SLURM_ARRAY_TASK_ID', -1))
+if task_id == -1:
+    submit_job()
+else:
+    run_task(task_id)
+
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 5770493cd59e15962cc1fb80022d0920b561cf63)
+++ sasmodels/compare.py	(revision af7a97c2acc05ec1c5d537fa21912dab44ff661a)
@@ -107,8 +107,17 @@
     -title="note" adds note to the plot title, after the model name
     -weights shows weights plots for the polydisperse parameters
+    -profile shows the sld profile if the model has a plottable sld profile
 
     === output options ===
     -edit starts the parameter explorer
     -help/-html shows the model docs instead of running the model
+
+    === environment variables ===
+    -DSAS_MODELPATH=path sets directory containing custom models
+    -DSAS_OPENCL=vendor:device|none sets the target OpenCL device
+    -DXDG_CACHE_HOME=~/.cache sets the pyopencl cache root (linux only)
+    -DSAS_COMPILER=tinycc|msvc|mingw|unix sets the DLL compiler
+    -DSAS_OPENMP=1 turns on OpenMP for the DLLs
+    -DSAS_DLL_PATH=path sets the path to the compiled modules
 
 The interpretation of quad precision depends on architecture, and may
@@ -645,5 +654,5 @@
 
 def make_data(opts):
-    # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray]
+    # type: (Dict[str, Any], float) -> Tuple[Data, np.ndarray]
     """
     Generate an empty dataset, used with the model to set Q points
@@ -667,5 +676,8 @@
         if opts['zero']:
             q = np.hstack((0, q))
-        data = empty_data1D(q, resolution=res)
+        # TODO: provide command line control of lambda and Delta lambda/lambda
+        #L, dLoL = 5, 0.14/np.sqrt(6)  # wavelength and 14% triangular FWHM
+        L, dLoL = 0, 0
+        data = empty_data1D(q, resolution=res, L=L, dL=L*dLoL)
         index = slice(None, None)
     return data, index
@@ -772,4 +784,21 @@
             dim = base._kernel.dim
             plot_weights(model_info, get_mesh(model_info, base_pars, dim=dim))
+        if opts['show_profile']:
+            import pylab
+            base, comp = opts['engines']
+            base_pars, comp_pars = opts['pars']
+            have_base = base._kernel.info.profile is not None
+            have_comp = (
+                comp is not None
+                and comp._kernel.info.profile is not None
+                and base_pars != comp_pars
+            )
+            if have_base or have_comp:
+                pylab.figure()
+            if have_base:
+                plot_profile(base._kernel.info, **base_pars)
+            if have_comp:
+                plot_profile(comp._kernel.info, label='comp', **comp_pars)
+                pylab.legend()
     if opts['plot']:
         import matplotlib.pyplot as plt
@@ -777,4 +806,33 @@
     return limits
 
+def plot_profile(model_info, label='base', **args):
+    # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None
+    """
+    Plot the profile returned by the model profile method.
+
+    *model_info* defines model parameters, etc.
+
+    *mesh* is a list of tuples containing (*value*, *dispersity*, *weights*)
+    for each parameter, where (*dispersity*, *weights*) pairs are the
+    distributions to be plotted.
+    """
+    import pylab
+
+    args = dict((k, v) for k, v in args.items()
+                if "_pd" not in k
+                   and ":" not in k
+                   and k not in ("background", "scale", "theta", "phi", "psi"))
+    args = args.copy()
+
+    args.pop('scale', 1.)
+    args.pop('background', 0.)
+    z, rho = model_info.profile(**args)
+    #pylab.interactive(True)
+    pylab.plot(z, rho, '-', label=label)
+    pylab.grid(True)
+    #pylab.show()
+
+
+
 def run_models(opts, verbose=False):
     # type: (Dict[str, Any]) -> Dict[str, Any]
@@ -786,5 +844,5 @@
     base_n, comp_n = opts['count']
     base_pars, comp_pars = opts['pars']
-    data = opts['data']
+    base_data, comp_data = opts['data']
 
     comparison = comp is not None
@@ -800,5 +858,5 @@
             print("%s t=%.2f ms, intensity=%.0f"
                   % (base.engine, base_time, base_value.sum()))
-        _show_invalid(data, base_value)
+        _show_invalid(base_data, base_value)
     except ImportError:
         traceback.print_exc()
@@ -812,5 +870,5 @@
                 print("%s t=%.2f ms, intensity=%.0f"
                       % (comp.engine, comp_time, comp_value.sum()))
-            _show_invalid(data, comp_value)
+            _show_invalid(base_data, comp_value)
         except ImportError:
             traceback.print_exc()
@@ -866,9 +924,10 @@
     have_base, have_comp = (base_value is not None), (comp_value is not None)
     base, comp = opts['engines']
-    data = opts['data']
+    base_data, comp_data = opts['data']
     use_data = (opts['datafile'] is not None) and (have_base ^ have_comp)
 
     # Plot if requested
     view = opts['view']
+    #view = 'log'
     if limits is None:
         vmin, vmax = np.inf, -np.inf
@@ -884,5 +943,5 @@
         if have_comp:
             plt.subplot(131)
-        plot_theory(data, base_value, view=view, use_data=use_data, limits=limits)
+        plot_theory(base_data, base_value, view=view, use_data=use_data, limits=limits)
         plt.title("%s t=%.2f ms"%(base.engine, base_time))
         #cbar_title = "log I"
@@ -891,6 +950,6 @@
             plt.subplot(132)
         if not opts['is2d'] and have_base:
-            plot_theory(data, base_value, view=view, use_data=use_data, limits=limits)
-        plot_theory(data, comp_value, view=view, use_data=use_data, limits=limits)
+            plot_theory(comp_data, base_value, view=view, use_data=use_data, limits=limits)
+        plot_theory(comp_data, comp_value, view=view, use_data=use_data, limits=limits)
         plt.title("%s t=%.2f ms"%(comp.engine, comp_time))
         #cbar_title = "log I"
@@ -908,5 +967,8 @@
             err[err > cutoff] = cutoff
         #err,errstr = base/comp,"ratio"
-        plot_theory(data, None, resid=err, view=errview, use_data=use_data)
+        # Note: base_data only since base and comp have same q values (though
+        # perhaps different resolution), and we are plotting the difference
+        # at each q
+        plot_theory(base_data, None, resid=err, view=errview, use_data=use_data)
         plt.xscale('log' if view == 'log' and not opts['is2d'] else 'linear')
         plt.legend(['P%d'%(k+1) for k in range(setnum+1)], loc='best')
@@ -942,5 +1004,6 @@
 OPTIONS = [
     # Plotting
-    'plot', 'noplot', 'weights',
+    'plot', 'noplot',
+    'weights', 'profile',
     'linear', 'log', 'q4',
     'rel', 'abs',
@@ -1058,6 +1121,7 @@
 
     invalid = [o[1:] for o in flags
-               if o[1:] not in NAME_OPTIONS
-               and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)]
+               if not (o[1:] in NAME_OPTIONS
+                       or any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)
+                       or o.startswith('-D'))]
     if invalid:
         print("Invalid options: %s"%(", ".join(invalid)))
@@ -1075,5 +1139,5 @@
         'qmax'      : 0.05,
         'nq'        : 128,
-        'res'       : 0.0,
+        'res'       : '0.0',
         'noise'     : 0.0,
         'accuracy'  : 'Low',
@@ -1096,4 +1160,5 @@
         'count'     : '1',
         'show_weights' : False,
+        'show_profile' : False,
         'sphere'    : 0,
         'ngauss'    : '0',
@@ -1115,5 +1180,5 @@
         elif arg.startswith('-q='):
             opts['qmin'], opts['qmax'] = [float(v) for v in arg[3:].split(':')]
-        elif arg.startswith('-res='):      opts['res'] = float(arg[5:])
+        elif arg.startswith('-res='):      opts['res'] = arg[5:]
         elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:])
         elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:])
@@ -1156,6 +1221,10 @@
         elif arg == '-default': opts['use_demo'] = False
         elif arg == '-weights': opts['show_weights'] = True
+        elif arg == '-profile': opts['show_profile'] = True
         elif arg == '-html':    opts['html'] = True
         elif arg == '-help':    opts['html'] = True
+        elif arg.startswith('-D'):
+            var, val = arg[2:].split('=')
+            os.environ[var] = val
     # pylint: enable=bad-whitespace,C0321
 
@@ -1173,8 +1242,4 @@
     if opts['qmin'] is None:
         opts['qmin'] = 0.001*opts['qmax']
-    if opts['datafile'] is not None:
-        data = load_data(os.path.expanduser(opts['datafile']))
-    else:
-        data, _ = make_data(opts)
 
     comparison = any(PAR_SPLIT in v for v in values)
@@ -1216,8 +1281,29 @@
         opts['cutoff'] = [float(opts['cutoff'])]*2
 
-    base = make_engine(model_info[0], data, opts['engine'][0],
+    if PAR_SPLIT in opts['res']:
+        opts['res'] = [float(k) for k in opts['res'].split(PAR_SPLIT, 2)]
+        comparison = True
+    else:
+        opts['res'] = [float(opts['res'])]*2
+
+    if opts['datafile'] is not None:
+        data = load_data(os.path.expanduser(opts['datafile']))
+    else:
+        # Hack around the fact that make_data doesn't take a pair of resolutions
+        res = opts['res']
+        opts['res'] = res[0]
+        data0, _ = make_data(opts)
+        if res[0] != res[1]:
+            opts['res'] = res[1]
+            data1, _ = make_data(opts)
+        else:
+            data1 = data0
+        opts['res'] = res
+        data = data0, data1
+
+    base = make_engine(model_info[0], data[0], opts['engine'][0],
                        opts['cutoff'][0], opts['ngauss'][0])
     if comparison:
-        comp = make_engine(model_info[1], data, opts['engine'][1],
+        comp = make_engine(model_info[1], data[1], opts['engine'][1],
                            opts['cutoff'][1], opts['ngauss'][1])
     else:
@@ -1376,5 +1462,5 @@
     path = os.path.dirname(info.filename)
     url = "file://" + path.replace("\\", "/")[2:] + "/"
-    rst2html.view_html_qtapp(html, url)
+    rst2html.view_html_wxapp(html, url)
 
 def explore(opts):
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision 3221de0ed9b3480163246ceb0503432126bd5134)
+++ sasmodels/core.py	(revision 4341dd4940a5037121414aaef66772bc0bd76640)
@@ -37,6 +37,6 @@
 if CUSTOM_MODEL_PATH == "":
     CUSTOM_MODEL_PATH = joinpath(os.path.expanduser("~"), ".sasmodels", "custom_models")
-    if not os.path.isdir(CUSTOM_MODEL_PATH):
-        os.makedirs(CUSTOM_MODEL_PATH)
+    #if not os.path.isdir(CUSTOM_MODEL_PATH):
+    #    os.makedirs(CUSTOM_MODEL_PATH)
 
 # TODO: refactor composite model support
Index: sasmodels/data.py
===================================================================
--- sasmodels/data.py	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
+++ sasmodels/data.py	(revision 1a8c11c9d9d7223aecf2c56d6742e33302bffda3)
@@ -36,4 +36,5 @@
 
 import numpy as np  # type: ignore
+from numpy import sqrt, sin, cos, pi
 
 # pylint: disable=unused-import
@@ -301,10 +302,13 @@
 
 
-def empty_data1D(q, resolution=0.0):
+def empty_data1D(q, resolution=0.0, L=0., dL=0.):
     # type: (np.ndarray, float) -> Data1D
-    """
+    r"""
     Create empty 1D data using the given *q* as the x value.
 
-    *resolution* dq/q defaults to 5%.
+    rms *resolution* $\Delta q/q$ defaults to 0%.  If wavelength *L* and rms
+    wavelength divergence *dL* are defined, then *resolution* defines
+    rms $\Delta \theta/\theta$ for the lowest *q*, with $\theta$ derived from
+    $q = 4\pi/\lambda \sin(\theta)$.
     """
 
@@ -313,5 +317,18 @@
     Iq, dIq = None, None
     q = np.asarray(q)
-    data = Data1D(q, Iq, dx=resolution * q, dy=dIq)
+    if L != 0 and resolution != 0:
+        theta = np.arcsin(q*L/(4*pi))
+        dtheta = theta[0]*resolution
+        ## Solving Gaussian error propagation from
+        ##   Dq^2 = (dq/dL)^2 DL^2 + (dq/dtheta)^2 Dtheta^2
+        ## gives
+        ##   (Dq/q)^2 = (DL/L)**2 + (Dtheta/tan(theta))**2
+        ## Take the square root and multiply by q, giving
+        ##   Dq = (4*pi/L) * sqrt((sin(theta)*dL/L)**2 + (cos(theta)*dtheta)**2)
+        dq = (4*pi/L) * sqrt((sin(theta)*dL/L)**2 + (cos(theta)*dtheta)**2)
+    else:
+        dq = resolution * q
+
+    data = Data1D(q, Iq, dx=dq, dy=dIq)
     data.filename = "fake data"
     return data
@@ -486,9 +503,11 @@
             # Note: masks merge, so any masked theory points will stay masked,
             # and the data mask will be added to it.
-            mtheory = masked_array(theory, data.mask.copy())
+            #mtheory = masked_array(theory, data.mask.copy())
+            theory_x = data.x[~data.mask]
+            mtheory = masked_array(theory)
             mtheory[~np.isfinite(mtheory)] = masked
             if view is 'log':
                 mtheory[mtheory <= 0] = masked
-            plt.plot(data.x, scale*mtheory, '-')
+            plt.plot(theory_x, scale*mtheory, '-')
             all_positive = all_positive and (mtheory > 0).all()
             some_present = some_present or (mtheory.count() > 0)
@@ -526,5 +545,6 @@
 
     if use_resid:
-        mresid = masked_array(resid, data.mask.copy())
+        theory_x = data.x[~data.mask]
+        mresid = masked_array(resid)
         mresid[~np.isfinite(mresid)] = masked
         some_present = (mresid.count() > 0)
@@ -532,5 +552,5 @@
         if num_plots > 1:
             plt.subplot(1, num_plots, use_calc + 2)
-        plt.plot(data.x, mresid, '.')
+        plt.plot(theory_x, mresid, '.')
         plt.xlabel("$q$/A$^{-1}$")
         plt.ylabel('residuals')
Index: sasmodels/direct_model.py
===================================================================
--- sasmodels/direct_model.py	(revision d18d6ddf640c8cff5c76cfac900c1940989612ca)
+++ sasmodels/direct_model.py	(revision 1a8c11c9d9d7223aecf2c56d6742e33302bffda3)
@@ -250,5 +250,5 @@
             qmax = getattr(data, 'qmax', np.inf)
             accuracy = getattr(data, 'accuracy', 'Low')
-            index = ~data.mask & (q >= qmin) & (q <= qmax)
+            index = (data.mask == 0) & (q >= qmin) & (q <= qmax)
             if data.data is not None:
                 index &= ~np.isnan(data.data)
@@ -263,4 +263,7 @@
         elif self.data_type == 'Iq':
             index = (data.x >= data.qmin) & (data.x <= data.qmax)
+            mask = getattr(data, 'mask', None)
+            if mask is not None:
+                index &= (mask == 0)
             if data.y is not None:
                 index &= ~np.isnan(data.y)
Index: sasmodels/jitter.py
===================================================================
--- sasmodels/jitter.py	(revision b3703f51bafe67afb8acefe986ad7466720f755b)
+++ sasmodels/jitter.py	(revision 1198f905febb85403d3b86dca3d2bec81af13029)
@@ -774,5 +774,5 @@
         # set small jitter as 0 if multiple pd dims
         dims = sum(v > 0 for v in jitter)
-        limit = [0, 0, 0.5, 5][dims]
+        limit = [0, 0.5, 5][dims]
         jitter = [0 if v < limit else v for v in jitter]
         axes.cla()
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision dc6f601ca60f7d630a484562e056a94ed48bae82)
+++ sasmodels/kernel_iq.c	(revision 7c35fda803d68561ee8433f7c9355dcc9678049e)
@@ -83,11 +83,19 @@
   in_spin = clip(in_spin, 0.0, 1.0);
   out_spin = clip(out_spin, 0.0, 1.0);
-  // Note: sasview 3.1 scaled all slds by sqrt(weight) and assumed that
+  // Previous version of this function took the square root of the weights,
+  // under the assumption that 
+  //
   //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...)
-  // which is likely to be the case for simple models.
-  weight[0] = sqrt((1.0-in_spin) * (1.0-out_spin)); // dd
-  weight[1] = sqrt((1.0-in_spin) * out_spin);       // du.real
-  weight[2] = sqrt(in_spin * (1.0-out_spin));       // ud.real
-  weight[3] = sqrt(in_spin * out_spin);             // uu
+  //
+  // However, since the weights are applied to the final intensity and
+  // are not interned inside the I(q) function, we want the full
+  // weight and not the square root.  Any function using
+  // set_spin_weights as part of calculating an amplitude will need to
+  // manually take that square root, but there is currently no such
+  // function.
+  weight[0] = (1.0-in_spin) * (1.0-out_spin); // dd
+  weight[1] = (1.0-in_spin) * out_spin;       // du
+  weight[2] = in_spin * (1.0-out_spin);       // ud
+  weight[3] = in_spin * out_spin;             // uu
   weight[4] = weight[1]; // du.imag
   weight[5] = weight[2]; // ud.imag
Index: sasmodels/kerneldll.py
===================================================================
--- sasmodels/kerneldll.py	(revision bf94e6e96340e0e655c2e0a16ff33bde3736ac5f)
+++ sasmodels/kerneldll.py	(revision 1a3559f14a5729a3c0f491a6669b4bd4429697c7)
@@ -99,4 +99,10 @@
     pass
 # pylint: enable=unused-import
+
+if "SAS_DLL_PATH" in os.environ:
+    SAS_DLL_PATH = os.environ["SAS_DLL_PATH"]
+else:
+    # Assume the default location of module DLLs is in .sasmodels/compiled_models.
+    SAS_DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models")
 
 if "SAS_COMPILER" in os.environ:
@@ -123,5 +129,8 @@
     # add openmp support if not running on a mac
     if sys.platform != "darwin":
-        CC.append("-fopenmp")
+        # OpenMP seems to be broken on gcc 5.4.0 (ubuntu 16.04.9)
+        # Shut it off for all unix until we can investigate.
+        #CC.append("-fopenmp")
+        pass
     def compile_command(source, output):
         """unix compiler command"""
@@ -158,7 +167,4 @@
         return CC + [source, "-o", output, "-lm"]
 
-# Assume the default location of module DLLs is in .sasmodels/compiled_models.
-DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models")
-
 ALLOW_SINGLE_PRECISION_DLLS = True
 
@@ -197,5 +203,5 @@
         return path
 
-    return joinpath(DLL_PATH, basename)
+    return joinpath(SAS_DLL_PATH, basename)
 
 
@@ -206,5 +212,5 @@
     exist yet if it hasn't been compiled.
     """
-    return os.path.join(DLL_PATH, dll_name(model_info, dtype))
+    return os.path.join(SAS_DLL_PATH, dll_name(model_info, dtype))
 
 
@@ -225,5 +231,6 @@
     models are not allowed as DLLs.
 
-    Set *sasmodels.kerneldll.DLL_PATH* to the compiled dll output path.
+    Set *sasmodels.kerneldll.SAS_DLL_PATH* to the compiled dll output path.
+    Alternatively, set the environment variable *SAS_DLL_PATH*.
     The default is in ~/.sasmodels/compiled_models.
     """
@@ -244,6 +251,6 @@
     if need_recompile:
         # Make sure the DLL path exists
-        if not os.path.exists(DLL_PATH):
-            os.makedirs(DLL_PATH)
+        if not os.path.exists(SAS_DLL_PATH):
+            os.makedirs(SAS_DLL_PATH)
         basename = splitext(os.path.basename(dll))[0] + "_"
         system_fd, filename = tempfile.mkstemp(suffix=".c", prefix=basename)
Index: sasmodels/models/core_shell_parallelepiped.c
===================================================================
--- sasmodels/models/core_shell_parallelepiped.c	(revision e0772316ffb087990c2e1b341c744594e3087097)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision dbf1a60059f2daf9ef240db7a7fbf91cf5de9b9e)
@@ -59,17 +59,19 @@
 
     // outer integral (with gauss points), integration limits = 0, 1
+    // substitute d_cos_alpha for sin_alpha d_alpha
     double outer_sum = 0; //initialize integral
     for( int i=0; i<GAUSS_N; i++) {
         const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 );
         const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha);
-
-        // inner integral (with gauss points), integration limits = 0, pi/2
         const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q);
         const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q);
+
+        // inner integral (with gauss points), integration limits = 0, 1
+        // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta)
         double inner_sum = 0.0;
         for(int j=0; j<GAUSS_N; j++) {
-            const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 );
+            const double u = 0.5 * ( GAUSS_Z[j] + 1.0 );
             double sin_beta, cos_beta;
-            SINCOS(M_PI_2*beta, sin_beta, cos_beta);
+            SINCOS(M_PI_2*u, sin_beta, cos_beta);
             const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta);
             const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta);
@@ -91,8 +93,10 @@
             inner_sum += GAUSS_W[j] * f * f;
         }
+        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5
         inner_sum *= 0.5;
         // now sum up the outer integral
         outer_sum += GAUSS_W[i] * inner_sum;
     }
+    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5
     outer_sum *= 0.5;
 
Index: sasmodels/models/core_shell_parallelepiped.py
===================================================================
--- sasmodels/models/core_shell_parallelepiped.py	(revision 97be87700f140bc9317b791ad4691cc31fa30981)
+++ sasmodels/models/core_shell_parallelepiped.py	(revision f89ec96583ec917dee91d159ad08cdc6c4222b10)
@@ -4,6 +4,12 @@
 
 Calculates the form factor for a rectangular solid with a core-shell structure.
-The thickness and the scattering length density of the shell or
-"rim" can be different on each (pair) of faces.
+The thickness and the scattering length density of the shell or "rim" can be
+different on each (pair) of faces. The three dimensions of the core of the
+parallelepiped (strictly here a cuboid) may be given in *any* size order as
+long as the particles are randomly oriented (i.e. take on all possible
+orientations see notes on 2D below). To avoid multiple fit solutions,
+especially with Monte-Carlo fit methods, it may be advisable to restrict their
+ranges. There may be a number of closely similar "best fits", so some trial and
+error, or fixing of some dimensions at expected values, may help.
 
 The form factor is normalized by the particle volume $V$ such that
@@ -11,47 +17,60 @@
 .. math::
 
-    I(q) = \text{scale}\frac{\langle f^2 \rangle}{V} + \text{background}
+    I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle
+    + \text{background}
 
 where $\langle \ldots \rangle$ is an average over all possible orientations
-of the rectangular solid.
-
-The function calculated is the form factor of the rectangular solid below.
-The core of the solid is defined by the dimensions $A$, $B$, $C$ such that
-$A < B < C$.
-
-.. image:: img/core_shell_parallelepiped_geometry.jpg
+of the rectangular solid, and the usual $\Delta \rho^2 \ V^2$ term cannot be
+pulled out of the form factor term due to the multiple slds in the model.
+
+The core of the solid is defined by the dimensions $A$, $B$, $C$ here shown
+such that $A < B < C$.
+
+.. figure:: img/parallelepiped_geometry.jpg
+
+   Core of the core shell parallelepiped with the corresponding definition
+   of sides.
+
 
 There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension
 (on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$
-$(=t_C)$ faces. The projection in the $AB$ plane is then
-
-.. image:: img/core_shell_parallelepiped_projection.jpg
-
-The volume of the solid is
+$(=t_C)$ faces. The projection in the $AB$ plane is
+
+.. figure:: img/core_shell_parallelepiped_projection.jpg
+
+   AB cut through the core-shell parallelipiped showing the cross secion of
+   four of the six shell slabs. As can be seen, this model leaves **"gaps"**
+   at the corners of the solid.
+
+
+The total volume of the solid is thus given as
 
 .. math::
 
     V = ABC + 2t_ABC + 2t_BAC + 2t_CAB
-
-**meaning that there are "gaps" at the corners of the solid.**
 
 The intensity calculated follows the :ref:`parallelepiped` model, with the
 core-shell intensity being calculated as the square of the sum of the
-amplitudes of the core and the slabs on the edges.
-
-the scattering amplitude is computed for a particular orientation of the
-core-shell parallelepiped with respect to the scattering vector and then
-averaged over all possible orientations, where $\alpha$ is the angle between
-the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is
-the angle between projection of the particle in the $xy$ detector plane
-and the $y$ axis.
-
-.. math::
-
-    F(Q)
+amplitudes of the core and the slabs on the edges. The scattering amplitude is
+computed for a particular orientation of the core-shell parallelepiped with
+respect to the scattering vector and then averaged over all possible
+orientations, where $\alpha$ is the angle between the $z$ axis and the $C$ axis
+of the parallelepiped, and $\beta$ is the angle between the projection of the
+particle in the $xy$ detector plane and the $y$ axis.
+
+.. math::
+
+    P(q)=\frac {\int_{0}^{\pi/2}\int_{0}^{\pi/2}F^2(q,\alpha,\beta) \ sin\alpha
+    \ d\alpha \ d\beta} {\int_{0}^{\pi/2} \ sin\alpha \ d\alpha \ d\beta}
+
+and
+
+.. math::
+
+    F(q,\alpha,\beta)
     &= (\rho_\text{core}-\rho_\text{solvent})
        S(Q_A, A) S(Q_B, B) S(Q_C, C) \\
     &+ (\rho_\text{A}-\rho_\text{solvent})
-        \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\
+        \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\
     &+ (\rho_\text{B}-\rho_\text{solvent})
         S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\
@@ -63,5 +82,5 @@
 .. math::
 
-    S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L}
+    S(Q_X, L) = L \frac{\sin (\tfrac{1}{2} Q_X L)}{\tfrac{1}{2} Q_X L}
 
 and
@@ -69,50 +88,65 @@
 .. math::
 
-    Q_A &= \sin\alpha \sin\beta \\
-    Q_B &= \sin\alpha \cos\beta \\
-    Q_C &= \cos\alpha
+    Q_A &= q \sin\alpha \sin\beta \\
+    Q_B &= q \sin\alpha \cos\beta \\
+    Q_C &= q \cos\alpha
 
 
 where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$
-are the scattering length of the parallelepiped core, and the rectangular
+are the scattering lengths of the parallelepiped core, and the rectangular
 slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$
 is the scattering length of the solvent.
 
+.. note:: 
+
+   the code actually implements two substitutions: $d(cos\alpha)$ is
+   substituted for -$sin\alpha \ d\alpha$ (note that in the
+   :ref:`parallelepiped` code this is explicitly implemented with
+   $\sigma = cos\alpha$), and $\beta$ is set to $\beta = u \pi/2$ so that
+   $du = \pi/2 \ d\beta$.  Thus both integrals go from 0 to 1 rather than 0
+   to $\pi/2$.
+
 FITTING NOTES
 ~~~~~~~~~~~~~
 
-If the scale is set equal to the particle volume fraction, $\phi$, the returned
-value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However,
-**no interparticle interference effects are included in this calculation.**
-
-There are many parameters in this model. Hold as many fixed as possible with
-known values, or you will certainly end up at a solution that is unphysical.
-
-The returned value is in units of |cm^-1|, on absolute scale.
-
-NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated
-based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$
-and length $(C+2t_C)$ values, after appropriately sorting the three dimensions
-to give an oblate or prolate particle, to give an effective radius,
-for $S(Q)$ when $P(Q) * S(Q)$ is applied.
-
-For 2d data the orientation of the particle is required, described using
-angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further
-details of the calculation and angular dispersions see :ref:`orientation`.
-The angle $\Psi$ is the rotational angle around the *long_c* axis. For example,
-$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector.
-
-For 2d, constraints must be applied during fitting to ensure that the
-inequality $A < B < C$ is not violated, and hence the correct definition
-of angles is preserved. The calculation will not report an error,
-but the results may be not correct.
+#. There are many parameters in this model. Hold as many fixed as possible with
+   known values, or you will certainly end up at a solution that is unphysical.
+
+#. The 2nd virial coefficient of the core_shell_parallelepiped is calculated
+   based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$
+   and length $(C+2t_C)$ values, after appropriately sorting the three
+   dimensions to give an oblate or prolate particle, to give an effective radius
+   for $S(q)$ when $P(q) * S(q)$ is applied.
+
+#. For 2d data the orientation of the particle is required, described using
+   angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, where $\theta$
+   and $\phi$ define the orientation of the director in the laboratry reference
+   frame of the beam direction ($z$) and detector plane ($x-y$ plane), while
+   the angle $\Psi$ is effectively the rotational angle around the particle
+   $C$ axis. For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the
+   $B$ axis oriented parallel to the y-axis of the detector with $A$ along
+   the x-axis. For other $\theta$, $\phi$ values, the order of rotations
+   matters. In particular, the parallelepiped must first be rotated $\theta$
+   degrees in the $x-z$ plane before rotating $\phi$ degrees around the $z$
+   axis (in the $x-y$ plane). Applying orientational distribution to the
+   particle orientation (i.e  `jitter` to one or more of these angles) can get
+   more confusing as `jitter` is defined **NOT** with respect to the laboratory
+   frame but the particle reference frame. It is thus highly recmmended to
+   read :ref:`orientation` for further details of the calculation and angular
+   dispersions.
+
+.. note:: For 2d, constraints must be applied during fitting to ensure that the
+   order of sides chosen is not altered, and hence that the correct definition
+   of angles is preserved. For the default choice shown here, that means
+   ensuring that the inequality $A < B < C$ is not violated,  The calculation
+   will not report an error, but the results may be not correct.
 
 .. figure:: img/parallelepiped_angle_definition.png
 
     Definition of the angles for oriented core-shell parallelepipeds.
-    Note that rotation $\theta$, initially in the $xz$ plane, is carried
+    Note that rotation $\theta$, initially in the $x-z$ plane, is carried
     out first, then rotation $\phi$ about the $z$ axis, finally rotation
-    $\Psi$ is now around the axis of the cylinder. The neutron or X-ray
-    beam is along the $z$ axis.
+    $\Psi$ is now around the $C$ axis of the particle. The neutron or X-ray
+    beam is along the $z$ axis and the detecotr defines the $x-y$ plane.
 
 .. figure:: img/parallelepiped_angle_projection.png
@@ -120,4 +154,13 @@
     Examples of the angles for oriented core-shell parallelepipeds against the
     detector plane.
+
+
+Validation
+----------
+
+Cross-checked against hollow rectangular prism and rectangular prism for equal
+thickness overlapping sides, and by Monte Carlo sampling of points within the
+shape for non-uniform, non-overlapping sides.
+
 
 References
@@ -135,9 +178,8 @@
 
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
-* **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016
+* **Converted to sasmodels by:** Miguel Gonzalez **Date:** February 26, 2016
 * **Last Modified by:** Paul Kienzle **Date:** October 17, 2017
-* Cross-checked against hollow rectangular prism and rectangular prism for
-  equal thickness overlapping sides, and by Monte Carlo sampling of points
-  within the shape for non-uniform, non-overlapping sides.
+* **Last Reviewed by:** Paul Butler **Date:** May 24, 2018 - documentation
+  updated
 """
 
Index: sasmodels/models/core_shell_sphere.py
===================================================================
--- sasmodels/models/core_shell_sphere.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/core_shell_sphere.py	(revision dc76240c1dc8dcc0235d0432f3a909bea6e0c6a2)
@@ -21,5 +21,5 @@
 .. math::
 
-    F^2(q) = \frac{3}{V_s}\left[
+    F(q) = \frac{3}{V_s}\left[
        V_c(\rho_c-\rho_s)\frac{\sin(qr_c)-qr_c\cos(qr_c)}{(qr_c)^3} +
        V_s(\rho_s-\rho_\text{solv})\frac{\sin(qr_s)-qr_s\cos(qr_s)}{(qr_s)^3}
Index: sasmodels/models/guinier.py
===================================================================
--- sasmodels/models/guinier.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/guinier.py	(revision c9fc87335cd86f10c8c9e3654927fdeeba6a70af)
@@ -7,5 +7,5 @@
 .. math::
 
-    I(q) = \text{scale} \cdot \exp{\left[ \frac{-Q^2R_g^2}{3} \right]}
+    I(q) = \text{scale} \cdot \exp{\left[ \frac{-Q^2 R_g^2 }{3} \right]}
             + \text{background}
 
@@ -19,4 +19,27 @@
 
 .. math:: q = \sqrt{q_x^2 + q_y^2}
+
+In scattering, the radius of gyration $R_g$ quantifies the objects's
+distribution of SLD (not mass density, as in mechanics) from the objects's
+SLD centre of mass. It is defined by
+
+.. math:: R_g^2 = \frac{\sum_i\rho_i\left(r_i-r_0\right)^2}{\sum_i\rho_i}
+
+where $r_0$ denotes the object's SLD centre of mass and $\rho_i$ is the SLD at
+a point $i$.
+
+Notice that $R_g^2$ may be negative (since SLD can be negative), which happens
+when a form factor $P(Q)$ is increasing with $Q$ rather than decreasing. This
+can occur for core/shell particles, hollow particles, or for composite
+particles with domains of different SLDs in a solvent with an SLD close to the
+average match point. (Alternatively, this might be regarded as there being an
+internal inter-domain "structure factor" within a single particle which gives
+rise to a peak in the scattering).
+
+To specify a negative value of $R_g^2$ in SasView, simply give $R_g$ a negative
+value ($R_g^2$ will be evaluated as $R_g |R_g|$). Note that the physical radius 
+of gyration, of the exterior of the particle, will still be large and positive. 
+It is only the apparent size from the small $Q$ data that will give a small or 
+negative value of $R_g^2$.
 
 References
@@ -42,8 +65,8 @@
 
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["rg", "Ang", 60.0, [0, inf], "", "Radius of Gyration"]]
+parameters = [["rg", "Ang", 60.0, [-inf, inf], "", "Radius of Gyration"]]
 
 Iq = """
-    double exponent = rg*rg*q*q/3.0;
+    double exponent = fabs(rg)*rg*q*q/3.0;
     double value = exp(-exponent);
     return value;
@@ -66,5 +89,5 @@
 
 # parameters for demo
-demo = dict(scale=1.0, rg=60.0)
+demo = dict(scale=1.0,  background=0.001, rg=60.0 )
 
 # parameters for unit tests
Index: sasmodels/models/mass_surface_fractal.py
===================================================================
--- sasmodels/models/mass_surface_fractal.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/mass_surface_fractal.py	(revision 7994359d6c078cf35c29184c8eb57657f44ef374)
@@ -39,14 +39,25 @@
     The surface ( $D_s$ ) and mass ( $D_m$ ) fractal dimensions are only
     valid if $0 < surface\_dim < 6$ , $0 < mass\_dim < 6$ , and
-    $(surface\_dim + mass\_dim ) < 6$ .
-
+    $(surface\_dim + mass\_dim ) < 6$ . 
+    Older versions of sasview may have the default primary particle radius
+    larger than the cluster radius, this was an error, also present in the 
+    Schmidt review paper below. The primary particle should be the smaller 
+    as described in the original Hurd et.al. who also point out that 
+    polydispersity in the primary particle sizes may affect their 
+    apparent surface fractal dimension.
+    
 
 References
 ----------
 
-P Schmidt, *J Appl. Cryst.*, 24 (1991) 414-435 Equation(19)
+.. [#] P Schmidt, *J Appl. Cryst.*, 24 (1991) 414-435 Equation(19)
+.. [#] A J Hurd, D W Schaefer, J E Martin, *Phys. Rev. A*,
+   35 (1987) 2361-2364 Equation(2)
 
-A J Hurd, D W Schaefer, J E Martin, *Phys. Rev. A*,
-35 (1987) 2361-2364 Equation(2)
+Authorship and Verification
+----------------------------
+
+* **Converted to sasmodels by:** Piotr Rozyczko **Date:** Jan 20, 2016
+* **Last Reviewed by:** Richard Heenan **Date:** May 30, 2018
 """
 
@@ -67,5 +78,4 @@
         rg_primary    =  rg
         background   =  background
-        Ref: Schmidt, J Appl Cryst, eq(19), (1991), 24, 414-435
         Hurd, Schaefer, Martin, Phys Rev A, eq(2),(1987),35, 2361-2364
         Note that 0 < Ds< 6 and 0 < Dm < 6.
@@ -78,6 +88,6 @@
     ["fractal_dim_mass", "",      1.8, [0.0, 6.0], "", "Mass fractal dimension"],
     ["fractal_dim_surf", "",      2.3, [0.0, 6.0], "", "Surface fractal dimension"],
-    ["rg_cluster",       "Ang",  86.7, [0.0, inf], "", "Cluster radius of gyration"],
-    ["rg_primary",       "Ang", 4000., [0.0, inf], "", "Primary particle radius of gyration"],
+    ["rg_cluster",       "Ang", 4000., [0.0, inf], "", "Cluster radius of gyration"],
+    ["rg_primary",       "Ang",  86.7, [0.0, inf], "", "Primary particle radius of gyration"],
 ]
 # pylint: enable=bad-whitespace, line-too-long
@@ -107,11 +117,11 @@
             fractal_dim_mass=1.8,
             fractal_dim_surf=2.3,
-            rg_cluster=86.7,
-            rg_primary=4000.0)
+            rg_cluster=4000.0,
+            rg_primary=86.7)
 
 tests = [
 
-    # Accuracy tests based on content in test/utest_other_models.py
-    [{'fractal_dim_mass':      1.8,
+    # Accuracy tests based on content in test/utest_other_models.py  All except first, changed so rg_cluster is the larger, RKH 30 May 2018
+    [{'fractal_dim_mass':   1.8,
       'fractal_dim_surf':   2.3,
       'rg_cluster':   86.7,
@@ -123,21 +133,21 @@
     [{'fractal_dim_mass':      3.3,
       'fractal_dim_surf':   1.0,
-      'rg_cluster':   90.0,
-      'rg_primary': 4000.0,
-     }, 0.001, 0.18562699016],
+      'rg_cluster': 4000.0,
+      'rg_primary':   90.0,
+     }, 0.001, 0.0932516614456],
 
     [{'fractal_dim_mass':      1.3,
-      'fractal_dim_surf':   1.0,
-      'rg_cluster':   90.0,
-      'rg_primary': 2000.0,
+      'fractal_dim_surf':   2.0,
+      'rg_cluster': 2000.0,
+      'rg_primary':   90.0,
       'background':    0.8,
-     }, 0.001, 1.16539753641],
+     }, 0.001, 1.28296431786],
 
     [{'fractal_dim_mass':      2.3,
-      'fractal_dim_surf':   1.0,
-      'rg_cluster':   90.0,
-      'rg_primary': 1000.0,
+      'fractal_dim_surf':   3.1,
+      'rg_cluster':  1000.0,
+      'rg_primary':  30.0,
       'scale':        10.0,
       'background':    0.0,
-     }, 0.051, 0.000169548800377],
+     }, 0.051, 0.00333804044899],
     ]
Index: sasmodels/models/parallelepiped.c
===================================================================
--- sasmodels/models/parallelepiped.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
+++ sasmodels/models/parallelepiped.c	(revision dbf1a60059f2daf9ef240db7a7fbf91cf5de9b9e)
@@ -38,4 +38,5 @@
             inner_total += GAUSS_W[j] * square(si1 * si2);
         }
+        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5
         inner_total *= 0.5;
 
@@ -43,4 +44,5 @@
         outer_total += GAUSS_W[i] * inner_total * si * si;
     }
+    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5
     outer_total *= 0.5;
 
Index: sasmodels/models/parallelepiped.py
===================================================================
--- sasmodels/models/parallelepiped.py	(revision ef07e950965796e2a3a8156fcade98ce5e403001)
+++ sasmodels/models/parallelepiped.py	(revision f89ec96583ec917dee91d159ad08cdc6c4222b10)
@@ -2,14 +2,12 @@
 # Note: model title and parameter table are inserted automatically
 r"""
-The form factor is normalized by the particle volume.
-For information about polarised and magnetic scattering, see
-the :ref:`magnetism` documentation.
-
 Definition
 ----------
 
- This model calculates the scattering from a rectangular parallelepiped
- (\:numref:`parallelepiped-image`\).
- If you need to apply polydispersity, see also :ref:`rectangular-prism`.
+This model calculates the scattering from a rectangular solid
+(:numref:`parallelepiped-image`).
+If you need to apply polydispersity, see also :ref:`rectangular-prism`. For
+information about polarised and magnetic scattering, see
+the :ref:`magnetism` documentation.
 
 .. _parallelepiped-image:
@@ -21,10 +19,13 @@
 
 The three dimensions of the parallelepiped (strictly here a cuboid) may be
-given in *any* size order. To avoid multiple fit solutions, especially
-with Monte-Carlo fit methods, it may be advisable to restrict their ranges.
-There may be a number of closely similar "best fits", so some trial and
-error, or fixing of some dimensions at expected values, may help.
-
-The 1D scattering intensity $I(q)$ is calculated as:
+given in *any* size order as long as the particles are randomly oriented (i.e.
+take on all possible orientations see notes on 2D below). To avoid multiple fit
+solutions, especially with Monte-Carlo fit methods, it may be advisable to
+restrict their ranges. There may be a number of closely similar "best fits", so
+some trial and error, or fixing of some dimensions at expected values, may
+help.
+
+The form factor is normalized by the particle volume and the 1D scattering
+intensity $I(q)$ is then calculated as:
 
 .. Comment by Miguel Gonzalez:
@@ -39,14 +40,16 @@
 
     I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2
-           \left< P(q, \alpha) \right> + \text{background}
+           \left< P(q, \alpha, \beta) \right> + \text{background}
 
 where the volume $V = A B C$, the contrast is defined as
-$\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$,
-$P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented
-at an angle $\alpha$ (angle between the long axis C and $\vec q$),
-and the averaging $\left<\ldots\right>$ is applied over all orientations.
+$\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, $P(q, \alpha, \beta)$
+is the form factor corresponding to a parallelepiped oriented
+at an angle $\alpha$ (angle between the long axis C and $\vec q$), and $\beta$
+(the angle between the projection of the particle in the $xy$ detector plane
+and the $y$ axis) and the averaging $\left<\ldots\right>$ is applied over all
+orientations.
 
 Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the
-form factor is given by (Mittelbach and Porod, 1961)
+form factor is given by (Mittelbach and Porod, 1961 [#Mittelbach]_)
 
 .. math::
@@ -66,75 +69,13 @@
     \mu &= qB
 
-The scattering intensity per unit volume is returned in units of |cm^-1|.
-
-NB: The 2nd virial coefficient of the parallelepiped is calculated based on
-the averaged effective radius, after appropriately sorting the three
-dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and
-length $(= C)$ values, and used as the effective radius for
-$S(q)$ when $P(q) \cdot S(q)$ is applied.
-
-For 2d data the orientation of the particle is required, described using
-angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
-of the calculation and angular dispersions see :ref:`orientation` .
-
-.. Comment by Miguel Gonzalez:
-   The following text has been commented because I think there are two
-   mistakes. Psi is the rotational angle around C (but I cannot understand
-   what it means against the q plane) and psi=0 corresponds to a||x and b||y.
-
-   The angle $\Psi$ is the rotational angle around the $C$ axis against
-   the $q$ plane. For example, $\Psi = 0$ when the $B$ axis is parallel
-   to the $x$-axis of the detector.
-
-The angle $\Psi$ is the rotational angle around the $C$ axis.
-For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the $B$ axis
-oriented parallel to the y-axis of the detector with $A$ along the x-axis.
-For other $\theta$, $\phi$ values, the parallelepiped has to be first rotated
-$\theta$ degrees in the $z-x$ plane and then $\phi$ degrees around the $z$ axis,
-before doing a final rotation of $\Psi$ degrees around the resulting $C$ axis
-of the particle to obtain the final orientation of the parallelepiped.
-
-.. _parallelepiped-orientation:
-
-.. figure:: img/parallelepiped_angle_definition.png
-
-    Definition of the angles for oriented parallelepiped, shown with $A<B<C$.
-
-.. figure:: img/parallelepiped_angle_projection.png
-
-    Examples of the angles for an oriented parallelepiped against the
-    detector plane.
-
-On introducing "Orientational Distribution" in the angles, "distribution of
-theta" and "distribution of phi" parameters will appear. These are actually
-rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped,
-perpendicular to the $a$ x $c$ and $b$ x $c$ faces. (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,
-perpendicular to the $a$ x $b$ face. Some experimentation may be required to
-understand the 2d patterns fully as discussed in :ref:`orientation` .
-
-For a given orientation of the parallelepiped, the 2D form factor is
-calculated as
-
-.. math::
-
-    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2
-                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2
-                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2
-
-with
-
-.. math::
-
-    \cos\alpha &= \hat A \cdot \hat q, \\
-    \cos\beta  &= \hat B \cdot \hat q, \\
-    \cos\gamma &= \hat C \cdot \hat q
-
-and the scattering intensity as:
-
-.. math::
-
-    I(q_x, q_y) = \frac{\text{scale}}{V} V^2 \Delta\rho^2 P(q_x, q_y)
+where substitution of $\sigma = cos\alpha$ and $\beta = \pi/2 \ u$ have been
+applied.
+
+For **oriented** particles, the 2D scattering intensity, $I(q_x, q_y)$, is
+given as:
+
+.. math::
+
+    I(q_x, q_y) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 P(q_x, q_y)
             + \text{background}
 
@@ -148,4 +89,81 @@
    with scale being the volume fraction.
 
+Where $P(q_x, q_y)$ for a given orientation of the form factor is calculated as
+
+.. math::
+
+    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}
+                   {2}qA\cos\alpha)}\right]^2
+                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}
+                   {2}qB\cos\beta)}\right]^2
+                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}
+                   {2}qC\cos\gamma)}\right]^2
+
+with
+
+.. math::
+
+    \cos\alpha &= \hat A \cdot \hat q, \\
+    \cos\beta  &= \hat B \cdot \hat q, \\
+    \cos\gamma &= \hat C \cdot \hat q
+
+
+FITTING NOTES
+~~~~~~~~~~~~~
+
+#. The 2nd virial coefficient of the parallelepiped is calculated based on
+   the averaged effective radius, after appropriately sorting the three
+   dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and
+   length $(= C)$ values, and used as the effective radius for
+   $S(q)$ when $P(q) \cdot S(q)$ is applied.
+
+#. For 2d data the orientation of the particle is required, described using
+   angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, where $\theta$
+   and $\phi$ define the orientation of the director in the laboratry reference
+   frame of the beam direction ($z$) and detector plane ($x-y$ plane), while
+   the angle $\Psi$ is effectively the rotational angle around the particle
+   $C$ axis. For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the
+   $B$ axis oriented parallel to the y-axis of the detector with $A$ along
+   the x-axis. For other $\theta$, $\phi$ values, the order of rotations
+   matters. In particular, the parallelepiped must first be rotated $\theta$
+   degrees in the $x-z$ plane before rotating $\phi$ degrees around the $z$
+   axis (in the $x-y$ plane). Applying orientational distribution to the
+   particle orientation (i.e  `jitter` to one or more of these angles) can get
+   more confusing as `jitter` is defined **NOT** with respect to the laboratory
+   frame but the particle reference frame. It is thus highly recmmended to
+   read :ref:`orientation` for further details of the calculation and angular
+   dispersions.
+
+.. note:: For 2d, constraints must be applied during fitting to ensure that the
+   order of sides chosen is not altered, and hence that the correct definition
+   of angles is preserved. For the default choice shown here, that means
+   ensuring that the inequality $A < B < C$ is not violated,  The calculation
+   will not report an error, but the results may be not correct.
+   
+.. _parallelepiped-orientation:
+
+.. figure:: img/parallelepiped_angle_definition.png
+
+    Definition of the angles for oriented parallelepiped, shown with $A<B<C$.
+
+.. figure:: img/parallelepiped_angle_projection.png
+
+    Examples of the angles for an oriented parallelepiped against the
+    detector plane.
+
+.. Comment by Paul Butler
+   I am commenting this section out as we are trying to minimize the amount of
+   oritentational detail here and encourage the user to go to the full
+   orientation documentation so that changes can be made in just one place.
+   below is the commented paragrah:
+   On introducing "Orientational Distribution" in the angles, "distribution of
+   theta" and "distribution of phi" parameters will appear. These are actually
+   rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped,
+   perpendicular to the $a$ x $c$ and $b$ x $c$ faces. (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,
+   perpendicular to the $a$ x $b$ face. Some experimentation may be required to
+   understand the 2d patterns fully as discussed in :ref:`orientation` .
+
 
 Validation
@@ -156,11 +174,10 @@
 angles.
 
-
 References
 ----------
 
-P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
-
-R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
+.. [#Mittelbach] P Mittelbach and G Porod, *Acta Physica Austriaca*,
+   14 (1961) 185-211
+.. [#] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
 
 Authorship and Verification
@@ -169,5 +186,6 @@
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
 * **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017
-* **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017
+* **Last Reviewed by:**  Miguel Gonzales and Paul Butler **Date:** May 24,
+  2018 - documentation updated
 """
 
Index: sasmodels/resolution.py
===================================================================
--- sasmodels/resolution.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/resolution.py	(revision 0b9c6df46513d1c375b80d4d720164332391302a)
@@ -20,4 +20,5 @@
 MINIMUM_RESOLUTION = 1e-8
 MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data
+PINHOLE_N_SIGMA = 2.5 # From: Barker & Pedersen 1995 JAC
 
 class Resolution(object):
@@ -65,6 +66,9 @@
     *q_calc* is the list of points to calculate, or None if this should
     be estimated from the *q* and *q_width*.
-    """
-    def __init__(self, q, q_width, q_calc=None, nsigma=3):
+
+    *nsigma* is the width of the resolution function.  Should be 2.5.
+    See :func:`pinhole_resolution` for details.
+    """
+    def __init__(self, q, q_width, q_calc=None, nsigma=PINHOLE_N_SIGMA):
         #*min_step* is the minimum point spacing to use when computing the
         #underlying model.  It should be on the order of
@@ -82,12 +86,15 @@
 
         # Protect against models which are not defined for very low q.  Limit
-        # the smallest q value evaluated (in absolute) to 0.02*min
+        # the smallest q value evaluated to 0.02*min.  Note that negative q
+        # values are trimmed even for broad resolution.  Although not possible
+        # from the geometry, they may appear since we are using a truncated
+        # gaussian to represent resolution rather than a skew distribution.
         cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q)
-        self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff]
+        self.q_calc = self.q_calc[self.q_calc >= cutoff]
 
         # Build weight matrix from calculated q values
         self.weight_matrix = pinhole_resolution(
-            self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION))
-        self.q_calc = abs(self.q_calc)
+            self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION),
+            nsigma=nsigma)
 
     def apply(self, theory):
@@ -101,7 +108,7 @@
     *q* points at which the data is measured.
 
-    *dqx* slit width in qx
-
-    *dqy* slit height in qy
+    *qx_width* slit width in qx
+
+    *qy_width* slit height in qy
 
     *q_calc* is the list of points to calculate, or None if this should
@@ -154,6 +161,6 @@
 
 
-def pinhole_resolution(q_calc, q, q_width):
-    """
+def pinhole_resolution(q_calc, q, q_width, nsigma=PINHOLE_N_SIGMA):
+    r"""
     Compute the convolution matrix *W* for pinhole resolution 1-D data.
 
@@ -162,5 +169,16 @@
     *W*, the resolution smearing can be computed using *dot(W,q)*.
 
+    Note that resolution is limited to $\pm 2.5 \sigma$.[1]  The true resolution
+    function is a broadened triangle, and does not extend over the entire
+    range $(-\infty, +\infty)$.  It is important to impose this limitation
+    since some models fall so steeply that the weighted value in gaussian
+    tails would otherwise dominate the integral.
+
     *q_calc* must be increasing.  *q_width* must be greater than zero.
+
+    [1] Barker, J. G., and J. S. Pedersen. 1995. Instrumental Smearing Effects
+    in Radially Symmetric Small-Angle Neutron Scattering by Numerical and
+    Analytical Methods. Journal of Applied Crystallography 28 (2): 105--14.
+    https://doi.org/10.1107/S0021889894010095.
     """
     # The current algorithm is a midpoint rectangle rule.  In the test case,
@@ -170,4 +188,10 @@
     cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :])
     weights = cdf[1:] - cdf[:-1]
+    # Limit q range to +/- 2.5 sigma
+    qhigh = q + nsigma*q_width
+    #qlow = q - nsigma*q_width  # linear limits
+    qlow = q*q/qhigh  # log limits
+    weights[q_calc[:, None] < qlow[None, :]] = 0.
+    weights[q_calc[:, None] > qhigh[None, :]] = 0.
     weights /= np.sum(weights, axis=0)[None, :]
     return weights
@@ -494,11 +518,14 @@
 
 
-def gaussian(q, q0, dq):
-    """
-    Return the Gaussian resolution function.
+def gaussian(q, q0, dq, nsigma=2.5):
+    """
+    Return the truncated Gaussian resolution function.
 
     *q0* is the center, *dq* is the width and *q* are the points to evaluate.
     """
-    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)
+    # Calculate the density of the tails; the resulting gaussian needs to be
+    # scaled by this amount in ordere to integrate to 1.0
+    two_tail_density = 2 * (1 + erf(-nsigma/sqrt(2)))/2
+    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)/(1-two_tail_density)
 
 
@@ -558,5 +585,5 @@
 
 
-def romberg_pinhole_1d(q, q_width, form, pars, nsigma=5):
+def romberg_pinhole_1d(q, q_width, form, pars, nsigma=2.5):
     """
     Romberg integration for pinhole resolution.
@@ -678,4 +705,7 @@
         np.testing.assert_equal(output, self.y)
 
+    # TODO: turn pinhole/slit demos into tests
+
+    @unittest.skip("suppress comparison with old version; pinhole calc changed")
     def test_pinhole(self):
         """
@@ -686,8 +716,11 @@
         theory = 12.0-1000.0*resolution.q_calc
         output = resolution.apply(theory)
+        # Note: answer came from output of previous run.  Non-integer
+        # values at ends come from the fact that q_calc does not
+        # extend beyond q, and so the weights don't balance.
         answer = [
-            10.44785079, 9.84991299, 8.98101708,
-            7.99906585, 6.99998311, 6.00001689,
-            5.00093415, 4.01898292, 3.15008701, 2.55214921,
+            10.47037734, 9.86925860,
+            9., 8., 7., 6., 5., 4.,
+            3.13074140, 2.52962266,
             ]
         np.testing.assert_allclose(output, answer, atol=1e-8)
@@ -732,4 +765,5 @@
         self._compare(q, output, answer, 1e-6)
 
+    @unittest.skip("suppress comparison with old version; pinhole calc changed")
     def test_pinhole(self):
         """
@@ -746,4 +780,5 @@
         self._compare(q, output, answer, 3e-4)
 
+    @unittest.skip("suppress comparison with old version; pinhole calc changed")
     def test_pinhole_romberg(self):
         """
@@ -761,5 +796,5 @@
         #                     2*np.pi/pars['radius']/200)
         #tol = 0.001
-        ## The default 3 sigma and no extra points gets 1%
+        ## The default 2.5 sigma and no extra points gets 1%
         q_calc = None  # type: np.ndarray
         tol = 0.01
@@ -1080,5 +1115,5 @@
 
     if isinstance(resolution, Slit1D):
-        width, height = resolution.dqx, resolution.dqy
+        width, height = resolution.qx_width, resolution.qy_width
         Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars)
     else:
Index: sasmodels/rst2html.py
===================================================================
--- sasmodels/rst2html.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/rst2html.py	(revision 1fbadb20289a15dcd2b077eb8ef47e0dbda5eef1)
@@ -56,5 +56,7 @@
     # others don't work properly with math_output!
     if math_output == "mathjax":
-        settings = {"math_output": math_output}
+        # TODO: this is copied from docs/conf.py; there should be only one
+        mathjax_path = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-MML-AM_CHTML"
+        settings = {"math_output": math_output + " " + mathjax_path}
     else:
         settings = {"math-output": math_output}
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision 05df1ded9ed7edca9420c90b0d39d936b8e1c97e)
+++ sasmodels/sasview_model.py	(revision d533590d9857d7aef962616cfdaa6ecfc2e6751d)
@@ -686,5 +686,5 @@
         self._intermediate_results = getattr(calculator, 'results', None)
         calculator.release()
-        self._model.release()
+        #self._model.release()
         return result
 
@@ -719,5 +719,5 @@
 
     def set_dispersion(self, parameter, dispersion):
-        # type: (str, weights.Dispersion) -> Dict[str, Any]
+        # type: (str, weights.Dispersion) -> None
         """
         Set the dispersion object for a model parameter
@@ -742,5 +742,6 @@
             self.dispersion[parameter] = dispersion.get_pars()
         else:
-            raise ValueError("%r is not a dispersity or orientation parameter")
+            raise ValueError("%r is not a dispersity or orientation parameter"
+                             % parameter)
 
     def _dispersion_mesh(self):
Index: tup_py2exe.py
===================================================================
--- setup_py2exe.py	(revision 87985ca5989a70a24a7e9ac09a971d544fe16dfe)
+++ 	(revision )
@@ -1,322 +1,0 @@
-#!/usr/bin/env python
-
-# Copyright (C) 2006-2010, University of Maryland
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/ or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in
-# all copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-# THE SOFTWARE.
-
-# Author: James Krycka
-
-"""
-This script uses py2exe to create dist\bumps.exe and dist\bumps_gui.exe for
-running the Bumps application in either CLI or GUI mode.
-
-These executables start the application and import the rest of the application
-code stored in library.zip.  The python interpreter and other required python
-packages and dlls are also placed in the zip file.  Additional resource files
-that are needed when Bumps is run are copied to the dist directory tree.  On
-completion, the contents of the dist directory tree can be used by the Inno
-Setup Compiler (via a separate script) to build a Windows installer/uninstaller
-for deployment of the Bumps application.  For testing purposes, bumps.exe or
-bumps_gui.exe can be run from the dist directory.
-"""
-
-import os
-import sys
-
-#sys.dont_write_bytecode = True
-
-# Force build before continuing
-os.system('"%s" setup.py build' % sys.executable)
-
-# Remove the current directory from the python path
-here = os.path.abspath(os.path.dirname(__file__))
-sys.path = [p for p in sys.path if os.path.abspath(p) != here]
-
-import glob
-
-from distutils.core import setup
-from distutils.util import get_platform
-
-# Augment the setup interface with the py2exe command and make sure the py2exe
-# option is passed to setup.
-import py2exe
-
-if len(sys.argv) == 1:
-    sys.argv.append('py2exe')
-
-# Put the build lib on the start of the path.
-# For packages with binary extensions, need platform.  If it is a pure
-# script library, use an empty platform string.
-platform = '.%s-%s' % (get_platform(), sys.version[:3])
-#platform = ''
-build_lib = os.path.abspath('build/lib' + platform)
-sys.path.insert(0, build_lib)
-
-# print "\n".join(sys.path)
-
-import wx
-import matplotlib
-matplotlib.use('WXAgg')
-import periodictable
-
-# Retrieve the application version string.
-import bumps
-version = bumps.__version__
-from bumps.gui.resources import resources as gui_resources
-
-# A manifest is required to be included in a py2exe image (or accessible as a
-# file in the image directory) when wxPython is included so that the Windows XP
-# theme is used when rendering wx widgets.  The manifest must be matched to the
-# version of Python that is being used.
-#
-# Create a manifest for use with Python 2.5 on Windows XP or Vista.  It is
-# adapted from the Python manifest file (C:\Python25\pythonw.exe.manifest).
-
-manifest_for_python25 = """
-<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
-<assembly xmlns="urn:schemas-microsoft-com:asm.v1" manifestVersion="1.0">
-<assemblyIdentity
-    version="1.0.0.0"
-    processorArchitecture="x86"
-    name="%(prog)s"
-    type="win32"
-/>
-<description>%(prog)s</description>
-<dependency>
-    <dependentAssembly>
-        <assemblyIdentity
-            type="win32"
-            name="Microsoft.Windows.Common-Controls"
-            version="6.0.0.0"
-            processorArchitecture="X86"
-            publicKeyToken="6595b64144ccf1df"
-            language="*"
-        />
-    </dependentAssembly>
-</dependency>
-</assembly>
-"""
-
-# Create a manifest for use with Python 2.6 or 2.7 on Windows XP or Vista.
-
-manifest_for_python26 = """
-<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
-<assembly xmlns="urn:schemas-microsoft-com:asm.v1" manifestVersion="1.0">
-  <assemblyIdentity
-    version="5.0.0.0"
-    processorArchitecture="x86"
-    name="%(prog)s"
-    type="win32">
-  </assemblyIdentity>
-  <description>%(prog)s</description>
-  <trustInfo xmlns="urn:schemas-microsoft-com:asm.v3">
-    <security>
-      <requestedPrivileges>
-        <requestedExecutionLevel
-          level="asInvoker"
-          uiAccess="false">
-        </requestedExecutionLevel>
-      </requestedPrivileges>
-    </security>
-  </trustInfo>
-  <dependency>
-    <dependentAssembly>
-      <assemblyIdentity
-        type="win32"
-        name="Microsoft.VC90.CRT"
-        version="9.0.21022.8"
-        processorArchitecture="x86"
-        publicKeyToken="1fc8b3b9a1e18e3b">
-      </assemblyIdentity>
-    </dependentAssembly>
-  </dependency>
-  <dependency>
-    <dependentAssembly>
-      <assemblyIdentity
-        type="win32"
-        name="Microsoft.Windows.Common-Controls"
-        version="6.0.0.0"
-        processorArchitecture="x86"
-        publicKeyToken="6595b64144ccf1df"
-        language="*">
-      </assemblyIdentity>
-    </dependentAssembly>
-  </dependency>
-</assembly>
-"""
-
-# Select the appropriate manifest to use.
-if sys.version_info >= (3, 0) or sys.version_info < (2, 5):
-    print("*** This script only works with Python 2.5, 2.6, or 2.7.")
-    sys.exit()
-elif sys.version_info >= (2, 6):
-    manifest = manifest_for_python26
-elif sys.version_info >= (2, 5):
-    manifest = manifest_for_python25
-
-# Create a list of all files to include along side the executable being built
-# in the dist directory tree.  Each element of the data_files list is a tuple
-# consisting of a path (relative to dist\) and a list of files in that path.
-data_files = []
-
-# Add resource files that need to reside in the same directory as the image.
-data_files.append(('.', [os.path.join('.', 'LICENSE.txt')]))
-data_files.append(('.', [os.path.join('.', 'README.txt')]))
-data_files.append(('.', [os.path.join('.', 'bin', 'bumps_launch.bat')]))
-
-# Add application specific data files from the bumps\bumps-data folder.
-data_files += gui_resources.data_files()
-
-# Add data files from the matplotlib\mpl-data folder and its subfolders.
-# For matploblib prior to version 0.99 see the examples at the end of the file.
-data_files += matplotlib.get_py2exe_datafiles()
-
-# Add data files from the periodictable\xsf folder.
-data_files += periodictable.data_files()
-
-# Add example directories and their files.  An empty directory is ignored.
-# Note that Inno Setup will determine where these files will be placed such as
-# C:\My Documents\... instead of the installation folder.
-for path in glob.glob(os.path.join('examples', '*')):
-    if os.path.isdir(path):
-        for file in glob.glob(os.path.join(path, '*.*')):
-            data_files.append((path, [file]))
-    else:
-        data_files.append(('examples', [path]))
-
-for path in glob.glob(os.path.join('doc', 'examples', '*')):
-    if os.path.isdir(path):
-        for file in glob.glob(os.path.join(path, '*.*')):
-            data_files.append((path, [file]))
-    else:
-        data_files.append(('doc', [path]))
-
-# Add PDF documentation to the dist staging directory.
-pdf = os.path.join('doc', 'Bumps.pdf')
-if os.path.isfile(pdf):
-    data_files.append(('doc', [pdf]))
-else:
-    print("*** %s not found - building frozen image without it ***" % pdf)
-
-# Add the Microsoft Visual C++ 2008 redistributable kit if we are building with
-# Python 2.6 or 2.7.  This kit will be installed on the target system as part
-# of the installation process for the frozen image.  Note that the Python 2.5
-# interpreter requires msvcr71.dll which is included in the Python25 package,
-# however, Python 2.6 and 2.7 require the msvcr90.dll but they do not bundle it
-# with the Python26 or Python27 package.  Thus, for Python 2.6 and later, the
-# appropriate dll must be present on the target system at runtime.
-if sys.version_info >= (2, 6):
-    pypath = os.path.dirname(sys.executable)
-    data_files.append(('.', [os.path.join(pypath, 'vcredist_x86.exe')]))
-
-# Specify required packages to bundle in the executable image.
-packages = ['numpy', 'scipy', 'matplotlib', 'pytz', 'pyparsing',
-            'periodictable', 'bumps', 'sasmodels', 'pyopencl',
-            ]
-
-# Specify files to include in the executable image.
-includes = []
-
-# Specify files to exclude from the executable image.
-# - We can safely exclude Tk/Tcl and Qt modules because our app uses wxPython.
-# - We do not use ssl services so they are omitted.
-# - We can safely exclude the TkAgg matplotlib backend because our app uses
-#   "matplotlib.use('WXAgg')" to override the default matplotlib configuration.
-# - On the web it is widely recommended to exclude certain lib*.dll modules
-#   but this does not seem necessary any more (but adding them does not hurt).
-# - Python25 requires mscvr71.dll, however, Win XP includes this file.
-# - Since we do not support Win 9x systems, w9xpopen.dll is not needed.
-# - For some reason cygwin1.dll gets included by default, but it is not needed.
-
-excludes = ['Tkinter', 'PyQt4', '_ssl', '_tkagg', 'numpy.distutils.test']
-
-dll_excludes = ['libgdk_pixbuf-2.0-0.dll',
-                'libgobject-2.0-0.dll',
-                'libgdk-win32-2.0-0.dll',
-                'tcl84.dll',
-                'tk84.dll',
-                'QtGui4.dll',
-                'QtCore4.dll',
-                'msvcr71.dll',
-                'msvcp90.dll',
-                'w9xpopen.exe',
-                'cygwin1.dll']
-
-
-class Target(object):
-
-    """This class stores metadata about the distribution in a dictionary."""
-
-    def __init__(self, **kw):
-        self.__dict__.update(kw)
-        self.version = version
-
-clientCLI = Target(
-    name='Bumps',
-    description='Bumps CLI application',
-    # module to run on application start
-    script=os.path.join('bin', 'bumps_cli.py'),
-    dest_base='bumps',  # file name part of the exe file to create
-    # also need to specify in data_files
-    icon_resources=[
-        (1, os.path.join('bumps', 'gui', 'resources', 'bumps.ico'))],
-    bitmap_resources=[],
-    other_resources=[(24, 1, manifest % dict(prog='Bumps'))])
-
-clientGUI = Target(
-    name='Bumps',
-    description='Bumps GUI application',
-    # module to run on application start
-    script=os.path.join('bin', 'bumps_gui.py'),
-    dest_base='bumps_gui',  # file name part of the exe file to create
-    # also need to specify in data_files
-    icon_resources=[
-        (1, os.path.join('bumps', 'gui', 'resources', 'bumps.ico'))],
-    bitmap_resources=[],
-    other_resources=[(24, 1, manifest % dict(prog='Bumps'))])
-
-# Now we do the work to create a standalone distribution using py2exe.
-#
-# When the application is run in console mode, a console window will be created
-# to receive any logging or error messages and the application will then create
-# a separate GUI application window.
-#
-# When the application is run in windows mode, it will create a GUI application
-# window and no console window will be provided.  Output to stderr will be
-# written to <app-image-name>.log.
-setup(
-    console=[clientCLI],
-    windows=[clientGUI],
-    options={'py2exe': {
-        'packages': packages,
-        'includes': includes,
-        'excludes': excludes,
-        'dll_excludes': dll_excludes,
-        'compressed': 1,   # standard compression
-        'optimize': 0,     # no byte-code optimization
-        'dist_dir': "dist",  # where to put py2exe results
-        'xref': False,     # display cross reference (as html doc)
-        'bundle_files': 1,  # bundle python25.dll in library.zip
-    }
-    },
-    # Since we are building two exe's, do not put the shared library in each
-    # of them.  Instead create a single, separate library.zip file.
-    # zipfile=None,               # bundle library.zip in exe
-    data_files=data_files           # list of files to copy to dist directory
-)
