Index: doc/guide/fitting_sq.rst
===================================================================
--- doc/guide/fitting_sq.rst	(revision 694c6d05ce7883894bb6e42d0f49b330f5bbb2c6)
+++ doc/guide/fitting_sq.rst	(revision 34483013ffb02edbec666039a8c85e18f605b690)
@@ -18,9 +18,9 @@
 of the form factor. For the theory behind this, see :ref:`PStheory` later.
 
-**If writing your own** $P@S$ **models, DO NOT give your model parameters**
-**these names!**
-
 Parameters
 ^^^^^^^^^^
+
+**Except for volfraction, when writing your own** $P@S$ **models, DO NOT give**
+**your model parameters these names!**
 
 Many parameters are common amongst $P@S$ models, but take on specific meanings:
@@ -37,6 +37,7 @@
     Structure factor models $S(Q)$ contain **volfraction**. In $P@S$ models
     this is *also* used as the volume fraction for the form factor model
-    $P(Q)$, *replacing* any **volfraction** parameter in $P(Q)$. This means
-    that $P@S$ models can also leave **scale** at 1.0.
+    $P(Q)$, so these models can also leave **scale** at 1.0.  If $P(Q)$ already
+    has a **volfraction** parameter, it is tied to the **volfraction** for
+    $S(Q)$.
 
     If the volume fraction required for $S(Q)$ is *not* the volume fraction
@@ -71,9 +72,9 @@
     length.
 
-    In use, it may be sensible to tie or constrain **radius_effective**
-    to one or other of the "size" parameters describing the form of the shape.
-
-    **radius_effective** may also be specified directly, independent of the
-    estimate from $P(Q)$.
+    If **radius_effective_mode = 0** (see below) it may be sensible to tie or
+    constrain **radius_effective** to one or other of the "size" parameters
+    describing the form of the shape (although the parameter cannot then be
+    polydisperse). But **radius_effective** may also be specified directly,
+    independent of the estimate from $P(Q)$.
 
     If **radius_effective** is calculated by $P(Q)$, it will be the
@@ -86,5 +87,5 @@
 .. note::
 
-   The following additional parameters are only available in SasView 5.0 and
+   The following additional parameters are only available in SasView 4.3 and
    later.
 
@@ -94,9 +95,14 @@
     be computed from the parameters of the shape.
 
-    When **radius_effective_mode = 0** then unconstrained **radius_effective**
-    parameter in the $S(Q)$ model is used. *This is the default in SasView*
-    *versions 4.x and earlier*. Otherwise, in SasView 5.0 and later,
-    **radius_effective_mode = k** represents an index in a list of alternative
-    **radius_effective** calculations which will appear in a drop-down box.
+    When **radius_effective_mode = 0** then the unconstrained
+    **radius_effective** parameter in the $S(Q)$ model is used. *This is the*
+    *default in SasView versions 4.x and earlier*. Otherwise, in SasView 4.3
+    and later, **radius_effective_mode = k** represents an index in a list of
+    alternative **radius_effective** calculations.
+
+    In SasView 4.3 and later **k** must be entered as an integer (and it will
+    be necessary to read the source code file to discover what calculations the
+    modes represent), but in SasView 5.0 and later the options appear in a
+    drop-down box.
 
     For example, the *ellipsoid* model defines the following
@@ -144,5 +150,5 @@
     oscillations in the normal (local monodisperse) $S(Q)$. When $\beta(Q) = 1$
     the local monodisperse approximation is recovered. *This mode is only*
-    *available in SasView 5.0 and later*.
+    *available in SasView 4.3 and later*.
 
     More mode options may appear in future as more complicated operations are
Index: doc/guide/pd/polydispersity.rst
===================================================================
--- doc/guide/pd/polydispersity.rst	(revision d089a00eb2e0e6bbb032f5e258e960009a322466)
+++ doc/guide/pd/polydispersity.rst	(revision a5cb9bcbc504127c414602837e600378ebda6add)
@@ -11,11 +11,11 @@
 --------------------------------------------
 
