source: sasmodels/sasmodels/models/squarewell.py @ 0507e09

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0507e09 was 0507e09, checked in by smk78, 5 months ago

Added link to source code to each model. Closes #883

  • Property mode set to 100644
File size: 5.0 KB
Line 
1# Note: model title and parameter table are inserted automatically
2r"""
3This calculates the interparticle structure factor for a square well fluid
4spherical particles. The mean spherical approximation (MSA) closure was
5used for this calculation, and is not the most appropriate closure for
6an attractive interparticle potential. This solution has been compared
7to Monte Carlo simulations for a square well fluid, showing this calculation
8to be limited in applicability to well depths $\epsilon < 1.5$ kT and
9volume fractions $\phi < 0.08$.
10
11Positive well depths correspond to an attractive potential well. Negative
12well depths correspond to a potential "shoulder", which may or may not be
13physically reasonable. The stickyhardsphere model may be a better choice in
14some circumstances. Computed values may behave badly at extremely small $qR$.
15
16The well width $(\lambda)$ is defined as multiples of the particle diameter
17$(2 R)$.
18
19The interaction potential is:
20
21  .. image:: img/squarewell.png
22
23.. math::
24
25    U(r) = \begin{cases}
26    \infty & r < 2R \\
27    -\epsilon & 2R \leq r < 2R\lambda \\
28    0 & r \geq 2R\lambda
29    \end{cases}
30
31where $r$ is the distance from the center of the sphere of a radius $R$.
32
33In sasview the effective radius may be calculated from the parameters
34used in the form factor $P(q)$ that this $S(q)$ is combined with.
35
36For 2D data: The 2D scattering intensity is calculated in the same way as 1D,
37where the $q$ vector is defined as
38
39.. math::
40
41    q = \sqrt{q_x^2 + q_y^2}
42
43References
44----------
45
46.. [#] R V Sharma, K C Sharma, *Physica*, 89A (1977) 213
47
48Source
49------
50
51`squarewell.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/squarewell.py>`_
52
53Authorship and Verification
54----------------------------
55
56* **Author:**
57* **Last Modified by:**
58* **Last Reviewed by:**
59* **Source added by :** Steve King **Date:** March 25, 2019
60"""
61
62import numpy as np
63from numpy import inf
64
65name = "squarewell"
66title = "Square well structure factor, with MSA closure"
67description = """\
68    [Square well structure factor, with MSA closure]
69        Interparticle structure factor S(Q)for a hard sphere fluid with
70        a narrow attractive well. Fits are prone to deliver non-physical
71        parameters, use with care and read the references in the full manual.
72        In sasview the effective radius will be calculated from the
73        parameters used in P(Q).
74"""
75category = "structure-factor"
76structure_factor = True
77single = False
78
79#single = False
80#             ["name", "units", default, [lower, upper], "type","description"],
81parameters = [
82    #   [ "name", "units", default, [lower, upper], "type",
83    #     "description" ],
84    ["radius_effective", "Ang", 50.0, [0, inf], "volume",
85     "effective radius of hard sphere"],
86    ["volfraction", "", 0.04, [0, 0.08], "",
87     "volume fraction of spheres"],
88    ["welldepth", "kT", 1.5, [0.0, 1.5], "",
89     "depth of well, epsilon"],
90    ["wellwidth", "diameters", 1.2, [1.0, inf], "",
91     "width of well in diameters (=2R) units, must be > 1"],
92    ]
93
94# No volume normalization despite having a volume parameter
95# This should perhaps be volume normalized?
96form_volume = """
97    return 1.0;
98    """
99
100Iq = """
101// single precision is very poor at extreme small Q, would need a Taylor series
102        double req,phis,edibkb,lambda,struc;
103        double sigma,eta,eta2,eta3,eta4,etam1,etam14,alpha,beta,gamm;
104        double x,sk,sk2,sk3,sk4,t1,t2,t3,t4,ck;
105        double S,C,SL,CL;
106        x= q;
107
108        req = radius_effective;
109        phis = volfraction;
110        edibkb = welldepth;
111        lambda = wellwidth;
112
113        sigma = req*2.;
114        eta = phis;
115        eta2 = eta*eta;
116        eta3 = eta*eta2;
117        eta4 = eta*eta3;
118        etam1 = 1. - eta;
119        etam14 = etam1*etam1*etam1*etam1;
120        // temp borrow sk for an intermediate calc
121        sk = 1.0 +2.0*eta;
122        sk *= sk;
123        alpha = (  sk + eta3*( eta-4.0 )  )/etam14;
124        beta = -(eta/3.0) * ( 18. + 20.*eta - 12.*eta2 + eta4 )/etam14;
125        gamm = 0.5*eta*( sk + eta3*(eta-4.) )/etam14;
126
127        //  calculate the structure factor
128
129        sk = x*sigma;
130        sk2 = sk*sk;
131        sk3 = sk*sk2;
132        sk4 = sk3*sk;
133        SINCOS(sk,S,C);
134        SINCOS(lambda*sk,SL,CL);
135        t1 = alpha * sk3 * ( S - sk * C );
136        t2 = beta * sk2 * 2.0*( sk*S - (0.5*sk2 - 1.)*C - 1.0 );
137        t3 = gamm*( ( 4.0*sk3 - 24.*sk ) * S - ( sk4 - 12.0*sk2 + 24.0 )*C + 24.0 );
138        t4 = -edibkb*sk3*(SL +sk*(C - lambda*CL) - S );
139        ck =  -24.0*eta*( t1 + t2 + t3 + t4 )/sk3/sk3;
140        struc  = 1./(1.-ck);
141
142        return(struc);
143"""
144
145def random():
146    """Return a random parameter set for the model."""
147    pars = dict(
148        scale=1, background=0,
149        radius_effective=10**np.random.uniform(1, 4.7),
150        volfraction=np.random.uniform(0.00001, 0.08),
151        welldepth=np.random.uniform(0, 1.5),
152        wellwidth=np.random.uniform(1, 1.2),
153    )
154    return pars
155
156demo = dict(radius_effective=50, volfraction=0.04, welldepth=1.5,
157            wellwidth=1.2, radius_effective_pd=0, radius_effective_pd_n=0)
158#
159tests = [
160    [{'scale': 1.0, 'background': 0.0, 'radius_effective': 50.0,
161      'volfraction': 0.04, 'welldepth': 1.5, 'wellwidth': 1.2,
162      'radius_effective_pd': 0}, [0.001], [0.97665742]],
163    ]
164# ADDED by: converting from sasview RKH  ON: 16Mar2016
Note: See TracBrowser for help on using the repository browser.