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
RevLine 
[301e096]1# Note: model title and parameter table are inserted automatically
[529b8b4]2r"""Calculate the interparticle structure factor for monodisperse
3spherical particles interacting through hard sphere (excluded volume)
4interactions.
[40a87fa]5May be a reasonable approximation for other shapes of particles that
6freely rotate, and for moderately polydisperse systems. Though strictly
[d529d93]7the maths needs to be modified (no \Beta(Q) correction yet in sasview).
[301e096]8
[d529d93]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
[40a87fa]15For numerical stability the computation uses a Taylor series expansion
[d529d93]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
[529b8b4]20potential is
[301e096]21
[eb69cce]22.. math::
[529b8b4]23
24    U(r) = \begin{cases}
25    \infty & r < 2R \\
26    0 & r \geq 2R
27    \end{cases}
[301e096]28
[eb69cce]29where $r$ is the distance from the center of the sphere of a radius $R$.
[301e096]30
31For a 2D plot, the wave transfer is defined as
32
33.. math::
34
[529b8b4]35    q = \sqrt{q_x^2 + q_y^2}
[301e096]36
37
[eb69cce]38References
39----------
[301e096]40
41J K Percus, J Yevick, *J. Phys. Rev.*, 110, (1958) 1
42"""
43
[2d81cfe]44import numpy as np
[3c56da87]45from numpy import inf
[301e096]46
[034e19a]47name = "hardsphere"
48title = "Hard sphere structure factor, with Percus-Yevick closure"
[301e096]49description = """\
[3e428ec]50    [Hard sphere structure factor, with Percus-Yevick closure]
[301e096]51        Interparticle S(Q) for random, non-interacting spheres.
[3e428ec]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.
[d529d93]56    radius_effective is the hard sphere radius
[3e428ec]57    volfraction is the volume fraction occupied by the spheres.
[301e096]58"""
[a5d0d00]59category = "structure-factor"
[8e45182]60structure_factor = True
[7f1ee79]61single = False # TODO: check
[301e096]62
[3e428ec]63#             ["name", "units", default, [lower, upper], "type","description"],
[71b751d]64parameters = [["radius_effective", "Ang", 50.0, [0, inf], "",
[3e428ec]65               "effective radius of hard sphere"],
66              ["volfraction", "", 0.2, [0, 0.74], "",
67               "volume fraction of hard spheres"],
68             ]
[301e096]69
[9eb3632]70Iq = r"""
[f0fb9fe]71      double D,A,B,G,X,X2,X4,S,C,FF,HARDSPH;
[8f04da4]72      // these are c compiler instructions, can also put normal code inside the "if else" structure
[934f906]73      #if FLOAT_SIZE > 4
[2d81cfe]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
[934f906]78      #define CUTOFFHS 0.05
79      #else
80      // 0.1 bad, 0.2 OK, 0.3 good, 0.4 better, 0.8 no good
[8f04da4]81      #define CUTOFFHS 0.4
[934f906]82      #endif
[f0fb9fe]83
[d529d93]84      if(fabs(radius_effective) < 1.E-12) {
[f0fb9fe]85               HARDSPH=1.0;
[9eb3632]86//printf("HS1 %g: %g\n",q,HARDSPH);
[f0fb9fe]87               return(HARDSPH);
88      }
[2d81cfe]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)
[97e6d3c]91      X = 1.0/( 1.0 -volfraction);
92      D= X*X;
93      A= (1.+2.*volfraction)*D;
94      A *=A;
[d529d93]95      X=fabs(q*radius_effective*2.0);
[f0fb9fe]96
97      if(X < 5.E-06) {
98                 HARDSPH=1./A;
[9eb3632]99//printf("HS2 %g: %g\n",q,HARDSPH);
[f0fb9fe]100                 return(HARDSPH);
101      }
[97e6d3c]102      X2 =X*X;
103      B = (1.0 +0.5*volfraction)*D;
104      B *= B;
105      B *= -6.*volfraction;
[f0fb9fe]106      G=0.5*volfraction*A;
107
[934f906]108      if(X < CUTOFFHS) {
109      // RKH Feb 2016, use Taylor series expansion for small X
[2d81cfe]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
[8f04da4]114      // for 5 or 6 significant figures, but I put the X^4 one in anyway
[97e6d3c]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
[f0fb9fe]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
[97e6d3c]124            //FF = (8 +2.*volfraction + ( volfraction/4. -0.8 +(volfraction/100. -1./35.)*X2 )*X2 )*A  + (3.0 -X2/3. +X4/40.)*2.*B;
[f0fb9fe]125            HARDSPH= 1./(1. + volfraction*FF );
[9eb3632]126//printf("HS3 %g: %g\n",q,HARDSPH);
[f0fb9fe]127            return(HARDSPH);
128      }
[97e6d3c]129      X4=X2*X2;
[f0fb9fe]130      SINCOS(X,S,C);
131
[2d81cfe]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.
[97e6d3c]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 ;
[f0fb9fe]138      HARDSPH= 1./(1. + 24.*volfraction*FF/X2 );
139
[97e6d3c]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
[f0fb9fe]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 );
[9eb3632]152//printf("HS4 %g: %g\n",q,HARDSPH);
[f0fb9fe]153      return(HARDSPH);
[301e096]154   """
155
[8f04da4]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
[2d81cfe]164# Q=0.001 is in the Taylor series, low Q part, so add Q=0.1,
165# assuming double precision sasview is correct
[7f47777]166tests = [
[2d81cfe]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.