diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 6e7a07ef..ff0441e5 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -59,3 +59,12 @@ jobs:
       if: ${{ matrix.os == 'ubuntu-latest' }}
       run: |
         make -j 4 -C doc SPHINXOPTS="-W --keep-going -n" html
+
+    - name: Publish samodels docs
+      if: ${{ matrix.os == 'ubuntu-latest' && matrix.python-version == '3.10'}}
+      uses: actions/upload-artifact@v3
+      with:
+        name: sasmodels-docs-${{ matrix.os }}-${{ matrix.python-version }}
+        path: |
+          doc/_build/html
+        if-no-files-found: error
diff --git a/doc/guide/theory.rst b/doc/guide/theory.rst
index 0bd3f371..50786723 100644
--- a/doc/guide/theory.rst
+++ b/doc/guide/theory.rst
@@ -1,3 +1,4 @@
+.. currentmodule:: sasmodels
 .. theory.rst
 
 .. Much of the following text was scraped from fitting_sq.py
diff --git a/sasmodels/direct_model.py b/sasmodels/direct_model.py
index e172de6b..0bff6600 100644
--- a/sasmodels/direct_model.py
+++ b/sasmodels/direct_model.py
@@ -121,8 +121,10 @@ def get_mesh(model_info, values, dim='1d', mono=False):
     values = values.copy()
     mesh = [_pop_par_weights(p, values, active(p.name))
             for p in parameters.call_parameters]
+
     if values:
         raise TypeError(f"Unused parameters in call: {', '.join(values.keys())}")
+
     return mesh
 
 
@@ -561,16 +563,17 @@ def near(value, target):
         return np.allclose(value, target, rtol=1e-6, atol=0, equal_nan=True)
     # Note: target values taken from running main() on parameters.
     # Resolution was 5% dq/q.
-    pars = dict(radius=200)
+    pars = dict(radius=200, background=0)  # default background=1e-3, scale=1
     # simple sphere in 1D (perfect, pinhole, slit)
-    assert near(Iq('sphere', [0.1], **pars), [0.6200146273894904])
-    assert near(Iq('sphere', [0.1], dq=[0.005], **pars), [2.3019224683980215])
-    assert near(Iq('sphere', [0.1], qw=[0.005], ql=[1.0], **pars), [0.3673431784535172])
+    perfect_target = 0.6190146273894904
+    assert near(Iq('sphere', [0.1], **pars), [perfect_target])
+    assert near(Iq('sphere', [0.1], dq=[0.005], **pars), [2.3009224683980215])
+    assert near(Iq('sphere', [0.1], qw=[0.005], ql=[1.0], **pars), [0.3663431784535172])
     # simple sphere in 2D (perfect, pinhole)
-    assert near(Iqxy('sphere', [0.1], [0.1], **pars), [1.1781532874802199])
+    assert near(Iqxy('sphere', [0.1], [0.1], **pars), [1.1771532874802199])
     assert near(Iqxy('sphere', [0.1], [0.1], dqx=[0.005], dqy=[0.005], **pars), 
-        [0.8177780778578667])
-    # sesans
+        [0.8167780778578667])
+    # sesans (no background or scale)
     assert near(Gxi('sphere', [100], **pars), [-0.19146959126623486])
     # Check that single point sesans matches value in an array
     xi = np.logspace(1, 3, 100)
@@ -587,6 +590,11 @@ def near(value, target):
         radius=200, radius_pd=0.1, radius_pd_n=15, radius_pd_nsigma=2.5,
         radius_pd_type="uniform")
     assert near(Iq('sphere', [0.1], **pars), [2.703169824954617])
+    # background and scale
+    background, scale = 1e-4, 0.1
+    pars = dict(radius=200, background=background, scale=scale)
+    assert near(Iq('sphere', [0.1], **pars), [perfect_target*scale + background])
+
 
 if __name__ == "__main__":
     import logging
diff --git a/sasmodels/jitter.py b/sasmodels/jitter.py
index cdaf64ba..b8429c0d 100755
--- a/sasmodels/jitter.py
+++ b/sasmodels/jitter.py
@@ -881,7 +881,7 @@ def build_model(model_name, n=150, qmax=0.5, **pars):
 
     # stuff the values for non-orientation parameters into the calculator
     calculator.pars = pars.copy()