-For some models we can calculate the average intensity for a population of 
-particles that possess size and/or orientational (ie, angular) distributions. 
-In SasView we call the former *polydispersity* but use the parameter *PD* to 
-parameterise both. In other words, the meaning of *PD* in a model depends on 
+For some models we can calculate the average intensity for a population of
+particles that possess size and/or orientational (ie, angular) distributions.
+In SasView we call the former *polydispersity* but use the parameter *PD* to
+parameterise both. In other words, the meaning of *PD* in a model depends on
 the actual parameter it is being applied too.
 
-The resultant intensity is then normalized by the average particle volume such 
+The resultant intensity is then normalized by the average particle volume such
 that
 
@@ -24,23 +24,26 @@
   P(q) = \text{scale} \langle F^* F \rangle / V + \text{background}
 
-where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an 
+where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an
 average over the distribution $f(x; \bar x, \sigma)$, giving
 
 .. math::
 
-  P(q) = \frac{\text{scale}}{V} \int_\mathbb{R} 
+  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
 $x_\text{med}$, a width parameter $\sigma$ (note this is *not necessarily*
-the standard deviation, so read the description of the distribution carefully), 
-the number of 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 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* parameters is just 
-$\sigma = \mathrm{PD}$.
+the standard deviation, so read the description carefully), the number of
+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
+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)
+parameters is just $\sigma = \mathrm{PD}$.
 
 $N_\sigma$ determines how far into the tails to evaluate the distribution,
@@ -52,7 +55,7 @@
 
 Users should note that the averaging computation is very intensive. Applying
-polydispersion and/or orientational distributions to multiple parameters at 
-the same time, or increasing the number of points in the distribution, will 
-require patience! However, the calculations are generally more robust with 
+polydispersion and/or orientational distributions to multiple parameters at
+the same time, or increasing the number of points in the distribution, will
+require patience! However, the calculations are generally more robust with
 more data points or more angles.
 
@@ -66,8 +69,8 @@
 *  *Schulz Distribution*
 *  *Array Distribution*
+*  *User-defined Distributions*
 
 These are all implemented as *number-average* distributions.
 
-Additional distributions are under consideration.
 
 **Beware: when the Polydispersity & Orientational Distribution panel in SasView is**
@@ -75,10 +78,10 @@
 **This may not be suitable. See Suggested Applications below.**
 
-.. 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 chemical 
-           properties unambiguous. However, these terms are unrelated to the 
-           proportional size distributions and orientational distributions used in 
+.. 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 chemical
+           properties unambiguous. However, these terms are unrelated to the
+           proportional size distributions and orientational distributions used in
            SasView models.
 
@@ -92,9 +95,11 @@
 or angular orientations, consider using 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
@@ -334,14 +339,102 @@
 .. 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, put it in a file such as *distname.py*
+in the *SAS_WEIGHTS_PATH* folder.  This is defined with an environment
+variable, defaulting to::
+
+    SAS_WEIGHTS_PATH=~/.sasview/weights
+
+The weights path is loaded on startup.  To update the distribution definition
+in a running application you will need to enter the following python commands::
+
+    import sasmodels.weights
+    sasmodels.weights.load_weights('path/to/distname.py')
+
+.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
+
 Note about DLS polydispersity
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-Several measures of polydispersity abound in Dynamic Light Scattering (DLS) and 
-it should not be assumed that any of the following can be simply equated with 
+Several measures of polydispersity abound in Dynamic Light Scattering (DLS) and
+it should not be assumed that any of the following can be simply equated with
 the polydispersity *PD* parameter used in SasView.
 
