source: sasmodels/sasmodels/models/hardsphere.py @ c1799d3

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since c1799d3 was c1799d3, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

Merge branch 'beta_approx' into ticket-1157

  • Property mode set to 100644
File size: 6.7 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
41J K Percus, J Yevick, *J. Phys. Rev.*, 110, (1958) 1
42"""
43
44import numpy as np
45from numpy import inf
46
47name = "hardsphere"
48title = "Hard sphere structure factor, with Percus-Yevick closure"
49description = """\
50    [Hard sphere structure factor, with Percus-Yevick closure]
51        Interparticle S(Q) for random, non-interacting spheres.
52    May be a reasonable approximation for other shapes of
53    particles that freely rotate, and for moderately polydisperse
54        systems. Though strictly the maths needs to be modified -
55    which sasview does not do yet.
56    radius_effective is the hard sphere radius
57    volfraction is the volume fraction occupied by the spheres.
58"""
59category = "structure-factor"
60structure_factor = True
61single = False # TODO: check
62
63#             ["name", "units", default, [lower, upper], "type","description"],
64parameters = [["radius_effective", "Ang", 50.0, [0, inf], "",
65               "effective radius of hard sphere"],
66              ["volfraction", "", 0.2, [0, 0.74], "",
67               "volume fraction of hard spheres"],
68             ]
69
70Iq = r"""
71      double D,A,B,G,X,X2,X4,S,C,FF,HARDSPH;
72      // these are c compiler instructions, can also put normal code inside the "if else" structure
73      #if FLOAT_SIZE > 4
74      // double precision
75      // orig had 0.2, don't call the variable cutoff as PAK already has one called that!
76      // Must use UPPERCASE name please.
77      // 0.05 better, 0.1 OK
78      #define CUTOFFHS 0.05
79      #else
80      // 0.1 bad, 0.2 OK, 0.3 good, 0.4 better, 0.8 no good
81      #define CUTOFFHS 0.4
82      #endif
83
84      if(fabs(radius_effective) < 1.E-12) {
85               HARDSPH=1.0;
86//printf("HS1 %g: %g\n",q,HARDSPH);
87               return(HARDSPH);
88      }
89      // removing use of pow(xxx,2) and rearranging the calcs
90      // of A, B & G cut ~40% off execution time ( 0.5 to 0.3 msec)
91      X = 1.0/( 1.0 -volfraction);
92      D= X*X;
93      A= (1.+2.*volfraction)*D;
94      A *=A;
95      X=fabs(q*radius_effective*2.0);
96
97      if(X < 5.E-06) {
98                 HARDSPH=1./A;
99//printf("HS2 %g: %g\n",q,HARDSPH);
100                 return(HARDSPH);
101      }
102      X2 =X*X;
103      B = (1.0 +0.5*volfraction)*D;
104      B *= B;
105      B *= -6.*volfraction;
106      G=0.5*volfraction*A;
107
108      if(X < CUTOFFHS) {
109      // RKH Feb 2016, use Taylor series expansion for small X
110      // else no obvious way to rearrange the equations to avoid
111      // needing a very high number of significant figures.
112      // Series expansion found using Mathematica software. Numerical test
113      // in .xls showed terms to X^2 are sufficient
114      // for 5 or 6 significant figures, but I put the X^4 one in anyway
115            //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;
116            // refactoring the polynomial makes it very slightly faster (0.5 not 0.6 msec)
117            //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;
118
119            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;
120
121            // combining the terms makes things worse at smallest Q in single precision
122            //FF = (8-0.8*X2)*A +(3.0-X2/3.)*2*B + (4+0.5*X2)*G +(A/35. +B/40. +G/50.)*X4;
123            // note that G = -volfraction*A/2, combining this makes no further difference at smallest Q
124            //FF = (8 +2.*volfraction + ( volfraction/4. -0.8 +(volfraction/100. -1./35.)*X2 )*X2 )*A  + (3.0 -X2/3. +X4/40.)*2.*B;
125            HARDSPH= 1./(1. + volfraction*FF );
126//printf("HS3 %g: %g\n",q,HARDSPH);
127            return(HARDSPH);
128      }
129      X4=X2*X2;
130      SINCOS(X,S,C);
131
132// RKH Feb 2016, use version FISH code as is better than original sasview one
133// at small Q in single precision, and more than twice as fast in double.
134      //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;
135      // refactoring the polynomial here & above makes it slightly faster
136
137      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 ;
138      HARDSPH= 1./(1. + 24.*volfraction*FF/X2 );
139
140      // changing /X and /X2 to *MX1 and *MX2, no significantg difference?
141      //MX=1.0/X;
142      //MX2=MX*MX;
143      //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)) ;
144      //HARDSPH= 1./(1. + 24.*volfraction*FF*MX2*MX );
145
146// grouping the terms, was about same as sasmodels for single precision issues
147//     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 );
148//     HARDSPH= 1./(1. + 24.*volfraction*FF/X2 );
149// remove 1/X2 from final line, take more powers of X inside the brackets, stil bad
150//      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;
151//      HARDSPH= 1./(1. + 24.*volfraction*FF );
152//printf("HS4 %g: %g\n",q,HARDSPH);
153      return(HARDSPH);
154   """
155
156def random():
157    pars = dict(
158        scale=1, background=0,
159        radius_effective=10**np.random.uniform(1, 4),
160        volfraction=10**np.random.uniform(-2, 0),  # high volume fraction
161    )
162    return pars
163
164# Q=0.001 is in the Taylor series, low Q part, so add Q=0.1,
165# assuming double precision sasview is correct
166tests = [
167    [{'scale': 1.0, 'background' : 0.0, 'radius_effective' : 50.0,
168      'volfraction' : 0.2, 'radius_effective_pd' : 0},
169     [0.001, 0.1], [0.209128, 0.930587]],
170]
171# ADDED by: RKH  ON: 16Mar2016  using equations from FISH as better than
172# orig sasview, see notes above. Added Taylor expansions at small Q.
Note: See TracBrowser for help on using the repository browser.