source: sasmodels/sasmodels/models/pringle.c @ c047acf

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

much faster pringle model which can run on single precision GPU

  • Property mode set to 100644
File size: 3.5 KB
Line 
1double form_volume(double radius, double thickness, double alpha, double beta);
2
3double Iq(double q,
4          double radius,
5          double thickness,
6          double alpha,
7          double beta,
8          double sld_pringle,
9          double sld_solvent);
10
11
12static
13void _integrate_bessel(
14    double radius,
15    double alpha,
16    double beta,
17    double q_sin_psi,
18    double q_cos_psi,
19    double n,
20    double *Sn,
21    double *Cn)
22{
23    // translate gauss point z in [-1,1] to a point in [0, radius]
24    const double zm = 0.5*radius;
25    const double zb = 0.5*radius;
26
27    // evaluate at Gauss points
28    double sumS = 0.0;          // initialize integral
29    double sumC = 0.0;          // initialize integral
30    double r;
31    for (int i=0; i < 76; i++) {
32        r = Gauss76Z[i]*zm + zb;
33
34        const double qrs = r*q_sin_psi;
35        const double qrrc = r*r*q_cos_psi;
36
37        double y = Gauss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs);
38        double S, C;
39        SINCOS(alpha*qrrc, S, C);
40        sumS += y*S;
41        sumC += y*C;
42    }
43
44    #if 1
45    // TODO: should the normalization be to R^2 or the last value, r^2
46    // the gaussian window does not go all the way from 0 to 1.
47    //radius = Gauss76Z[75] * zm + zb;
48    *Sn = zm*sumS / (r*r);
49    *Cn = zm*sumC / (r*r);
50    #else
51    *Sn = zm*sumS / (radius*radius);
52    *Cn = zm*sumC / (radius*radius);
53    #endif
54}
55
56static
57double _sum_bessel_orders(
58    double radius,
59    double alpha,
60    double beta,
61    double q_sin_psi,
62    double q_cos_psi)
63{
64    //calculate sum term from n = -3 to 3
65    //Note 1:
66    //    S_n(-x) = (-1)^S_n(x)
67    //    => S_n^2(-x) = S_n^2(x),
68    //    => sum_-k^k Sk = S_0^2 + 2*sum_1^kSk^2
69    //Note 2:
70    //    better precision to sum terms from smaller to larger
71    //    though it doesn't seem to make a difference in this case.
72    double Sn, Cn, sum;
73    sum = 0.0;
74    for (int n=3; n>0; n--) {
75      _integrate_bessel(radius, alpha, beta, q_sin_psi, q_cos_psi, n, &Sn, &Cn);
76      sum += 2.0*(Sn*Sn + Cn*Cn);
77    }
78    _integrate_bessel(radius, alpha, beta, q_sin_psi, q_cos_psi, 0, &Sn, &Cn);
79    sum += Sn*Sn+ Cn*Cn;
80    return sum;
81}
82
83static
84double _integrate_psi(
85    double q,
86    double radius,
87    double thickness,
88    double alpha,
89    double beta)
90{
91    // translate gauss point z in [-1,1] to a point in [0, pi/2]
92    const double zm = M_PI_4;
93    const double zb = M_PI_4;
94
95    double sum = 0.0;
96    for (int i = 0; i < 76; i++) {
97        double psi = Gauss76Z[i]*zm + zb;
98        double sin_psi, cos_psi;
99        SINCOS(psi, sin_psi, cos_psi);
100        double bessel_term = _sum_bessel_orders(radius, alpha, beta, q*sin_psi, q*cos_psi);
101        double sinc_term = square(sinc(q * thickness * cos_psi / 2.0));
102        double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term;
103        sum += Gauss76Wt[i] * pringle_kernel;
104    }
105
106    return zm * sum;
107}
108
109double form_volume(double radius, double thickness, double alpha, double beta)
110{
111    // TODO: Normalize by form volume
112    //return M_PI*radius*radius*thickness;
113    return 1.0;
114}
115
116double Iq(
117    double q,
118    double radius,
119    double thickness,
120    double alpha,
121    double beta,
122    double sld_pringle,
123    double sld_solvent)
124{
125    double form = _integrate_psi(q, radius, thickness, alpha, beta);
126    double contrast = sld_pringle - sld_solvent;
127    double volume = M_PI*radius*radius*thickness;
128    // TODO: If normalize by form volume, need an extra volume here
129    //return 1.0e-4*form * square(contrast * volume);
130    return 1.0e-4*form * square(contrast) * volume;
131}
Note: See TracBrowser for help on using the repository browser.