-The dimensionless **Polydispersity Index (PI)** is a measure of the width of the 
-distribution of autocorrelation function decay rates (*not* the distribution of 
-particle sizes itself, though the two are inversely related) and is defined by 
+The dimensionless **Polydispersity Index (PI)** is a measure of the width of the
+distribution of autocorrelation function decay rates (*not* the distribution of
+particle sizes itself, though the two are inversely related) and is defined by
 ISO 22412:2017 as
 
@@ -350,5 +443,5 @@
     PI = \mu_{2} / \bar \Gamma^2
 
-where $\mu_\text{2}$ is the second cumulant, and $\bar \Gamma^2$ is the 
+where $\mu_\text{2}$ is the second cumulant, and $\bar \Gamma^2$ is the
 intensity-weighted average value, of the distribution of decay rates.
 
@@ -359,5 +452,5 @@
     PI = \sigma^2 / 2\bar \Gamma^2
 
-where $\sigma$ is the standard deviation, allowing a **Relative Polydispersity (RP)** 
+where $\sigma$ is the standard deviation, allowing a **Relative Polydispersity (RP)**
 to be defined as
 
@@ -366,9 +459,9 @@
     RP = \sigma / \bar \Gamma = \sqrt{2 \cdot PI}
 
-PI values smaller than 0.05 indicate a highly monodisperse system. Values 
+PI values smaller than 0.05 indicate a highly monodisperse system. Values
 greater than 0.7 indicate significant polydispersity.
 
-The **size polydispersity P-parameter** is defined as the relative standard 
-deviation coefficient of variation  
+The **size polydispersity P-parameter** is defined as the relative standard
+deviation coefficient of variation
 
 .. math::
@@ -377,5 +470,5 @@
 
 where $\nu$ is the variance of the distribution and $\bar R$ is the mean
-value of $R$. Here, the product $P \bar R$ is *equal* to the standard 
+value of $R$. Here, the product $P \bar R$ is *equal* to the standard
 deviation of the Lognormal distribution.
 
