Index: .travis.yml
===================================================================
--- .travis.yml	(revision 947a61e810a04157de1476d697ab71e62f108a72)
+++ .travis.yml	(revision ce8c38870db377607bac6fe577399e8ec56f634c)
@@ -24,6 +24,6 @@
     packages: opencl-headers
 before_install:
-- openssl aes-256-cbc -K $encrypted_fe6026add10a_key -iv $encrypted_fe6026add10a_iv
-  -in .travis/travis_rsa.enc -out .travis/travis_rsa -d
+- if [[ $encrypted_cb04388797b6_iv ]]; then openssl aes-256-cbc -K $encrypted_cb04388797b6_key -iv $encrypted_cb04388797b6_iv
+  -in .travis/travis_rsa.enc -out .travis/travis_rsa -d; fi;
 - echo $TRAVIS_OS_NAME
 - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
@@ -46,9 +46,9 @@
 - echo -e "Host danse.chem.utk.edu\n\tStrictHostKeyChecking no\n" >> ~/.ssh/config
 deploy:
-    skip_cleanup: true
-    provider: script
-    script: /bin/sh -ex ./deploy.sh
-    on:
-        branch: master
+  skip_cleanup: true
+  provider: script
+  script: "/bin/sh -ex ./deploy.sh"
+  on:
+    branch: master
 notifications:
   slack:
Index: deploy.sh
===================================================================
--- deploy.sh	(revision d9d77b031eaf23342c57552f4c95488b3d9ffb24)
+++ deploy.sh	(revision 9d1e3e4d7de512edd3a384aa735b57475b17a551)
@@ -1,5 +1,8 @@
-eval "$(ssh-agent -s)"
-chmod 600 .travis/travis_rsa
-ssh-add .travis/travis_rsa
-git remote add deploy git@danse.chem.utk.edu:/home/git/sasmodels
-git push deploy master
+if [[ $encrypted_cb04388797b6_iv ]]
+then
+    eval "$(ssh-agent -s)"
+    chmod 600 .travis/travis_rsa
+    ssh-add .travis/travis_rsa
+    git remote add deploy git@danse.chem.utk.edu:/home/git/sasmodels
+    git push deploy master
+fi
Index: doc/conf.py
===================================================================
--- doc/conf.py	(revision 8ae853293022b5081a1e1a8c26f74eaa2b3f8a13)
+++ doc/conf.py	(revision 30b60d21e384fd9fda8e281e787ff2869925e10d)
@@ -211,4 +211,6 @@
 #latex_preamble = ''
 LATEX_PREAMBLE=r"""
+\newcommand{\lt}{<}
+\newcommand{\gt}{>}
 \renewcommand{\AA}{\text{\r{A}}} % Allow \AA in math mode
 \usepackage[utf8]{inputenc}      % Allow unicode symbols in text
Index: doc/guide/magnetism/magnetism.rst
===================================================================
--- doc/guide/magnetism/magnetism.rst	(revision 990d8dfe5ac399564bc0bb3b873c5a1b56d4e573)
+++ doc/guide/magnetism/magnetism.rst	(revision 990d8dfe5ac399564bc0bb3b873c5a1b56d4e573)
@@ -16,5 +16,5 @@
 
 .. figure::
-    mag_img/mag_vector.bmp
+    mag_img/mag_vector.png
 
 The magnetic scattering length density is then
@@ -36,5 +36,5 @@
 
 .. figure::
-    mag_img/M_angles_pic.bmp
+    mag_img/M_angles_pic.png
 
 If the angles of the $Q$ vector and the spin-axis $x'$ to the $x$ - axis are
Index: doc/guide/plugin.rst
===================================================================
--- doc/guide/plugin.rst	(revision 870a2f45e993d96ff75196552953ec259f2ba197)
+++ doc/guide/plugin.rst	(revision 30b60d21e384fd9fda8e281e787ff2869925e10d)
@@ -117,9 +117,4 @@
 Models that do not conform to these requirements will *never* be incorporated
 into the built-in library.
-
-More complete documentation for the sasmodels package can be found at
-`<http://www.sasview.org/sasmodels>`_. In particular,
-`<http://www.sasview.org/sasmodels/api/generate.html#module-sasmodels.generate>`_
-describes the structure of a model.
 
 
@@ -613,5 +608,5 @@
 
     sas_gamma(x):
-        Gamma function $\text{sas_gamma}(x) = \Gamma(x)$.
+        Gamma function sas_gamma\ $(x) = \Gamma(x)$.
 
         The standard math function, tgamma(x) is unstable for $x < 1$
@@ -623,7 +618,7 @@
     sas_erf(x), sas_erfc(x):
         Error function
-        $\text{sas_erf}(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt$
+        sas_erf\ $(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt$
         and complementary error function
-        $\text{sas_erfc}(x) = \frac{2}{\sqrt\pi}\int_x^{\infty} e^{-t^2}\,dt$.
+        sas_erfc\ $(x) = \frac{2}{\sqrt\pi}\int_x^{\infty} e^{-t^2}\,dt$.
 
         The standard math functions erf(x) and erfc(x) are slower and broken
@@ -634,5 +629,5 @@
 
     sas_J0(x):
-        Bessel function of the first kind $\text{sas_J0}(x)=J_0(x)$ where
+        Bessel function of the first kind sas_J0\ $(x)=J_0(x)$ where
         $J_0(x) = \frac{1}{\pi}\int_0^\pi \cos(x\sin(\tau))\,d\tau$.
 
@@ -643,5 +638,5 @@
 
     sas_J1(x):
-        Bessel function of the first kind  $\text{sas_J1}(x)=J_1(x)$ where
+        Bessel function of the first kind  sas_J1\ $(x)=J_1(x)$ where
         $J_1(x) = \frac{1}{\pi}\int_0^\pi \cos(\tau - x\sin(\tau))\,d\tau$.
 
@@ -652,8 +647,8 @@
 
     sas_JN(n, x):
-        Bessel function of the first kind and integer order $n$:
-        $\text{sas_JN}(n, x)=J_n(x)$ where
+        Bessel function of the first kind and integer order $n$,
+        sas_JN\ $(n, x) =J_n(x)$ where
         $J_n(x) = \frac{1}{\pi}\int_0^\pi \cos(n\tau - x\sin(\tau))\,d\tau$.
-        If $n$ = 0 or 1, it uses sas_J0(x) or sas_J1(x), respectively.
+        If $n$ = 0 or 1, it uses sas_J0($x$) or sas_J1($x$), respectively.
 
         The standard math function jn(n, x) is not available on all platforms.
@@ -663,5 +658,5 @@
 
     sas_Si(x):
-        Sine integral $\text{Si}(x) = \int_0^x \tfrac{\sin t}{t}\,dt$.
+        Sine integral Si\ $(x) = \int_0^x \tfrac{\sin t}{t}\,dt$.
 
         This function uses Taylor series for small and large arguments:
@@ -688,5 +683,5 @@
     sas_3j1x_x(x):
         Spherical Bessel form
-        $\text{sph_j1c}(x) = 3 j_1(x)/x = 3 (\sin(x) - x \cos(x))/x^3$,
+        sph_j1c\ $(x) = 3 j_1(x)/x = 3 (\sin(x) - x \cos(x))/x^3$,
         with a limiting value of 1 at $x=0$, where $j_1(x)$ is the spherical
         Bessel function of the first kind and first order.
@@ -699,5 +694,5 @@
 
     sas_2J1x_x(x):
-        Bessel form $\text{sas_J1c}(x) = 2 J_1(x)/x$, with a limiting value
+        Bessel form sas_J1c\ $(x) = 2 J_1(x)/x$, with a limiting value
         of 1 at $x=0$, where $J_1(x)$ is the Bessel function of first kind
         and first order.
Index: doc/guide/resolution.rst
===================================================================
--- doc/guide/resolution.rst	(revision f8a2baa4aebd863d4c966fc129f4f4354fc63d5b)
+++ doc/guide/resolution.rst	(revision 30b60d21e384fd9fda8e281e787ff2869925e10d)
@@ -212,5 +212,5 @@
 elliptical Gaussian distribution. The $A$ is a normalization factor.
 
-.. figure:: resolution_2d_rotation.gif
+.. figure:: resolution_2d_rotation.png
 
     Coordinate axis rotation for 2D resolution calculation.
Index: doc/rst_prolog
===================================================================
--- doc/rst_prolog	(revision a0fb06a2ed023b879d041098faa161e70cffcd07)
+++ doc/rst_prolog	(revision 30b60d21e384fd9fda8e281e787ff2869925e10d)
@@ -1,45 +1,4 @@
 .. Set up some substitutions to make life easier...
-.. Remove |biggamma|, etc. when they are no longer needed.
 
-
-.. |alpha| unicode:: U+03B1
-.. |beta| unicode:: U+03B2
-.. |gamma| unicode:: U+03B3
-.. |delta| unicode:: U+03B4
-.. |epsilon| unicode:: U+03B5
-.. |zeta| unicode:: U+03B6
-.. |eta| unicode:: U+03B7
-.. |theta| unicode:: U+03B8
-.. |iota| unicode:: U+03B9
-.. |kappa| unicode:: U+03BA
-.. |lambda| unicode:: U+03BB
-.. |mu| unicode:: U+03BC
-.. |nu| unicode:: U+03BD
-.. |xi| unicode:: U+03BE
-.. |omicron| unicode:: U+03BF
-.. |pi| unicode:: U+03C0
-.. |rho| unicode:: U+03C1
-.. |sigma| unicode:: U+03C3
-.. |tau| unicode:: U+03C4
-.. |upsilon| unicode:: U+03C5
-.. |phi| unicode:: U+03C6
-.. |chi| unicode:: U+03C7
-.. |psi| unicode:: U+03C8
-.. |omega| unicode:: U+03C9
-
-
-.. |biggamma| unicode:: U+0393
-.. |bigdelta| unicode:: U+0394
-.. |bigzeta| unicode:: U+039E
-.. |bigpsi| unicode:: U+03A8
-.. |bigphi| unicode:: U+03A6
-.. |bigsigma| unicode:: U+03A3
-.. |Gamma| unicode:: U+0393
-.. |Delta| unicode:: U+0394
-.. |Zeta| unicode:: U+039E
-.. |Psi| unicode:: U+03A8
-
-
-.. |drho| replace:: |Delta|\ |rho|
 .. |Ang| unicode:: U+212B
 .. |Ang^-1| replace:: |Ang|\ :sup:`-1`
@@ -57,12 +16,4 @@
 .. |cm^-3| replace:: cm\ :sup:`-3`
 .. |sr^-1| replace:: sr\ :sup:`-1`
-.. |P0| replace:: P\ :sub:`0`\
-.. |A2| replace:: A\ :sub:`2`\
-
-
-.. |equiv| unicode:: U+2261
-.. |noteql| unicode:: U+2260
-.. |TM| unicode:: U+2122
-
 
 .. |cdot| unicode:: U+00B7
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision d9ec8f9cf8c3196125717d175f8a281d82eed86f)
+++ sasmodels/compare.py	(revision 765eb0ed12773d554fa60f4cd5952b5aa5da7917)
@@ -88,4 +88,5 @@
     -magnetic/-nonmagnetic* suppress magnetism
     -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
+    -neval=1 sets the number of evals for more accurate timing
 
     === precision options ===
@@ -1013,4 +1014,5 @@
     'magnetic', 'nonmagnetic',
     'accuracy=',
+    'neval=',  # for timing...
 
     # Precision options
@@ -1018,5 +1020,4 @@
     'half', 'fast', 'single', 'double', 'single!', 'double!', 'quad!',
     'sasview',  # TODO: remove sasview 3.x support
-    'timing=',
 
     # Output options
@@ -1068,5 +1069,5 @@
 
     # Plug in values given in demo
-    if use_demo:
+    if use_demo and model_info.demo:
         pars.update(model_info.demo)
     return pars
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision bb39b4a56cfca3e54321e008ce4dfecbde90f716)
+++ sasmodels/core.py	(revision 60335cc2ad4a8f97f40247b0bd0df4b5ec65f0fe)
@@ -10,4 +10,5 @@
 
 import os
+import re
 from os.path import basename, dirname, join as joinpath
 from glob import glob
@@ -21,4 +22,5 @@
 from . import kernelpy
 from . import kerneldll
+from . import custom
 
 if os.environ.get("SAS_OPENCL", "").lower() == "none":
@@ -30,4 +32,11 @@
     except Exception:
         HAVE_OPENCL = False
+
+CUSTOM_MODEL_PATH = os.environ.get('SAS_MODELPATH', "")
+if CUSTOM_MODEL_PATH == "":
+    path = joinpath(os.path.expanduser("~"), ".sasmodels", "custom_models")
+    if not os.path.isdir(path):
+        os.makedirs(path)
+    CUSTOM_MODEL_PATH = path
 
 try:
@@ -125,30 +134,60 @@
                        dtype=dtype, platform=platform)
 
-
-def load_model_info(model_name):
+def load_model_info(model_string):
     # type: (str) -> modelinfo.ModelInfo
     """
     Load a model definition given the model name.
 
