Changeset 35d2300 in sasmodels


Ignore:
Timestamp:
Sep 13, 2018 8:08:42 PM (2 months ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
ticket-608-user-defined-weights
Children:
b50e28b
Parents:
a5cb9bc
Message:

Relate order parameter P_2 to 'a' value in Maier-Saupe dist

Location:
example/weights
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • example/weights/cyclic_gaussian.py

    r55e82f0 r35d2300  
    1919 
    2020    This is eqivalent to a Maier-Saupe distribution with order 
    21     parameter $P_2 = 1/(2 \sigma^2)$, with $\sigma$ in radians. 
     21    parameter $a = 1/(2 \sigma^2)$, with $\sigma$ in radians. 
    2222    """ 
    2323    type = "cyclic_gaussian" 
  • example/weights/maier_saupe.py

    r55e82f0 r35d2300  
    1111    .. math: 
    1212 
    13         w(\theta) = e^{P_2{\cos^2 \theta}} 
     13        w(\theta) = e^{a{\cos^2 \theta}} 
    1414 
    1515    This provides a close match to the gaussian distribution for 
    16     low angles, but the tails are limited to $\pm 90^\circ$.  For $P_2 \ll 1$ 
     16    low angles, but the tails are limited to $\pm 90^\circ$.  For $a \ll 1$ 
    1717    the distribution is approximately uniform.  The usual polar coordinate 
    1818    projection applies, with $\theta$ weights scaled by $\cos \theta$ 
     
    2020 
    2121    This is equivalent to a cyclic gaussian distribution 
    22     $w(\theta) = e^{-sin^2(\theta)/(2\P_2^2)}. 
     22    $w(\theta) = e^{-sin^2(\theta)/(2\a^2)}. 
     23 
     24    The order parameter $P_2$ is defined as 
     25 
     26    .. math: 
     27 
     28        P(a, \beta) = \frac{e^{a \cos^2 \beta}}{ 
     29                       4\pi \int_0^{\pi/2} e^{a \cos^2 \beta} \sin \beta\,d\beta} 
     30 
     31        P_2 = 4\pi \int_0^{\pi/2} \frac{3 \cos^2 \beta - 1)}{2} 
     32                        P(a, \beta) \sin \beta\,d\beta 
     33 
     34    where $a$ is the distribution width $\sigma$ handed to the weights function. 
     35    There is a somewhat complicated closed form solution 
     36 
     37    .. math: 
     38 
     39        P_2 = \frac{3e^a}{2\sqrt{\pi a} E_a} - \frac{3}{4a} - \frac{1}{2} 
     40 
     41    where $E_a = \mathrm{erfi}(\sqrt{a})$ is the imaginary error function, 
     42    which can be coded in python as:: 
     43 
     44        from numpy import pi, sqrt, exp 
     45        from scipy.special import erfi 
     46 
     47        def P_2(a): 
     48            # protect against overflow with a Taylor series at infinity 
     49            if a <= 700: 
     50                r = exp(a)/sqrt(pi*a)/erfi(sqrt(a)) 
     51            else: 
     52                r = 1/((((6.5525/a + 1.875)/a + 0.75)/a + 0.5)/a + 1) 
     53            return 1.5*r - 0.75/a - 0.5 
     54 
     55    References 
     56    ---------- 
     57 
     58    [1] Hardouin, et al., 1995. SANS study of a semiflexible main chain 
     59    liquid crystalline polyether. *Macromolecules* 28, 5427–5433. 
    2360    """ 
    2461    type = "maier_saupe" 
     
    2764    # Note: center is always zero for orientation distributions 
    2865    def _weights(self, center, sigma, lb, ub): 
    29         # use the width parameter as the value for Maier-Saupe P_2 
     66        # use the width parameter as the value for Maier-Saupe "a" 
    3067        # and find the equivalent width sigma 
    31         P2 = sigma 
    32         sigma = 1./sqrt(2.*P2) 
     68        a = sigma 
     69        sigma = 1./sqrt(2.*a) 
    3370 
    3471        # Limit width to +/- 90 degrees 
     
    4380        # by an arbitrary scale factor c = exp(m) to get: 
    4481        #     w = exp(m*cos(x)**2)/c = exp(-m*sin(x)**2) 
    45         return degrees(x), exp(-P2*sin(x)**2) 
     82        return degrees(x), exp(-a*sin(x)**2) 
     83 
     84 
     85 
     86def P_2(a): 
     87    from numpy import pi, sqrt, exp 
     88    from scipy.special import erfi 
     89 
     90    # Double precision e^x overflows so do a Taylor expansion around infinity 
     91    # for large values of a.  With five terms it is accurate to 12 digits 
     92    # at a = 700, and 7 digits at a = 75. 
     93    if a <= 700: 
     94        r = exp(a)/sqrt(pi*a)/erfi(sqrt(a)) 
     95    else: 
     96        r = 1/((((6.5525/a + 1.875)/a + 0.75)/a + 0.5)/a + 1) 
     97    return 1.5*r - 0.75/a - 0.5 
     98 
     99def P_2_numerical(a): 
     100    from scipy.integrate import romberg 
     101    from numpy import cos, sin, pi, exp 
     102    def num(beta): 
     103        return (1.5 * cos(beta)**2 - 0.5) * exp(a * cos(beta)**2) * sin(beta) 
     104    def denom(beta): 
     105        return exp(a * cos(beta)**2) * sin(beta) 
     106    return romberg(num, 0, pi/2) / romberg(denom, 0, pi/2) 
     107 
     108if __name__ == "__main__": 
     109    import sys 
     110    a = float(sys.argv[1]) 
     111    #print("P_2", P_2(a), "difference from integral", P_2(a) - P_2_numerical(a)) 
     112    print("P_2", P_2(a)) 
  • example/weights/maier_saupe_eq.py

    r55e82f0 r35d2300  
    1010    .. math: 
    1111 
    12         w(\theta) = e^{P_2{\cos^2 \theta}} 
     12        w(\theta) = e^{a{\cos^2 \theta}} 
    1313 
    1414    This provides a close match to the gaussian distribution for 
    15     low angles, but the tails are limited to $\pm 90^\circ$.  For $P_2 \ll 1$ 
     15    low angles, but the tails are limited to $\pm 90^\circ$.  For $a \ll 1$ 
    1616    the distribution is approximately uniform.  The usual polar coordinate 
    1717    projection applies, with $\theta$ weights scaled by $\cos \theta$ 
     
    3636    # Note: center is always zero for orientation distributions 
    3737    def _weights(self, center, sigma, lb, ub): 
    38         # use the width parameter as the value for Maier-Saupe P_2 
    39         P2 = sigma 
    40         sigma = 1./sqrt(2.*P2) 
     38        # use the width parameter as the value for Maier-Saupe "a" 
     39        a = sigma 
     40        sigma = 1./sqrt(2.*a) 
    4141 
    4242        # Create a lookup table for finding n points equally spaced 
     
    4949        # we can scale by an arbitrary scale factor c = exp(m) to get: 
    5050        #     w = exp(m*cos(x)**2)/c = exp(-m*sin(x)**2) 
    51         yp = np.cumsum(exp(-P2*sin(xp)**2)) 
     51        yp = np.cumsum(exp(-a*sin(xp)**2)) 
    5252        yp /= yp[-1] 
    5353 
Note: See TracChangeset for help on using the changeset viewer.