Index: example/weights/cyclic_gaussian.py
===================================================================
--- example/weights/cyclic_gaussian.py	(revision 35d230086190fc31926c237440d996e7ce2c269b)
+++ example/weights/cyclic_gaussian.py	(revision 35d230086190fc31926c237440d996e7ce2c269b)
@@ -0,0 +1,39 @@
+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, 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 eqivalent to a Maier-Saupe distribution with order
+    parameter $a = 1/(2 \sigma^2)$, with $\sigma$ in radians.
+    """
+    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 radians
+        sigma = radians(sigma)
+
+        # Limit width to +/- 90 degrees
+        width = min(self.nsigmas*sigma, pi/2)
+        x = np.linspace(-width, width, self.npts)
+
+        # Truncate the distribution in case the parameter value is limited
+        x[(x >= radians(lb)) & (x <= radians(ub))]
+
+        # Return orientation in degrees with Maier-Saupe weights
+        return degrees(x), exp(-0.5*sin(x)**2/sigma**2)
Index: example/weights/gaussian_eq.py
===================================================================
--- example/weights/gaussian_eq.py	(revision 55e82f04d46e375f4f5aa858d1291b28da4d6219)
+++ example/weights/gaussian_eq.py	(revision 55e82f04d46e375f4f5aa858d1291b28da4d6219)
@@ -0,0 +1,42 @@
+import numpy as np
+import scipy.stats
+
+from sasmodels.weights import Dispersion as BaseDispersion
+
+class Dispersion(BaseDispersion):
+    r"""
+    Gaussian dispersion, with 1-$\sigma$ width.
+
+    .. math::
+
+        w = \exp\left(-\tfrac12 (x - c)^2/\sigma^2\right)
+
+    $x$ points are chosen such that each interval has equal weight.
+
+    This works surprisingly poorly.  Try::
+
+        $ sascomp cylinder -2d theta=45 phi=20 phi_pd_type=gaussian_eq \
+          phi_pd_n=100,1000 radius=50 length=2*radius -midq phi_pd=5
+
+    Leaving it here for others to improve.
+    """
+    type = "gaussian_eq"
+    default = dict(npts=35, width=0, nsigmas=3)
+
+    def _weights(self, center, sigma, lb, ub):
+        # Use the gaussian distribution from scipy.stats
+        dist = scipy.stats.norm(center, sigma)
+
+        # Find the mid-points of the cdf intervals
+        cdf = np.linspace(0, 1, self.npts+2)[1:-1]
+        x = dist.ppf(cdf)
+
+        # Since we are equally spaced in cdf, all weights are the same
+        wx = np.ones_like(x)
+
+        # Truncate the distribution in case the parameter value is limited
+        index = (x >= lb) & (x <= ub)
+        x, wx = x[index], wx[index]
+
+        return x, wx
+
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: example/weights/maier_saupe.py
===================================================================
--- example/weights/maier_saupe.py	(revision a5516b1b51923e98e6773ae919d62f40f1c6e567)
+++ example/weights/maier_saupe.py	(revision a5516b1b51923e98e6773ae919d62f40f1c6e567)
@@ -0,0 +1,132 @@
+from __future__ import print_function
+import numpy as np
+from numpy import exp, sin, degrees, radians, pi, sqrt
+
+from sasmodels.weights import Dispersion as BaseDispersion
+
+class Dispersion(BaseDispersion):
+    r"""
+    Maier-Saupe dispersion on orientation.
+
+    .. math:
+
+        w(\theta) = e^{a{\cos^2 \theta}}
+
+    This provides a close match to the gaussian distribution for
+    low angles, but the tails are limited to $\pm 90^\circ$.  For $a \ll 1$
+    the distribution is approximately uniform.  The usual polar coordinate
+    projection applies, with $\theta$ weights scaled by $\cos \theta$
+    and $\phi$ weights unscaled.
+
+    This is equivalent to a cyclic gaussian distribution
+    $w(\theta) = e^{-sin^2(\theta)/(2\a^2)}.
+
+    Note that the incorrect distribution is used for $a=0$.  With the
+    distribution parameter labelled as *width*, the value *width=0* is
+    is assumed to be completely oriented and no polydispersity distribution
+    is generated.  However a value of $a=0$ should be completely unoriented.
+    Any small value (e.g., 0.01 or lower) should suffice.
+
+    The order parameter $P_2$ is defined as
+
+    .. math:
+
+        P(a, \beta) = \frac{e^{a \cos^2 \beta}}{
+                       4\pi \int_0^{\pi/2} e^{a \cos^2 \beta} \sin \beta\,d\beta}
+
+        P_2 = 4\pi \int_0^{\pi/2} \frac{3 \cos^2 \beta - 1)}{2}
+                        P(a, \beta) \sin \beta\,d\beta
+
+    where $a$ is the distribution width $\sigma$ handed to the weights function.
+    There is a somewhat complicated closed form solution
+
+    .. math:
+
+        P_2 = \frac{3e^a}{2\sqrt{\pi a} E_a} - \frac{3}{4a} - \frac{1}{2}
+
+    where $E_a = \mathrm{erfi}(\sqrt{a})$ is the imaginary error function,
+    which can be coded in python as::
+
+        from numpy import pi, sqrt, exp
+        from scipy.special import erfi
+
+        def P_2(a):
+            # protect against overflow with a Taylor series at infinity
+            if a <= 700:
+                r = exp(a)/sqrt(pi*a)/erfi(sqrt(a))
+            else:
+                r = 1/((((6.5525/a + 1.875)/a + 0.75)/a + 0.5)/a + 1)
+            return 1.5*r - 0.75/a - 0.5
+
+    Given an order parameter $S = P_2(a)$, one can also solve for the
+    equivalent $a$:
+
+        from scipy.optimize import fsolve
+
+        def P_2_inv(S):
+            return fsolve(lambda x: P_2(x) - S, 1.0)[0]
+
+    References
+    ----------
+
+    [1] Hardouin, et al., 1995. SANS study of a semiflexible main chain
+    liquid crystalline polyether. *Macromolecules* 28, 5427-5433.
+    """
+    type = "maier_saupe"
+    default = dict(npts=35, width=1, nsigmas=None)
+
+    # Note: center is always zero for orientation distributions
+    def _weights(self, center, sigma, lb, ub):
+        # use the width parameter as the value for Maier-Saupe "a"
+        # and find the equivalent width sigma
+        a = sigma
+        sigma = 1./sqrt(2.*a)
+
+        # Limit width to +/- 90 degrees
+        width = min(self.nsigmas*sigma, pi/2)
+        x = np.linspace(-width, width, self.npts)
+
+        # Truncate the distribution in case the parameter value is limited
+        x[(x >= radians(lb)) & (x <= radians(ub))]
+
+        # Return orientation in degrees with Maier-Saupe weights
+        # Note: weights are normalized to sum to 1, so we can scale
+        # by an arbitrary scale factor c = exp(m) to get:
+        #     w = exp(m*cos(x)**2)/c = exp(-m*sin(x)**2)
+        return degrees(x), exp(-a*sin(x)**2)
+
+
+
+def P_2(a):
+    from numpy import pi, sqrt, exp
+    from scipy.special import erfi
+
+    # Double precision e^x overflows so do a Taylor expansion around infinity
+    # for large values of a.  With five terms it is accurate to 12 digits
+    # at a = 700, and 7 digits at a = 75.
+    if a <= 700:
+        r = exp(a)/sqrt(pi*a)/erfi(sqrt(a))
+    else:
+        r = 1/((((6.5525/a + 1.875)/a + 0.75)/a + 0.5)/a + 1)
+    return 1.5*r - 0.75/a - 0.5
+
+def P_2_inv(S):
+    from scipy.optimize import fsolve
+    return fsolve(lambda x: P_2(x) - S, 1.0)[0]
+
+def P_2_numerical(a):
+    from scipy.integrate import romberg
+    from numpy import cos, sin, pi, exp
+    def num(beta):
+        return (1.5 * cos(beta)**2 - 0.5) * exp(a * cos(beta)**2) * sin(beta)
+    def denom(beta):
+        return exp(a * cos(beta)**2) * sin(beta)
+    return romberg(num, 0, pi/2) / romberg(denom, 0, pi/2)
+
+if __name__ == "__main__":
+    import sys
+    a = float(sys.argv[1])
+    sigma = 1/(2*radians(a)**2)
+    #print("P_2", P_2(a), "difference from integral", P_2(a) - P_2_numerical(a))
+    print("a=%g, sigma=%g, P_2=%g, P_2_inv(P_2(a))-a=%g"
+          % (a, sigma, P_2(a), P_2_inv(P_2(a))-a))
Index: example/weights/maier_saupe_eq.py
===================================================================
--- example/weights/maier_saupe_eq.py	(revision 35d230086190fc31926c237440d996e7ce2c269b)
+++ example/weights/maier_saupe_eq.py	(revision 35d230086190fc31926c237440d996e7ce2c269b)
@@ -0,0 +1,63 @@
+import numpy as np
+from numpy import exp, sin, degrees, radians, pi, sqrt
+
+from sasmodels.weights import Dispersion as BaseDispersion
+
+class Dispersion(BaseDispersion):
+    r"""
+    Maier-Saupe dispersion on orientation (equal weights).
+
+    .. math:
+
+        w(\theta) = e^{a{\cos^2 \theta}}
+
+    This provides a close match to the gaussian distribution for
+    low angles, but the tails are limited to $\pm 90^\circ$.  For $a \ll 1$
+    the distribution is approximately uniform.  The usual polar coordinate
+    projection applies, with $\theta$ weights scaled by $\cos \theta$
+    and $\phi$ weights unscaled.
+
+    This is equivalent to a cyclic gaussian distribution
+    $w(\theta) = e^{-sin^2(\theta)/(2\sigma^2)}.
+
+    The $\theta$ points are spaced such that each interval has an
+    equal contribution to the distribution.
+
+    This works surprisingly poorly.  Try::
+
+        $ sascomp cylinder -2d theta=45 phi=20 phi_pd_type=maier_saupe_eq \
+          phi_pd_n=100,1000 radius=50 length=2*radius -midq phi_pd=5
+
+    Leaving it here for others to improve.
+    """
+    type = "maier_saupe_eq"
+    default = dict(npts=35, width=1, nsigmas=None)
+
+    # Note: center is always zero for orientation distributions
+    def _weights(self, center, sigma, lb, ub):
+        # use the width parameter as the value for Maier-Saupe "a"
+        a = sigma
+        sigma = 1./sqrt(2.*a)
+
+        # Create a lookup table for finding n points equally spaced
+        # in the cumulative density function.
+        # Limit width to +/-90 degrees.
+        width = min(self.nsigmas*sigma, pi/2)
+        xp = np.linspace(-width, width, max(self.npts*10, 100))
+
+        # Compute CDF. Since we normalized the sum of the weights to 1,
+        # we can scale by an arbitrary scale factor c = exp(m) to get:
+        #     w = exp(m*cos(x)**2)/c = exp(-m*sin(x)**2)
+        yp = np.cumsum(exp(-a*sin(xp)**2))
+        yp /= yp[-1]
+
+        # Find the mid-points of the equal-weighted intervals in the CDF
+        y = np.linspace(0, 1, self.npts+2)[1:-1]
+        x = np.interp(y, yp, xp)
+        wx = np.ones(self.npts)
+
+        # Truncate the distribution in case the parameter value is limited
+        index = (x >= lb) & (x <= ub)
+        x, wx = x[index], wx[index]
+
+        return degrees(x), wx
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision b297ba9e837aada12de87cd46ea7a01bfe549da7)
+++ sasmodels/compare.py	(revision b297ba9e837aada12de87cd46ea7a01bfe549da7)
@@ -39,4 +39,5 @@
 
 from . import core
+from . import weights
 from . import kerneldll
 from . import kernelcl
@@ -45,5 +46,4 @@
 from .direct_model import DirectModel, get_mesh
 from .generate import FLOAT_RE, set_integration_size
-from .weights import plot_weights
 
 # pylint: disable=unused-import
@@ -115,10 +115,11 @@
 
     === environment variables ===
-    -DSAS_MODELPATH=path sets directory containing custom models
+    -DSAS_MODELPATH=~/.sasmodels/custom_models sets path to custom models
+    -DSAS_WEIGHTS_PATH=~/.sasview/weights sets path to custom distributions
     -DSAS_OPENCL=vendor:device|cuda:device|none sets the target GPU 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
+    -DSAS_OPENMP=0 set to 1 to turn on OpenMP for the DLLs
+    -DSAS_DLL_PATH=~/.sasmodels/compiled_models sets the DLL cache
 
 The interpretation of quad precision depends on architecture, and may
@@ -784,5 +785,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
@@ -1441,4 +1442,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/models/__init__.py
===================================================================
--- sasmodels/models/__init__.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/__init__.py	(revision d827c5e8ff3afdce75ac38dbce632ac8eb18a386)
@@ -1,3 +1,98 @@
 """