-    *model_name* is the name of the model, or perhaps a model expression
-    such as sphere*hardsphere or sphere+cylinder.
+    *model_string* is the name of the model, or perhaps a model expression
+    such as sphere*cylinder or sphere+cylinder. Use '@' for a structure
+    factor product, e.g. sphere@hardsphere. Custom models can be specified by
+    prefixing the model name with 'custom.', e.g. 'custom.MyModel+sphere'.
 
     This returns a handle to the module defining the model.  This can be
     used with functions in generate to build the docs or extract model info.
     """
-    parts = model_name.split('+')
-    if len(parts) > 1:
-        model_info_list = [load_model_info(p) for p in parts]
-        return mixture.make_mixture_info(model_info_list)
-
-    parts = model_name.split('*')
-    if len(parts) > 1:
-        if len(parts) > 2:
-            raise ValueError("use P*S to apply structure factor S to model P")
-        P_info, Q_info = [load_model_info(p) for p in parts]
+    if '@' in model_string:
+        parts = model_string.split('@')
+        if len(parts) != 2:
+            raise ValueError("Use P@S to apply a structure factor S to model P")
+        P_info, Q_info = [load_model_info(part) for part in parts]
         return product.make_product_info(P_info, Q_info)
 
-    kernel_module = generate.load_kernel_module(model_name)
-    return modelinfo.make_model_info(kernel_module)
+    product_parts = []
+    addition_parts = []
+
+    addition_parts_names = model_string.split('+')
+    if len(addition_parts_names) >= 2:
+        addition_parts = [load_model_info(part) for part in addition_parts_names]
+    elif len(addition_parts_names) == 1:
+        product_parts_names = model_string.split('*')
+        if len(product_parts_names) >= 2:
+            product_parts = [load_model_info(part) for part in product_parts_names]
+        elif len(product_parts_names) == 1:
+            if "custom." in product_parts_names[0]:
+                # Extract ModelName from "custom.ModelName"
+                pattern = "custom.([A-Za-z0-9_-]+)"
+                result = re.match(pattern, product_parts_names[0])
+                if result is None:
+                    raise ValueError("Model name in invalid format: " + product_parts_names[0])
+                model_name = result.group(1)
+                # Use ModelName to find the path to the custom model file
+                model_path = joinpath(CUSTOM_MODEL_PATH, model_name + ".py")
+                if not os.path.isfile(model_path):
+                    raise ValueError("The model file {} doesn't exist".format(model_path))
+                kernel_module = custom.load_custom_kernel_module(model_path)
+                return modelinfo.make_model_info(kernel_module)
+            # Model is a core model
+            kernel_module = generate.load_kernel_module(product_parts_names[0])
+            return modelinfo.make_model_info(kernel_module)
+
+    model = None
+    if len(product_parts) > 1:
+        model = mixture.make_mixture_info(product_parts, operation='*')
+    if len(addition_parts) > 1:
+        if model is not None:
+            addition_parts.append(model)
+        model = mixture.make_mixture_info(addition_parts, operation='+')
+    return model
 
 
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision 573ffab723e1b8b0b093a98e7408853703cb65d4)
+++ sasmodels/generate.py	(revision 30b60d21e384fd9fda8e281e787ff2869925e10d)
@@ -210,5 +210,5 @@
 
 # Conversion from units defined in the parameter table for each model
-# to units displayed in the sphinx documentation. 
+# to units displayed in the sphinx documentation.
 # This section associates the unit with the macro to use to produce the LaTex
 # code.  The macro itself needs to be defined in sasmodels/doc/rst_prolog.
@@ -216,5 +216,5 @@
 # NOTE: there is an RST_PROLOG at the end of this file which is NOT
 # used for the bundled documentation. Still as long as we are defining the macros
-# in two places any new addition should define the macro in both places. 
+# in two places any new addition should define the macro in both places.
 RST_UNITS = {
     "Ang": "|Ang|",
@@ -898,9 +898,4 @@
 .. |cm^-3| replace:: cm\ :sup:`-3`
 .. |sr^-1| replace:: sr\ :sup:`-1`
-.. |P0| replace:: P\ :sub:`0`\
-
-.. |equiv| unicode:: U+2261
-.. |noteql| unicode:: U+2260
-.. |TM| unicode:: U+2122
 
 .. |cdot| unicode:: U+00B7
Index: sasmodels/mixture.py
===================================================================
--- sasmodels/mixture.py	(revision 6dc78e4926d9f0e39bbb55c42a8f11b256f166fd)
+++ sasmodels/mixture.py	(revision 765eb0ed12773d554fa60f4cd5952b5aa5da7917)
@@ -25,28 +25,87 @@
     pass
 
-def make_mixture_info(parts):
+def make_mixture_info(parts, operation='+'):
     # type: (List[ModelInfo]) -> ModelInfo
     """
     Create info block for mixture model.
     """
-    flatten = []
-    for part in parts:
-        if part.composition and part.composition[0] == 'mixture':
-            flatten.extend(part.composition[1])
-        else:
-            flatten.append(part)
-    parts = flatten
-
     # Build new parameter list
     combined_pars = []
-    demo = {}
-    for k, part in enumerate(parts):
+
+    model_num = 0
+    all_parts = copy(parts)
+    is_flat = False
+    while not is_flat:
+        is_flat = True
+        for part in all_parts:
+            if part.composition and part.composition[0] == 'mixture' and \
+                len(part.composition[1]) > 1:
+                all_parts += part.composition[1]
+                all_parts.remove(part)
+                is_flat = False
+
+    # When creating a mixture model that is a sum of product models (ie (1*2)+(3*4))
+    # the parameters for models 1 & 2 will be prefixed with A & B respectively,
+    # but so will the parameters for models 3 & 4. We need to rename models 3 & 4
+    # so that they are prefixed with C & D to avoid overlap of parameter names.
+    used_prefixes = []
+    for part in parts:
+        i = 0
+        if part.composition and part.composition[0] == 'mixture':
+            npars_list = [info.parameters.npars for info in part.composition[1]]
+            for npars in npars_list:
+                # List of params of one of the constituent models of part
+                submodel_pars = part.parameters.kernel_parameters[i:i+npars]
+                # Prefix of the constituent model
+                prefix = submodel_pars[0].name[0]
+                if prefix not in used_prefixes: # Haven't seen this prefix so far
+                    used_prefixes.append(prefix)
+                    i += npars
+                    continue
+                while prefix in used_prefixes:
+                    # This prefix has been already used, so change it to the
+                    # next letter that hasn't been used
+                    prefix = chr(ord(prefix) + 1)
+                used_prefixes.append(prefix)
+                prefix += "_"
+                # Update the parameters of this constituent model to use the
+                # new prefix
+                for par in submodel_pars:
+                    par.id = prefix + par.id[2:]
+                    par.name = prefix + par.name[2:]
+                    if par.length_control is not None:
+                        par.length_control = prefix + par.length_control[2:]
+                i += npars
+
+    for part in parts:
         # Parameter prefix per model, A_, B_, ...
         # Note that prefix must also be applied to id and length_control
         # to support vector parameters
-        prefix = chr(ord('A')+k) + '_'
-        scale =  Parameter(prefix+'scale', default=1.0,
-                           description="model intensity for " + part.name)
-        combined_pars.append(scale)
+        prefix = ''
+        if not part.composition:
+            # Model isn't a composition model, so it's parameters don't have a
+            # a prefix. Add the next available prefix
+            prefix = chr(ord('A')+len(used_prefixes))
+            used_prefixes.append(prefix)
+            prefix += '_'
+            
+        if operation == '+':
+            # If model is a sum model, each constituent model gets its own scale parameter
+            scale_prefix = prefix
+            if prefix == '' and part.operation == '*':
+                # `part` is a composition product model. Find the prefixes of
+                # it's parameters to form a new prefix for the scale, eg:
+                # a model with A*B*C will have ABC_scale
+                sub_prefixes = []
+                for param in part.parameters.kernel_parameters:
+                    # Prefix of constituent model
+                    sub_prefix = param.id.split('_')[0]
+                    if sub_prefix not in sub_prefixes:
+                        sub_prefixes.append(sub_prefix)
+                # Concatenate sub_prefixes to form prefix for the scale
+                scale_prefix = ''.join(sub_prefixes) + '_'
+            scale =  Parameter(scale_prefix + 'scale', default=1.0,
+                            description="model intensity for " + part.name)
+            combined_pars.append(scale)
         for p in part.parameters.kernel_parameters:
             p = copy(p)
@@ -56,13 +115,19 @@
                 p.length_control = prefix + p.length_control
             combined_pars.append(p)
-        demo.update((prefix+k, v) for k, v in part.demo.items()
-                    if k != "background")
-    #print("pars",combined_pars)
     parameters = ParameterTable(combined_pars)
     parameters.max_pd = sum(part.parameters.max_pd for part in parts)
 
+    def random():
+        combined_pars = {}
+        for k, part in enumerate(parts):
+            prefix = chr(ord('A')+k) + '_'
+            pars = part.random()
+            combined_pars.update((prefix+k, v) for k, v in pars.items())
+        return combined_pars
+
     model_info = ModelInfo()
-    model_info.id = '+'.join(part.id for part in parts)
-    model_info.name = ' + '.join(part.name for part in parts)
+    model_info.id = operation.join(part.id for part in parts)
+    model_info.operation = operation
+    model_info.name = '(' + operation.join(part.name for part in parts) + ')'
     model_info.filename = None
     model_info.title = 'Mixture model with ' + model_info.name
@@ -71,4 +136,5 @@
     model_info.category = "custom"
     model_info.parameters = parameters
+    model_info.random = random
     #model_info.single = any(part['single'] for part in parts)
     model_info.structure_factor = False
@@ -79,5 +145,4 @@
     # Remember the component info blocks so we can build the model
     model_info.composition = ('mixture', parts)
-    model_info.demo = demo
     return model_info
 
@@ -88,4 +153,5 @@
         self.info = model_info
         self.parts = parts
+        self.dtype = parts[0].dtype
 
     def make_kernel(self, q_vectors):
@@ -116,4 +182,5 @@
         self.kernels = kernels
         self.dtype = self.kernels[0].dtype
+        self.operation = model_info.operation
         self.results = []  # type: List[np.ndarray]
 
@@ -124,11 +191,17 @@
         # remember the parts for plotting later
         self.results = []  # type: List[np.ndarray]
-        offset = 2 # skip scale & background
         parts = MixtureParts(self.info, self.kernels, call_details, values)
         for kernel, kernel_details, kernel_values in parts:
             #print("calling kernel", kernel.info.name)
             result = kernel(kernel_details, kernel_values, cutoff, magnetic)
-            #print(kernel.info.name, result)
-            total += result
+            result = np.array(result).astype(kernel.dtype)
+            # print(kernel.info.name, result)
+            if self.operation == '+':
+                total += result
+            elif self.operation == '*':
+                if np.all(total) == 0.0:
+                    total = result
+                else:
+                    total *= result
             self.results.append(result)
 
@@ -171,5 +244,7 @@
 
         self.part_num += 1
-        self.par_index += info.parameters.npars + 1
+        self.par_index += info.parameters.npars
+        if self.model_info.operation == '+':
+            self.par_index += 1 # Account for each constituent model's scale param
         self.mag_index += 3 * len(info.parameters.magnetism_index)
 
@@ -182,9 +257,10 @@
         # which includes the initial scale and background parameters.
         # We want the index into the weight length/offset for each parameter.
-        # Exclude the initial scale and background, so subtract two, but each
-        # component has its own scale factor which we need to skip when
-        # constructing the details for the kernel, so add one, giving a
-        # net subtract one.
-        index = slice(par_index - 1, par_index - 1 + info.parameters.npars)
+        # Exclude the initial scale and background, so subtract two. If we're
+        # building an addition model, each component has its own scale factor
+        # which we need to skip when constructing the details for the kernel, so
+        # add one, giving a net subtract one.
+        diff = 1 if self.model_info.operation == '+' else 2
+        index = slice(par_index - diff, par_index - diff + info.parameters.npars)
         length = full.length[index]
         offset = full.offset[index]
@@ -196,7 +272,8 @@
     def _part_values(self, info, par_index, mag_index):
         # type: (ModelInfo, int, int) -> np.ndarray
-        #print(info.name, par_index, self.values[par_index:par_index + info.parameters.npars + 1])
-        scale = self.values[par_index]
-        pars = self.values[par_index + 1:par_index + info.parameters.npars + 1]
+        # Set each constituent model's scale to 1 if this is a multiplication model
+        scale = self.values[par_index] if self.model_info.operation == '+' else 1.0
+        diff = 1 if self.model_info.operation == '+' else 0 # Skip scale if addition model
+        pars = self.values[par_index + diff:par_index + info.parameters.npars + diff]
         nmagnetic = len(info.parameters.magnetism_index)
         if nmagnetic:
Index: sasmodels/model_test.py
===================================================================
--- sasmodels/model_test.py	(revision bedb9b0fed4f3f4bc2bbfa5878de6f2b6fdfbcc9)
+++ sasmodels/model_test.py	(revision 65314f75d3d74de18ada1d98440acf06ab8c4ce3)
@@ -201,6 +201,7 @@
                 ({}, 'VR', None),
                 ]
-
-            tests = smoke_tests + self.info.tests
+            tests = smoke_tests
+            if self.info.tests is not None:
+                tests += self.info.tests
             try:
                 model = build_model(self.info, dtype=self.dtype,
@@ -371,5 +372,4 @@
         stream.writeln(traceback.format_exc())
         return
-
     # Run the test suite
     suite.run(result)
Index: sasmodels/modelinfo.py
===================================================================
--- sasmodels/modelinfo.py	(revision 0bdddc2288f8fc5bc5e416af84e456e811c3b733)
+++ sasmodels/modelinfo.py	(revision 65314f75d3d74de18ada1d98440acf06ab8c4ce3)
@@ -727,4 +727,7 @@
     models when the model is first called, not when the model is loaded.
     """
