Index: doc/guide/pd/polydispersity.rst
===================================================================
--- doc/guide/pd/polydispersity.rst	(revision d712a0f4699e15d422e39b923be41331fe3922bf)
+++ doc/guide/pd/polydispersity.rst	(revision 9ce5bcbf7c22686c65d77b820be6660414e5a35b)
@@ -33,12 +33,12 @@
 sigmas $N_\sigma$ to include from the tails of the distribution, and the
 number of points used to compute the average. The center of the distribution
-is set by the value of the model parameter. The meaning of a polydispersity 
-parameter *PD* (not to be confused with a molecular weight distributions 
-in polymer science) in a model depends on the type of parameter it is being 
+is set by the value of the model parameter. The meaning of a polydispersity
+parameter *PD* (not to be confused with a molecular weight distributions
+in polymer science) in a model depends on the type of parameter it is being
 applied too.
 
-The distribution width applied to *volume* (ie, shape-describing) parameters 
-is relative to the center value such that $\sigma = \mathrm{PD} \cdot \bar x$. 
-However, the distribution width applied to *orientation* (ie, angle-describing) 
+The distribution width applied to *volume* (ie, shape-describing) parameters
+is relative to the center value such that $\sigma = \mathrm{PD} \cdot \bar x$.
+However, the distribution width applied to *orientation* (ie, angle-describing)
 parameters is just $\sigma = \mathrm{PD}$.
 
@@ -64,8 +64,8 @@
 *  *Schulz Distribution*
 *  *Array Distribution*
+*  *User-defined Distributions*
 
 These are all implemented as *number-average* distributions.
 
-Additional distributions are under consideration.
 
 .. note:: In 2009 IUPAC decided to introduce the new term 'dispersity' to replace 
@@ -87,9 +87,11 @@
 or angular orientations, use the Gaussian or Boltzmann distributions.
 
