Changeset a5a12ca in sasmodels


Ignore:
Timestamp:
Oct 28, 2017 2:02:56 AM (2 years ago)
Author:
dirk
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
75e4319
Parents:
6db17bd
Message:

implementation of Boltzmann distribution refs #1018

Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed
  • doc/guide/pd/polydispersity.rst

    r1f058ea ra5a12ca  
    4040calculations are generally more robust with more data points or more angles. 
    4141 
    42 The following five distribution functions are provided: 
     42The following six distribution functions are provided: 
    4343 
    4444*  *Rectangular Distribution* 
     
    4747*  *Schulz Distribution* 
    4848*  *Array Distribution* 
     49*  *Boltzmann Distribution* 
    4950 
    5051These are all implemented as *number-average* distributions. 
     
    181182^^^^^^^^^^^^^^^^^^ 
    182183 
    183 This user-definable distribution should be given as as a simple ASCII text 
     184This user-definable distribution should be given as a simple ASCII text 
    184185file where the array is defined by two columns of numbers: $x$ and $f(x)$. 
    185186The $f(x)$ will be normalized to 1 during the computation. 
     
    200201given for the model will have no affect, and will be ignored when computing 
    201202the average.  This means that any parameter with an array distribution will 
    202 not be fittable. 
     203not be fitable. 
     204 
     205.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     206 
     207Boltzmann Distribution 
     208^^^^^^^^^^^^^^^^^^^^^^ 
     209 
     210The Boltzmann Distribution is defined as 
     211 
     212.. math:: 
     213 
     214    f(x) = \frac{1}{\text{Norm}} 
     215           \exp\left(-\frac{ | x - \bar x | }{\sigma}\right) 
     216 
     217where $\bar x$ is the mean of the distribution and *Norm* is a normalization 
     218factor which is determined during the numerical calculation. 
     219The width is defined as 
     220 
     221.. math:: \sigma=\frac{k T}{E} 
     222 
     223which is the inverse Boltzmann factor, 
     224where $k$ is the Boltzmann constant, $T$ the temperature in Kelvin and $E$ a 
     225characteristic energy per particle. 
     226 
     227.. figure:: pd_boltzmann.jpg 
     228 
     229    Boltzmann distribution. 
    203230 
    204231.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
  • sasmodels/weights.py

    r41e7f2e ra5a12ca  
    5555        """ 
    5656        sigma = self.width * center if relative else self.width 
     57        if not relative: 
     58            # For orientation, the jitter is relative to 0 not the angle 
     59            center = 0 
     60            pass 
    5761        if sigma == 0 or self.npts < 2: 
    5862            if lb <= center <= ub: 
     
    186190        return x, px 
    187191 
     192class BoltzmannDispersion(Dispersion): 
     193    r""" 
     194    Boltzmann dispersion, with $\sigma=k T/E$. 
     195 
     196    .. math:: 
     197 
     198        w = \exp\left( -|x - c|/\sigma\right) 
     199    """ 
     200    type = "boltzmann" 
     201    default = dict(npts=35, width=0, nsigmas=3) 
     202    def _weights(self, center, sigma, lb, ub): 
     203        x = self._linspace(center, sigma, lb, ub) 
     204        px = np.exp(-np.abs(x-center) / np.abs(sigma)) 
     205        return x, px 
    188206 
    189207# dispersion name -> disperser lookup table. 
     
    196214    GaussianDispersion, 
    197215    SchulzDispersion, 
     216    BoltzmannDispersion 
    198217)) 
    199218 
     
    225244    obj = cls(n, width, nsigmas) 
    226245    v, w = obj.get_weights(value, limits[0], limits[1], relative) 
    227     return v, w 
    228  
    229  
    230 def plot_weights(model_info, pairs): 
    231     # type: (ModelInfo, List[Tuple[np.ndarray, np.ndarray]]) -> None 
     246    return v, w/np.sum(w) 
     247 
     248 
     249def plot_weights(model_info, mesh): 
     250    # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None 
    232251    """ 
    233252    Plot the weights returned by :func:`get_weights`. 
    234253 
    235     *model_info* is 
    236     :param model_info: 
    237     :param pairs: 
    238     :return: 
     254    *model_info* defines model parameters, etc. 
     255 
     256    *mesh* is a list of tuples containing (*value*, *dispersity*, *weights*) 
     257    for each parameter, where (*dispersity*, *weights*) pairs are the 
     258    distributions to be plotted. 
    239259    """ 
    240260    import pylab 
    241261 
    242     if any(len(values)>1 for values, weights in pairs): 
     262    if any(len(dispersity)>1 for value, dispersity, weights in mesh): 
    243263        labels = [p.name for p in model_info.parameters.call_parameters] 
    244         pylab.interactive(True) 
     264        #pylab.interactive(True) 
    245265        pylab.figure() 
    246         for (v,w), s in zip(pairs, labels): 
    247             if len(v) > 1: 
    248                 #print("weights for", s, v, w) 
    249                 pylab.plot(v, w, '-o', label=s) 
     266        for (v,x,w), s in zip(mesh, labels): 
     267            if len(x) > 1: 
     268                pylab.plot(x, w, '-o', label=s) 
    250269        pylab.grid(True) 
    251270        pylab.legend() 
Note: See TracChangeset for help on using the changeset viewer.