+    if hasattr(kernel_module, "model_info"):
+        # Custom sum/multi models
+        return kernel_module.model_info
     info = ModelInfo()
     #print("make parameter table", kernel_module.parameters)
Index: sasmodels/models/binary_hard_sphere.py
===================================================================
--- sasmodels/models/binary_hard_sphere.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
+++ sasmodels/models/binary_hard_sphere.py	(revision 30b60d21e384fd9fda8e281e787ff2869925e10d)
@@ -23,10 +23,10 @@
     :nowrap:
 
-    \begin{align}
+    \begin{align*}
     x &= \frac{(\phi_2 / \phi)\alpha^3}{(1-(\phi_2/\phi) + (\phi_2/\phi)
     \alpha^3)} \\
     \phi &= \phi_1 + \phi_2 = \text{total volume fraction} \\
     \alpha &= R_1/R_2 = \text{size ratio}
-    \end{align}
+    \end{align*}
 
 The 2D scattering intensity is the same as 1D, regardless of the orientation of
Index: sasmodels/models/core_shell_bicelle.py
===================================================================
--- sasmodels/models/core_shell_bicelle.py	(revision a151caa16f645cbb8ff93c008adabc6d6a7f5136)
+++ sasmodels/models/core_shell_bicelle.py	(revision 30b60d21e384fd9fda8e281e787ff2869925e10d)
@@ -41,11 +41,12 @@
 
     I(Q,\alpha) = \frac{\text{scale}}{V_t} \cdot
