source: sasmodels/sasmodels/models/pringle.c @ 30fbe2e

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 30fbe2e was 30fbe2e, checked in by butler, 7 years ago

Last of parameter normalization. Fixes #649

  • Property mode set to 100644
File size: 3.0 KB
RevLine 
[c047acf]1double form_volume(double radius, double thickness, double alpha, double beta);
[9da407c]2
3double Iq(double q,
4          double radius,
5          double thickness,
6          double alpha,
7          double beta,
[30fbe2e]8          double sld,
[9da407c]9          double sld_solvent);
10
11
12static
[c047acf]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;
[9da407c]26
27    // evaluate at Gauss points
[c047acf]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    }
[9da407c]43
[c047acf]44    *Sn = zm*sumS / (radius*radius);
45    *Cn = zm*sumC / (radius*radius);
[9da407c]46}
47
48static
[c047acf]49double _sum_bessel_orders(
50    double radius,
51    double alpha,
52    double beta,
53    double q_sin_psi,
54    double q_cos_psi)
55{
[9da407c]56    //calculate sum term from n = -3 to 3
[c047acf]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);
[9da407c]69    }
[c047acf]70    _integrate_bessel(radius, alpha, beta, q_sin_psi, q_cos_psi, 0, &Sn, &Cn);
71    sum += Sn*Sn+ Cn*Cn;
72    return sum;
[9da407c]73}
74
[c047acf]75static
76double _integrate_psi(
77    double q,
78    double radius,
79    double thickness,
80    double alpha,
81    double beta)
[9da407c]82{
[c047acf]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;
[9da407c]96    }
97
[c047acf]98    return zm * sum;
[9da407c]99}
100
[c047acf]101double form_volume(double radius, double thickness, double alpha, double beta)
102{
[58c3367]103    return M_PI*radius*radius*thickness;
[9da407c]104}
105
[c047acf]106double Iq(
107    double q,
108    double radius,
109    double thickness,
110    double alpha,
111    double beta,
[30fbe2e]112    double sld,
[c047acf]113    double sld_solvent)
[9da407c]114{
[c047acf]115    double form = _integrate_psi(q, radius, thickness, alpha, beta);
[30fbe2e]116    double contrast = sld - sld_solvent;
[c047acf]117    double volume = M_PI*radius*radius*thickness;
[58c3367]118    return 1.0e-4*form * square(contrast * volume);
[9da407c]119}
Note: See TracBrowser for help on using the repository browser.