-    calculator.pars.setdefault('backgound', 1e-3)
+    calculator.pars.setdefault('background', 1e-3)
 
     # fix the data limits so that we can see if the pattern fades
     # under rotation or angular dispersion
@@ -908,47 +908,53 @@ def select_calculator(model_name, n=150, size=(10, 40, 100)):
     a, b, c = size
     d_factor = 0.06  # for paracrystal models
     if model_name == 'sphere':
-        calculator = build_model('sphere', n=n, radius=c)
+        calculator = build_model(
+            'sphere', n=n, radius=c)
         a = b = c
     elif model_name == 'sc_paracrystal':
         a = b = c
         dnn = c
         radius = 0.5*c
-        calculator = build_model('sc_paracrystal', n=n, dnn=dnn,
-                                 d_factor=d_factor, radius=(1-d_factor)*radius,
-                                 background=0)
+        calculator = build_model(
+            'sc_paracrystal', n=n,
+            dnn=dnn, d_factor=d_factor, radius=(1-d_factor)*radius,
+            background=0)
     elif model_name == 'fcc_paracrystal':
         a = b = c
         # nearest neigbour distance dnn should be 2 radius, but I think the
         # model uses lattice spacing rather than dnn in its calculations
         dnn = 0.5*c
         radius = sqrt(2)/4 * c
-        calculator = build_model('fcc_paracrystal', n=n, dnn=dnn,
-                                 d_factor=d_factor, radius=(1-d_factor)*radius,
-                                 background=0)
+        calculator = build_model(
+            'fcc_paracrystal', n=n,
+            dnn=dnn, d_factor=d_factor, radius=(1-d_factor)*radius,
+            background=0)
     elif model_name == 'bcc_paracrystal':
         a = b = c
         # nearest neigbour distance dnn should be 2 radius, but I think the
         # model uses lattice spacing rather than dnn in its calculations
         dnn = 0.5*c
         radius = sqrt(3)/2 * c
-        calculator = build_model('bcc_paracrystal', n=n, dnn=dnn,
-                                 d_factor=d_factor, radius=(1-d_factor)*radius,
-                                 background=0)
+        calculator = build_model(
+            'bcc_paracrystal', n=n,
+            dnn=dnn, d_factor=d_factor, radius=(1-d_factor)*radius,
+            background=0)
     elif model_name == 'cylinder':
-        calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c)
+        calculator = build_model(
+            'cylinder', n=n, qmax=0.3, radius=b, length=c)
         a = b
     elif model_name == 'ellipsoid':
-        calculator = build_model('ellipsoid', n=n, qmax=1.0,
-                                 radius_polar=c, radius_equatorial=b)
+        calculator = build_model(
+            'ellipsoid', n=n, qmax=1.0,
+            radius_polar=c, radius_equatorial=b)
         a = b
     elif model_name == 'triaxial_ellipsoid':
-        calculator = build_model('triaxial_ellipsoid', n=n, qmax=0.5,
-                                 radius_equat_minor=a,
-                                 radius_equat_major=b,
-                                 radius_polar=c)
+        calculator = build_model(
+            'triaxial_ellipsoid', n=n, qmax=0.5,
+            radius_equat_minor=a, radius_equat_major=b, radius_polar=c)
     elif model_name == 'parallelepiped':
-        calculator = build_model('parallelepiped', n=n, a=a, b=b, c=c)
+        calculator = build_model(
+            'parallelepiped', n=n, length_a=a, length_b=b, length_c=c)
     else:
         raise ValueError("unknown model %s"%model_name)
 
diff --git a/sasmodels/models/lamellar_hg.py b/sasmodels/models/lamellar_hg.py
index ebdd183b..88767313 100644
--- a/sasmodels/models/lamellar_hg.py
+++ b/sasmodels/models/lamellar_hg.py
@@ -29,7 +29,7 @@
 and $\Delta\rho_T$ is tail contrast (*sld* $-$ *sld_solvent*).
 
 The total thickness of the lamellar sheet is
-a_H + \delta_T + \delta_T + \delta_H$. Note that in a non aqueous solvent
+$a_H + \delta_T + \delta_T + \delta_H$. Note that in a non aqueous solvent
 the chemical "head" group may be the "Tail region" and vice-versa.
 
 The 2D scattering intensity is calculated in the same way as 1D, where