-        F(Q,\alpha)^2.sin(\alpha) + \text{background}
+        F(Q,\alpha)^2 \cdot sin(\alpha) + \text{background}
 
 where
 
 .. math::
+    :nowrap:
 
-    \begin{align}
+    \begin{align*}
     F(Q,\alpha) = &\bigg[
     (\rho_c - \rho_f) V_c \frac{2J_1(QRsin \alpha)}{QRsin\alpha}\frac{sin(QLcos\alpha/2)}{Q(L/2)cos\alpha} \\
@@ -53,5 +54,5 @@
     &+(\rho_r - \rho_s) V_t \frac{2J_1(Q(R+t_r)sin\alpha)}{Q(R+t_r)sin\alpha}\frac{sin(Q(L/2+t_f)cos\alpha)}{Q(L/2+t_f)cos\alpha}
     \bigg]
-    \end{align}
+    \end{align*}
 
 where $V_t$ is the total volume of the bicelle, $V_c$ the volume of the core,
Index: sasmodels/models/core_shell_bicelle_elliptical.py
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
+++ sasmodels/models/core_shell_bicelle_elliptical.py	(revision 30b60d21e384fd9fda8e281e787ff2869925e10d)
@@ -42,11 +42,12 @@
 
     I(Q,\alpha,\psi) = \frac{\text{scale}}{V_t} \cdot
