source: sasmodels/sasmodels/models/pringle.c @ 58c3367

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 58c3367 was 58c3367, checked in by ajj, 5 years ago

Fixing volume normalisation and docs on pringle model.

Closes #604

  • Property mode set to 100644
File size: 3.0 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    *Sn = zm*sumS / (radius*radius);
45    *Cn = zm*sumC / (radius*radius);
46}
47
48static
49double _sum_bessel_orders(
50    double radius,
51    double alpha,
52    double beta,
53    double q_sin_psi,
54    double q_cos_psi)
55{
56    //calculate sum term from n = -3 to 3
57    //Note 1:
58    //    S_n(-x) = (-1)^S_n(x)
59    //    => S_n^2(-x) = S_n^2(x),
60    //    => sum_-k^k Sk = S_0^2 + 2*sum_1^kSk^2
61    //Note 2:
62    //    better precision to sum terms from smaller to larger
63    //    though it doesn't seem to make a difference in this case.
64    double Sn, Cn, sum;
65    sum = 0.0;
66    for (int n=3; n>0; n--) {
67      _integrate_bessel(radius, alpha, beta, q_sin_psi, q_cos_psi, n, &Sn, &Cn);
68      sum += 2.0*(Sn*Sn + Cn*Cn);
69    }
70    _integrate_bessel(radius, alpha, beta, q_sin_psi, q_cos_psi, 0, &Sn, &Cn);
71    sum += Sn*Sn+ Cn*Cn;
72    return sum;
73}
74
75static
76double _integrate_psi(
77    double q,
78    double radius,
79    double thickness,
80    double alpha,
81    double beta)
82{
83    // translate gauss point z in [-1,1] to a point in [0, pi/2]
84    const double zm = M_PI_4;
85    const double zb = M_PI_4;
86
87    double sum = 0.0;
88    for (int i = 0; i < 76; i++) {
89        double psi = Gauss76Z[i]*zm + zb;
90        double sin_psi, cos_psi;
91        SINCOS(psi, sin_psi, cos_psi);
92        double bessel_term = _sum_bessel_orders(radius, alpha, beta, q*sin_psi, q*cos_psi);
93        double sinc_term = square(sinc(q * thickness * cos_psi / 2.0));
94        double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term;
95        sum += Gauss76Wt[i] * pringle_kernel;
96    }
97
98    return zm * sum;
99}
100
101double form_volume(double radius, double thickness, double alpha, double beta)
102{
103    return M_PI*radius*radius*thickness;
104}
105
106double Iq(
107    double q,
108    double radius,
109    double thickness,
110    double alpha,
111    double beta,
112    double sld_pringle,
113    double sld_solvent)
114{
115    double form = _integrate_psi(q, radius, thickness, alpha, beta);
116    double contrast = sld_pringle - sld_solvent;
117    double volume = M_PI*radius*radius*thickness;
118    return 1.0e-4*form * square(contrast * volume);
119}
Note: See TracBrowser for help on using the repository browser.