source: sasmodels/sasmodels/models/hardsphere.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 years ago

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

  • Property mode set to 100644
File size: 7.0 KB
Line 
1# Note: model title and parameter table are inserted automatically
2r"""Calculate the interparticle structure factor for monodisperse
3spherical particles interacting through hard sphere (excluded volume)
4interactions.
5May be a reasonable approximation for other shapes of particles that
6freely rotate, and for moderately polydisperse systems. Though strictly
7the maths needs to be modified (no \Beta(Q) correction yet in sasview).
8
9radius_effective is the effective hard sphere radius.
10volfraction is the volume fraction occupied by the spheres.
11
12In sasview the effective radius may be calculated from the parameters
13used in the form factor $P(q)$ that this $S(q)$ is combined with.
14
15For numerical stability the computation uses a Taylor series expansion
16at very small $qR$, there may be a very minor glitch at the transition point
17in some circumstances.
18
19The S(Q) uses the Percus-Yevick closure where the interparticle
20potential is
21
22.. math::
23
24    U(r) = \begin{cases}
25    \infty & r < 2R \\
26    0 & r \geq 2R
27    \end{cases}
28
29where $r$ is the distance from the center of the sphere of a radius $R$.
30
31For a 2D plot, the wave transfer is defined as
32
33.. math::
34
35    q = \sqrt{q_x^2 + q_y^2}
36
37
38References
39----------
40
41.. [#] J K Percus, J Yevick, *J. Phys. Rev.*, 110, (1958) 1
42
43Source
44------
45
46`hardsphere.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/hardsphere.py>`_
47
48Authorship and Verification
49----------------------------
50
51* **Author:**
52* **Last Modified by:**
53* **Last Reviewed by:**
54* **Source added by :** Steve King **Date:** March 25, 2019
55"""
56
57import numpy as np
58from numpy import inf
59
60name = "hardsphere"
61title = "Hard sphere structure factor, with Percus-Yevick closure"
62description = """\
63    [Hard sphere structure factor, with Percus-Yevick closure]
64        Interparticle S(Q) for random, non-interacting spheres.
65    May be a reasonable approximation for other shapes of
66    particles that freely rotate, and for moderately polydisperse
67        systems. Though strictly the maths needs to be modified -
68    which sasview does not do yet.
69    radius_effective is the hard sphere radius
70    volfraction is the volume fraction occupied by the spheres.
71"""
72category = "structure-factor"
73structure_factor = True
74single = False # TODO: check
75
76#             ["name", "units", default, [lower, upper], "type","description"],
77parameters = [["radius_effective", "Ang", 50.0, [0, inf], "",
78               "effective radius of hard sphere"],
79              ["volfraction", "", 0.2, [0, 0.74], "",
80               "volume fraction of hard spheres"],
81             ]
82
83Iq = r"""
84      double D,A,B,G,X,X2,X4,S,C,FF,HARDSPH;
85      // these are c compiler instructions, can also put normal code inside the "if else" structure
86      #if FLOAT_SIZE > 4
87      // double precision
88      // orig had 0.2, don't call the variable cutoff as PAK already has one called that!
89      // Must use UPPERCASE name please.
90      // 0.05 better, 0.1 OK
91      #define CUTOFFHS 0.05
92      #else
93      // 0.1 bad, 0.2 OK, 0.3 good, 0.4 better, 0.8 no good
94      #define CUTOFFHS 0.4
95      #endif
96
97      if(fabs(radius_effective) < 1.E-12) {
98               HARDSPH=1.0;
99//printf("HS1 %g: %g\n",q,HARDSPH);
100               return(HARDSPH);
101      }
102      // removing use of pow(xxx,2) and rearranging the calcs
103      // of A, B & G cut ~40% off execution time ( 0.5 to 0.3 msec)
104      X = 1.0/( 1.0 -volfraction);
105      D= X*X;
106      A= (1.+2.*volfraction)*D;
107      A *=A;
108      X=fabs(q*radius_effective*2.0);
109
110      if(X < 5.E-06) {
111                 HARDSPH=1./A;
112//printf("HS2 %g: %g\n",q,HARDSPH);
113                 return(HARDSPH);
114      }
115      X2 =X*X;
116      B = (1.0 +0.5*volfraction)*D;
117      B *= B;
118      B *= -6.*volfraction;
119      G=0.5*volfraction*A;
120
121      if(X < CUTOFFHS) {
122      // RKH Feb 2016, use Taylor series expansion for small X
123      // else no obvious way to rearrange the equations to avoid
124      // needing a very high number of significant figures.
125      // Series expansion found using Mathematica software. Numerical test
126      // in .xls showed terms to X^2 are sufficient
127      // for 5 or 6 significant figures, but I put the X^4 one in anyway
128            //FF = 8*A +6*B + 4*G - (0.8*A +2.0*B/3.0 +0.5*G)*X2 +(A/35. +B/40. +G/50.)*X4;
129            // refactoring the polynomial makes it very slightly faster (0.5 not 0.6 msec)
130            //FF = 8*A +6*B + 4*G + ( -0.8*A -2.0*B/3.0 -0.5*G +(A/35. +B/40. +G/50.)*X2)*X2;
131
132            FF = 8.0*A +6.0*B + 4.0*G + ( -0.8*A -B/1.5 -0.5*G +(A/35. +0.0125*B +0.02*G)*X2)*X2;
133
134            // combining the terms makes things worse at smallest Q in single precision
135            //FF = (8-0.8*X2)*A +(3.0-X2/3.)*2*B + (4+0.5*X2)*G +(A/35. +B/40. +G/50.)*X4;
136            // note that G = -volfraction*A/2, combining this makes no further difference at smallest Q
137            //FF = (8 +2.*volfraction + ( volfraction/4. -0.8 +(volfraction/100. -1./35.)*X2 )*X2 )*A  + (3.0 -X2/3. +X4/40.)*2.*B;
138            HARDSPH= 1./(1. + volfraction*FF );
139//printf("HS3 %g: %g\n",q,HARDSPH);
140            return(HARDSPH);
141      }
142      X4=X2*X2;
143      SINCOS(X,S,C);
144
145// RKH Feb 2016, use version FISH code as is better than original sasview one
146// at small Q in single precision, and more than twice as fast in double.
147      //FF=A*(S-X*C)/X + B*(2.*X*S -(X2-2.)*C -2.)/X2 + G*( (4.*X2*X -24.*X)*S -(X4 -12.*X2 +24.)*C +24. )/X4;
148      // refactoring the polynomial here & above makes it slightly faster
149
150      FF=  (( G*( (4.*X2 -24.)*X*S -(X4 -12.*X2 +24.)*C +24. )/X2 + B*(2.*X*S -(X2-2.)*C -2.) )/X + A*(S-X*C))/X ;
151      HARDSPH= 1./(1. + 24.*volfraction*FF/X2 );
152
153      // changing /X and /X2 to *MX1 and *MX2, no significantg difference?
154      //MX=1.0/X;
155      //MX2=MX*MX;
156      //FF=  (( G*( (4.*X2 -24.)*X*S -(X4 -12.*X2 +24.)*C +24. )*MX2 + B*(2.*X*S -(X2-2.)*C -2.) )*MX + A*(S-X*C)) ;
157      //HARDSPH= 1./(1. + 24.*volfraction*FF*MX2*MX );
158
159// grouping the terms, was about same as sasmodels for single precision issues
160//     FF=A*(S/X-C) + B*(2.*S/X - C +2.0*(C-1.0)/X2) + G*( (4./X -24./X3)*S -(1.0 -12./X2 +24./X4)*C +24./X4 );
161//     HARDSPH= 1./(1. + 24.*volfraction*FF/X2 );
162// remove 1/X2 from final line, take more powers of X inside the brackets, stil bad
163//      FF=A*(S/X3-C/X2) + B*(2.*S/X3 - C/X2 +2.0*(C-1.0)/X4) + G*( (4./X -24./X3)*S -(1.0 -12./X2 +24./X4)*C +24./X4 )/X2;
164//      HARDSPH= 1./(1. + 24.*volfraction*FF );
165//printf("HS4 %g: %g\n",q,HARDSPH);
166      return(HARDSPH);
167   """
168
169def random():
170    """Return a random parameter set for the model."""
171    pars = dict(
172        scale=1, background=0,
173        radius_effective=10**np.random.uniform(1, 4),
174        volfraction=10**np.random.uniform(-2, 0),  # high volume fraction
175    )
176    return pars
177
178# Q=0.001 is in the Taylor series, low Q part, so add Q=0.1,
179# assuming double precision sasview is correct
180tests = [
181    [{'scale': 1.0, 'background' : 0.0, 'radius_effective' : 50.0,
182      'volfraction' : 0.2, 'radius_effective_pd' : 0},
183     [0.001, 0.1], [0.209128, 0.930587]],
184]
185# ADDED by: RKH  ON: 16Mar2016  using equations from FISH as better than
186# orig sasview, see notes above. Added Taylor expansions at small Q.
Note: See TracBrowser for help on using the repository browser.