-        F(Q,\alpha, \psi)^2.sin(\alpha) + \text{background}
+        F(Q,\alpha, \psi)^2 \cdot sin(\alpha) + \text{background}
 
-where a numerical integration of $F(Q,\alpha, \psi)^2.sin(\alpha)$ is carried out over \alpha and \psi for:
+where a numerical integration of $F(Q,\alpha, \psi)^2 \cdot sin(\alpha)$ is carried out over \alpha and \psi for:
 
 .. math::
+    :nowrap:
 
-        \begin{align}
+    \begin{align*}
     F(Q,\alpha,\psi) = &\bigg[
     (\rho_c - \rho_f) V_c \frac{2J_1(QR'sin \alpha)}{QR'sin\alpha}\frac{sin(QLcos\alpha/2)}{Q(L/2)cos\alpha} \\
@@ -54,5 +55,5 @@
     &+(\rho_r - \rho_s) V_t \frac{2J_1(Q(R'+t_r)sin\alpha)}{Q(R'+t_r)sin\alpha}\frac{sin(Q(L/2+t_f)cos\alpha)}{Q(L/2+t_f)cos\alpha}
     \bigg]
-    \end{align}
+    \end{align*}
 
 where
Index: sasmodels/models/core_shell_ellipsoid.py
===================================================================
--- sasmodels/models/core_shell_ellipsoid.py	(revision 9f6823bee65fae551d3909e8bcc7dbb030b8d5ed)
+++ sasmodels/models/core_shell_ellipsoid.py	(revision 30b60d21e384fd9fda8e281e787ff2869925e10d)
@@ -44,8 +44,10 @@
 
 .. math::
-    \begin{align}
+    :nowrap:
+
+    \begin{align*}
     F(q,\alpha) = &f(q,radius\_equat\_core,radius\_equat\_core.x\_core,\alpha) \\
     &+ f(q,radius\_equat\_core + thick\_shell,radius\_equat\_core.x\_core + thick\_shell.x\_polar\_shell,\alpha)
-    \end{align}
+    \end{align*}
 
 where
Index: sasmodels/models/fractal_core_shell.py
===================================================================
--- sasmodels/models/fractal_core_shell.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
+++ sasmodels/models/fractal_core_shell.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -31,5 +31,5 @@
 $\rho_{solv}$ are the scattering length densities of the core, shell, and
 solvent respectively, $r_c$ and $r_s$ are the radius of the core and the radius
-of the whole particle respectively, $D_f$ is the fractal dimension, and |xi| the
+of the whole particle respectively, $D_f$ is the fractal dimension, and $\xi$ the
 correlation length.
 
Index: sasmodels/models/hollow_rectangular_prism.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
+++ sasmodels/models/hollow_rectangular_prism.py	(revision 30b60d21e384fd9fda8e281e787ff2869925e10d)
@@ -31,5 +31,5 @@
   :nowrap:
 
-  \begin{align}
+  \begin{align*}
   A_{P\Delta}(q) & =  A B C
     \left[\frac{\sin \bigl( q \frac{C}{2} \cos\theta \bigr)}
@@ -47,5 +47,5 @@
     \left[ \frac{\sin \bigl[ q \bigl(\frac{B}{2}-\Delta\bigr) \sin\theta \cos\phi \bigr]}
     {q \bigl(\frac{B}{2}-\Delta\bigr) \sin\theta \cos\phi} \right]
-  \end{align}
+  \end{align*}
 
 where $A$, $B$ and $C$ are the external sides of the parallelepiped fulfilling
Index: sasmodels/models/multilayer_vesicle.py
===================================================================
--- sasmodels/models/multilayer_vesicle.py	(revision 870a2f45e993d96ff75196552953ec259f2ba197)
+++ sasmodels/models/multilayer_vesicle.py	(revision 870a2f45e993d96ff75196552953ec259f2ba197)
@@ -33,6 +33,6 @@
 .. math::
 
-     r_i &= r_c + (i-1)(t_s + t_w) && \text{ solvent radius before shell } i \\
-     R_i &= r_i + t_s && \text{ shell radius for shell } i
+     r_i &= r_c + (i-1)(t_s + t_w) \text{ solvent radius before shell } i \\
+     R_i &= r_i + t_s \text{ shell radius for shell } i
 
 $\phi$ is the volume fraction of particles, $V(r)$ is the volume of a sphere
Index: sasmodels/models/parallelepiped.py
===================================================================
--- sasmodels/models/parallelepiped.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
+++ sasmodels/models/parallelepiped.py	(revision 30b60d21e384fd9fda8e281e787ff2869925e10d)
@@ -20,11 +20,9 @@
    Parallelepiped with the corresponding definition of sides.
 
-.. note::
-
-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 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:
Index: sasmodels/models/rpa.py
===================================================================
--- sasmodels/models/rpa.py	(revision 4f9e288196c233a7668a135cc94bb86a11da7958)
+++ sasmodels/models/rpa.py	(revision 30b60d21e384fd9fda8e281e787ff2869925e10d)
@@ -30,10 +30,10 @@
     These case numbers are different from those in the NIST SANS package!
 
-The models are based on the papers by Akcasu et al. [#Akcasu]_ and by
-Hammouda [#Hammouda]_ assuming the polymer follows Gaussian statistics such
+The models are based on the papers by Akcasu *et al.* and by
+Hammouda assuming the polymer follows Gaussian statistics such
 that $R_g^2 = n b^2/6$ where $b$ is the statistical segment length and $n$ is
 the number of statistical segment lengths. A nice tutorial on how these are
 constructed and implemented can be found in chapters 28 and 39 of Boualem
-Hammouda's 'SANS Toolbox'[#toolbox]_.
+Hammouda's 'SANS Toolbox'.
 
 In brief the macroscopic cross sections are derived from the general forms
@@ -49,5 +49,5 @@
   are calculated with respect to component D).** So the scattering contrast
   for a C/D blend = [SLD(component C) - SLD(component D)]\ :sup:`2`.
-* Depending on which case is being used, the number of fitting parameters can 
+* Depending on which case is being used, the number of fitting parameters can
   vary.
 
@@ -57,15 +57,16 @@
       component are obtained from other methods and held fixed while The *scale*
       parameter should be held equal to unity.
-    * The variables are normally the segment lengths (b\ :sub:`a`, b\ :sub:`b`,
-      etc) and $\chi$ parameters (K\ :sub:`ab`, K\ :sub:`ac`, etc).
-
+    * The variables are normally the segment lengths ($b_a$, $b_b$,
+      etc.) and $\chi$ parameters ($K_{ab}$, $K_{ac}$, etc).
 
 References
 ----------
 
-.. [#Akcasu] A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993)
-   4136.
-.. [#Hammouda] B. Hammouda, *Advances in Polymer Science* 106 (1993) 87.
-.. [#toolbox] https://www.ncnr.nist.gov/staff/hammouda/the_sans_toolbox.pdf
+A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136.
+
+B. Hammouda, *Advances in Polymer Science* 106 (1993) 87.
+
+B. Hammouda, *SANS Toolbox*
+https://www.ncnr.nist.gov/staff/hammouda/the_sans_toolbox.pdf.
 
 Authorship and Verification
Index: sasmodels/models/spherical_sld.py
===================================================================
--- sasmodels/models/spherical_sld.py	(revision 63a7fe827229bb42f3d028392a52333e778ea890)
+++ sasmodels/models/spherical_sld.py	(revision 2ad5d3007889c5a943ccedaeada465c4be51317a)
@@ -18,9 +18,9 @@
 Interface shapes are as follows::
 
-    0: erf(|nu|*z)
-    1: Rpow(z^|nu|)
-    2: Lpow(z^|nu|)
-    3: Rexp(-|nu|z)
-    4: Lexp(-|nu|z)
+    0: erf($\nu z$)
+    1: Rpow($z^\nu$)
+    2: Lpow($z^\nu$)
+    3: Rexp($-\nu z$)
+    4: Lexp($-\nu z$)
 
 Definition
Index: sasmodels/models/star_polymer.py
===================================================================
--- sasmodels/models/star_polymer.py	(revision 5da1ac80e764ac02d7efd2aa7500dac229b0754b)
+++ sasmodels/models/star_polymer.py	(revision 5da1ac80e764ac02d7efd2aa7500dac229b0754b)
@@ -7,5 +7,5 @@
 emanating from a common central (in the case of this model) point.  It is
 derived as a special case of on the Benoit model for general branched
-polymers\ [#CITBenoit]_ as also used by Richter ''et. al.''\ [#CITRichter]_
+polymers\ [#CITBenoit]_ as also used by Richter *et al.*\ [#CITRichter]_
 
 For a star with $f$ arms the scattering intensity $I(q)$ is calculated as
Index: sasmodels/models/surface_fractal.py
===================================================================
--- sasmodels/models/surface_fractal.py	(revision 48462b091fcc67fde1f1d038e41eb58b73077f65)
+++ sasmodels/models/surface_fractal.py	(revision 30b60d21e384fd9fda8e281e787ff2869925e10d)
@@ -9,5 +9,7 @@
 
 .. math::
+    :nowrap:
 
+    \begin{align*}
     I(q) &= \text{scale} \times P(q)S(q) + \text{background} \\
     P(q) &= F(qR)^2 \\
@@ -15,6 +17,7 @@
     S(q) &= \Gamma(5-D_S)\xi^{\,5-D_S}\left[1+(q\xi)^2 \right]^{-(5-D_S)/2}
             \sin\left[-(5-D_S) \tan^{-1}(q\xi) \right] q^{-1} \\
-    \text{scale} &= \text{scale_factor}\, N V^2(\rho_\text{particle} - \rho_\text{solvent})^2 \\
+    \text{scale} &= \text{scale factor}\, N V^1(\rho_\text{particle} - \rho_\text{solvent})^2 \\
     V &= \frac{4}{3}\pi R^3
+    \end{align*}
 
 where $R$ is the radius of the building block, $D_S$ is the **surface** fractal
Index: sasmodels/product.py
===================================================================
--- sasmodels/product.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
+++ sasmodels/product.py	(revision 765eb0ed12773d554fa60f4cd5952b5aa5da7917)
@@ -68,15 +68,19 @@
     translate_name = dict((old.id, new.id) for old, new
                           in zip(s_pars.kernel_parameters[1:], s_list))
-    demo = {}
-    demo.update(p_info.demo.items())
-    demo.update((translate_name[k], v) for k, v in s_info.demo.items()
-                if k not in ("background", "scale") and not k.startswith(ER_ID))
     combined_pars = p_pars.kernel_parameters + s_list
     parameters = ParameterTable(combined_pars)
     parameters.max_pd = p_pars.max_pd + s_pars.max_pd
+    def random():
+        combined_pars = p_info.random()
+        s_names = set(par.id for par in s_pars.kernel_parameters[1:])
+        s = s_info.random()
+        combined_pars.update((translate_name[k], v)
+                    for k, v in s_info.random().items()
+                    if k in s_names)
+        return combined_pars
 
     model_info = ModelInfo()
-    model_info.id = '*'.join((p_id, s_id))
-    model_info.name = '*'.join((p_name, s_name))
+    model_info.id = '@'.join((p_id, s_id))
+    model_info.name = '@'.join((p_name, s_name))
     model_info.filename = None
     model_info.title = 'Product of %s and %s'%(p_name, s_name)
@@ -85,4 +89,5 @@
     model_info.category = "custom"
     model_info.parameters = parameters
+    model_info.random = random
     #model_info.single = p_info.single and s_info.single
     model_info.structure_factor = False
@@ -95,10 +100,9 @@
     # TODO: delegate random to p_info, s_info
     #model_info.random = lambda: {}
-    model_info.demo = demo
-
-    ## Show the parameter table with the demo values
+
+    ## Show the parameter table
     #from .compare import get_pars, parlist
     #print("==== %s ====="%model_info.name)
-    #values = get_pars(model_info, use_demo=True)
+    #values = get_pars(model_info)
     #print(parlist(model_info, values, is2d=True))
     return model_info
@@ -126,4 +130,5 @@
         self.P = P
         self.S = S
+        self.dtype = P.dtype
 
     def make_kernel(self, q_vectors):
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision 724257c92ca5f608f1dc18211ee76a0f7d2fe6f1)
+++ sasmodels/sasview_model.py	(revision bcdd6c9cc1061cc25272a94d1560f23cf0e31659)
@@ -120,5 +120,5 @@
     else:
         model_info = modelinfo.make_model_info(kernel_module)
-        model = _make_model_from_info(model_info)
+        model = make_model_from_info(model_info)
     model.timestamp = getmtime(path)
 
@@ -142,4 +142,18 @@
 
 
+def make_model_from_info(model_info):
+    # type: (ModelInfo) -> SasviewModelType
+    """
+    Convert *model_info* into a SasView model wrapper.
+    """
+    def __init__(self, multiplicity=None):
+        SasviewModel.__init__(self, multiplicity=multiplicity)
+    attrs = _generate_model_attributes(model_info)
+    attrs['__init__'] = __init__
+    attrs['filename'] = model_info.filename
+    ConstructedModel = type(model_info.name, (SasviewModel,), attrs) # type: SasviewModelType
+    return ConstructedModel
+
+
 def _make_standard_model(name):
     # type: (str) -> SasviewModelType
@@ -153,7 +167,7 @@
     kernel_module = generate.load_kernel_module(name)
     model_info = modelinfo.make_model_info(kernel_module)
-    return _make_model_from_info(model_info)
-
-
+    return make_model_from_info(model_info)
+
+    
 def _register_old_models():
     # type: () -> None
@@ -187,19 +201,7 @@
     model_info = product.make_product_info(form_factor._model_info,
                                            structure_factor._model_info)
-    ConstructedModel = _make_model_from_info(model_info)
+    ConstructedModel = make_model_from_info(model_info)
     return ConstructedModel()
 
-def _make_model_from_info(model_info):
-    # type: (ModelInfo) -> SasviewModelType
-    """
-    Convert *model_info* into a SasView model wrapper.
-    """
-    def __init__(self, multiplicity=None):
-        SasviewModel.__init__(self, multiplicity=multiplicity)
-    attrs = _generate_model_attributes(model_info)
-    attrs['__init__'] = __init__
-    attrs['filename'] = model_info.filename
-    ConstructedModel = type(model_info.name, (SasviewModel,), attrs) # type: SasviewModelType
-    return ConstructedModel
 
 def _generate_model_attributes(model_info):
@@ -603,6 +605,6 @@
         if hasattr(self._model_info, "composition") \
            and self._model_info.composition is not None:
-            p_model = _make_model_from_info(self._model_info.composition[1][0])()
-            s_model = _make_model_from_info(self._model_info.composition[1][1])()
+            p_model = make_model_from_info(self._model_info.composition[1][0])()
+            s_model = make_model_from_info(self._model_info.composition[1][1])()
         return p_model, s_model
 