-If applying polydispersion to parameters describing angles, use the Uniform 
-distribution. Beware of using distributions that are always positive (eg, the 
+If applying polydispersion to parameters describing angles, use the Uniform
+distribution. Beware of using distributions that are always positive (eg, the
 Lognormal) because angles can be negative!
 
-The array distribution allows a user-defined distribution to be applied.
+The array distribution provides a very simple means of implementing a user-
+defined distribution, but without any fittable parameters. Greater flexibility
+is conferred by the user-defined distribution. 
 
 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
@@ -329,4 +331,84 @@
 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
 
+User-defined Distributions
+^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+You can also define your own distribution by creating a python file defining a
+*Distribution* object with a *_weights* method.  The *_weights* method takes
+*center*, *sigma*, *lb* and *ub* as arguments, and can access *self.npts*
+and *self.nsigmas* from the distribution.  They are interpreted as follows:
+
+* *center* the value of the shape parameter (for size dispersity) or zero
+  if it is an angular dispersity.  This parameter may be fitted.
+
+* *sigma* the width of the distribution, which is the polydispersity parameter
+  times the center for size dispersity, or the polydispersity parameter alone
+  for angular dispersity.  This parameter may be fitted.
+
+* *lb*, *ub* are the parameter limits (lower & upper bounds) given in the model
+  definition file.  For example, a radius parameter has *lb* equal to zero.  A
+  volume fraction parameter would have *lb* equal to zero and *ub* equal to one.
+
+* *self.nsigmas* the distance to go into the tails when evaluating the
+  distribution.  For a two parameter distribution, this value could be
+  co-opted to use for the second parameter, though it will not be available
+  for fitting.
+
+* *self.npts* the number of points to use when evaluating the distribution.
+  The user will adjust this to trade calculation time for accuracy, but the
+  distribution code is free to return more or fewer, or use it for the third
+  parameter in a three parameter distribution.
+
+As an example, the code following wraps the Laplace distribution from scipy stats::
+
+    import numpy as np
+    from scipy.stats import laplace
+
+    from sasmodels import weights
+
+    class Dispersion(weights.Dispersion):
+        r"""
+        Laplace distribution
+
+        .. math::
+
+            w(x) = e^{-\sigma |x - \mu|}
+        """
+        type = "laplace"
+        default = dict(npts=35, width=0, nsigmas=3)  # default values
+        def _weights(self, center, sigma, lb, ub):
+            x = self._linspace(center, sigma, lb, ub)
+            wx = laplace.pdf(x, center, sigma)
+            return x, wx
+
+You can plot the weights for a given value and width using the following::
+
+    from numpy import inf
+    from matplotlib import pyplot as plt
+    from sasmodels import weights
+
+    # reload the user-defined weights
+    weights.load_weights()
+    x, wx = weights.get_weights('laplace', n=35, width=0.1, nsigmas=3, value=50,
+                                limits=[0, inf], relative=True)
+
+    # plot the weights
+    plt.interactive(True)
+    plt.plot(x, wx, 'x')
+
+The *self.nsigmas* and *self.npts* parameters are normally used to control
+the accuracy of the distribution integral. The *self._linspace* function
+uses them to define the *x* values (along with the *center*, *sigma*,
+*lb*, and *ub* which are passed as parameters).  If you repurpose npts or
+nsigmas you will need to generate your own *x*.  Be sure to honour the
+limits *lb* and *ub*, for example to disallow a negative radius or constrain
+the volume fraction to lie between zero and one.
+
+To activate a user-defined distribution, set the following environment variable:
+
+    SASMODELS_WEIGHTS=path/to/folder/name_of_distribution.py
+
+.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
+
 Note about DLS polydispersity
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Index: example/weights/cyclic_gaussian.py
===================================================================
--- example/weights/cyclic_gaussian.py	(revision 36873d13d5ba83beb0b0f046deab7473e08aec25)
+++ example/weights/cyclic_gaussian.py	(revision 36873d13d5ba83beb0b0f046deab7473e08aec25)
@@ -0,0 +1,42 @@
+import numpy as np
+from numpy import exp, sin, cos, pi, radians, degrees
+
+from sasmodels.weights import Dispersion as BaseDispersion
+
+class Dispersion(BaseDispersion):
+    r"""
+    Cyclic gaussian dispersion on orientation.
+    
+    .. math:
+
+        w(\theta) = e^{-\frac{\sin^2 \theta}{2 \sigma^2}}
+
+    This provides a close match to the gaussian distribution for
+    low angles (with $\sin \theta \approx \theta$), but the tails 
+    are limited to $\pm 90^\circ$.  For $\sigma$ large the
+    distribution is approximately uniform.  The usual polar coordinate
+    projection applies, with $\theta$ weights scaled by $\cos \theta$
+    and $\phi$ weights unscaled.
+
+    This is closely related to a Maier-Saupe distribution with order
+    parameter $P_2$ and appropriate scaling constants, and changes
+    between $\sin$ and $\cos$ as appropriate for the coordinate system
+    representation.
+    """
+    type = "cyclic_gaussian"
+    default = dict(npts=35, width=1, nsigmas=3)
+
+    # Note: center is always zero for orientation distributions
+    def _weights(self, center, sigma, lb, ub):
+        # Convert sigma in degrees to the approximately equivalent Maier-Saupe "a"
+        sigma = radians(sigma)
+        a = -0.5/sigma**2
+
+        # Limit width to +/-90 degrees; use an open interval since the
+        # pattern at +90 is the same as that at -90.
+        width = min(self.nsigmas*sigma, pi/2)
+        x = np.linspace(-width, width, self.npts+2)[1:-1]
+
+        # Return orientation in degrees with Maier-Saupe weights
+        return degrees(x), exp(a*sin(x)**2)
+
Index: example/weights/laplace.py
===================================================================
--- example/weights/laplace.py	(revision 342b3dd82d82df44b7021cc1cd6ee52766afdf36)
+++ example/weights/laplace.py	(revision 342b3dd82d82df44b7021cc1cd6ee52766afdf36)
@@ -0,0 +1,20 @@
+import numpy as np
+from scipy.stats import laplace
+
+from sasmodels import weights
+
+class Dispersion(weights.Dispersion):
+    r"""
+    Laplace distribution
+
+    .. math::
+
+        w(x) = e^{-\sigma |x - \mu|}
+    """
+    type = "laplace"
+    default = dict(npts=35, width=0, nsigmas=3)  # default values
+    def _weights(self, center, sigma, lb, ub):
+        x = self._linspace(center, sigma, lb, ub)
+        wx = laplace.pdf(x, center, sigma)
+        return x, wx
+
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision af7a97c2acc05ec1c5d537fa21912dab44ff661a)
+++ sasmodels/compare.py	(revision af7a97c2acc05ec1c5d537fa21912dab44ff661a)
@@ -41,8 +41,8 @@
 from . import kerneldll
 from . import kernelcl
