source: sasmodels/sasmodels/models/hardsphere.py @ 9eb3632

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 9eb3632 was 9eb3632, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

restructure kernels using fixed PD loops

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