-1D Modeling for SAS
+Model definition files
+----------------------
+
+The models below are grouped by type.  The list is a snapshot at a particular
+time and may be out of date.
+
+Models with pure form factor (all of which define *F(Q)*):
+
+    :mod:`barbell`
+    :mod:`capped_cylinder`
+    :mod:`core_multi_shell`
+    :mod:`core_shell_bicelle`
+    :mod:`core_shell_bicelle_elliptical`
+    :mod:`core_shell_bicelle_elliptical_belt_rough`
+    :mod:`core_shell_cylinder`
+    :mod:`core_shell_ellipsoid`
+    :mod:`core_shell_parallelepipied`
+    :mod:`core_shell_sphere`
+    :mod:`cylinder` [limiting conditions (long rods, thin disks) not implemented]
+    :mod:`ellipsoid`
+    :mod:`elliptical_cylinder`
+    :mod:`fuzzy_sphere`
+    :mod:`hollow_cylinder`
+    :mod:`hollow_rectangular_prism`
+    :mod:`hollow_rectangular_prism_thin_walls`
+    :mod:`multilayer_vesicle`
+    :mod:`onion`
+    :mod:`parallelepiped`
+    :mod:`rectangular_prism`
+    :mod:`sphere`
+    :mod:`spherical_sld`
+    :mod:`triaxial_ellipsoid`
+    :mod:`vesicle`
+
+Models with local structure factor:
+
+    :mod:`flexible_cylinder`
+    :mod:`flexible_cylinder_elliptical`
+    :mod:`linear_pearls`
+    :mod:`mono_gauss_coil`
+    :mod:`pearl_necklace`
+    :mod:`poly_gauss_coil`
+    :mod:`polymer_micelle`
+    :mod:`pringle`
+    :mod:`raspberry`
+    :mod:`stacked_disks`
+    :mod:`star_polymer`
+
+Models with long range structure factor:
+
+    :mod:`binary_hard_sphere`
+    :mod:`bcc_paracrystal`
+    :mod:`fcc_paracrystal`
+    :mod:`fractal`
+    :mod:`fractal_core_shell`
+    :mod:`lamellar`
+    :mod:`lamellar_hg`
+    :mod:`lamellar_hg_stack_caille`
+    :mod:`lamellar_stack_caille`
+    :mod:`lamellar_stack_paracrystal`
+    :mod:`mass_fractal`
+    :mod:`mass_surface_fractal`
+    :mod:`rpa`
+    :mod:`sc_paracrystal`
+    :mod:`surface_fractal`
+
+Models which are pure structure factors::
+
+    :mod:`hardsphere`
+    :mod:`hayter_msa`
+    :mod:`squarewell`
+    :mod:`stickyhardsphere`
+
+Other models:
+
+    :mod:`adsorbed_layer`
+    :mod:`be_polyelectrolyte`
+    :mod:`broad_peak`
+    :mod:`correlation_length`
+    :mod:`dab`
+    :mod:`gauss_lorentz_gel`
+    :mod:`gaussian_peak`
+    :mod:`gel_fit`
+    :mod:`guinier_porod`
+    :mod:`guinier`
+    :mod:`line`
+    :mod:`lorentz`
+    :mod:`peak_lorentz`
+    :mod:`polymer_excl_volume`
+    :mod:`porod`
+    :mod:`power_law`
+    :mod:`spinodal`
+    :mod:`teubner_strey`
+    :mod:`two_lorentzian`
+    :mod:`unified_power_Rg`
+
 """
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision a34b8119b864fa10b972c2a11f92a9c5e13207f1)
+++ sasmodels/sasview_model.py	(revision a34b8119b864fa10b972c2a11f92a9c5e13207f1)
@@ -31,4 +31,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 b297ba9e837aada12de87cd46ea7a01bfe549da7)
+++ sasmodels/weights.py	(revision b297ba9e837aada12de87cd46ea7a01bfe549da7)
@@ -230,4 +230,26 @@
 ))
 
+SAS_WEIGHTS_PATH = "~/.sasview/weights"
+def load_weights(pattern=None):
+    # type: (str) -> None
+    """
+    Load dispersion distributions matching the given glob pattern
+    """
+    import logging
+    import os
+    import os.path
+    import glob
+    import traceback
+    from .custom import load_custom_kernel_module
+    if pattern is None:
+        path = os.environ.get("SAS_WEIGHTS_PATH", SAS_WEIGHTS_PATH)
+        pattern = os.path.join(path, "*.py")
+    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):