+from . import weights
 from .data import plot_theory, empty_data1D, empty_data2D, load_data
 from .direct_model import DirectModel, get_mesh
 from .generate import FLOAT_RE, set_integration_size
-from .weights import plot_weights
 
 # pylint: disable=unused-import
@@ -783,5 +783,5 @@
             model_info = base._kernel.info
             dim = base._kernel.dim
-            plot_weights(model_info, get_mesh(model_info, base_pars, dim=dim))
+            weights.plot_weights(model_info, get_mesh(model_info, base_pars, dim=dim))
         if opts['show_profile']:
             import pylab
@@ -1439,4 +1439,12 @@
     #import pprint; pprint.pprint(model_info)
 
+    # Hack to load user-defined distributions; run through all parameters
+    # and make sure any pd_type parameter is a defined distribution.
+    if (any(p.endswith('pd_type') and v not in weights.MODELS
+            for p, v in pars.items()) or
+        any(p.endswith('pd_type') and v not in weights.MODELS
+            for p, v in pars2.items())):
+       weights.load_weights()
+
     if opts['show_pars']:
         if model_info.name != model_info2.name or pars != pars2:
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision d533590d9857d7aef962616cfdaa6ecfc2e6751d)
+++ sasmodels/sasview_model.py	(revision aa25fc74c6ba4d068e48cae523a18e719c2e5ba7)
@@ -30,4 +30,9 @@
 from . import modelinfo
 from .details import make_kernel_args, dispersion_mesh
+
+# Hack: load in any custom distributions
+# Uses ~/.sasview/weights/*.py unless SASMODELS_WEIGHTS is set in the environ.
+# Override with weights.load_weights(pattern="<weights_path>/*.py")
+weights.load_weights()
 
 # pylint: disable=unused-import
Index: sasmodels/weights.py
===================================================================
--- sasmodels/weights.py	(revision 3d5824792e80506c5c40f05fb2ddc2234a8acedf)
+++ sasmodels/weights.py	(revision 910c0f46f74b74ccfc5001d660712566db4813eb)
@@ -231,4 +231,25 @@
 ))
 
+SASMODELS_WEIGHTS = "~/.sasview/weights/*.py"
+def load_weights(pattern=None):
+    # type: (str) -> None
+    """
+    Load dispersion distributions matching the given pattern
+    """
+    import logging
+    import os
+    import os.path
+    import glob
+    import traceback
+    from .custom import load_custom_kernel_module
+    if pattern is None:
+        pattern = os.environ.get("SASMODELS_WEIGHTS", SASMODELS_WEIGHTS)
+    for filename in sorted(glob.glob(os.path.expanduser(pattern))):
+        try:
+            #print("loading weights from", filename)
+            module = load_custom_kernel_module(filename)
+            MODELS[module.Dispersion.type] = module.Dispersion
+        except Exception as exc:
+            logging.error(traceback.format_exc(exc))
 
 def get_weights(disperser, n, width, nsigmas, value, limits, relative):
