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

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since e5cb3df was 99658f6, checked in by grethevj, 5 years ago

updated ER functions including cylinder excluded volume, to match 4.x

  • Property mode set to 100644
File size: 3.6 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,
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 < GAUSS_N; i++) {
32        r = GAUSS_Z[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 = GAUSS_W[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 < GAUSS_N; i++) {
89        double psi = GAUSS_Z[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(sas_sinx_x(q * thickness * cos_psi / 2.0));
94        double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term;
95        sum += GAUSS_W[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
106static double
107radius_from_excluded_volume(double radius, double thickness)
108{
109    return 0.5*cbrt(0.75*radius*(2.0*radius*thickness + (radius + thickness)*(M_PI*radius + thickness)));
110}
111
112static double
113effective_radius(int mode, double radius, double thickness, double alpha, double beta)
114{
115    switch (mode) {
116    default:
117    case 1: // equivalent cylinder excluded volume
118        return radius_from_excluded_volume(radius, thickness);
119    case 2: // equivalent volume sphere
120        return cbrt(M_PI*radius*radius*thickness/M_4PI_3);
121    case 3: // radius
122        return radius;
123    }
124}
125
126double Iq(
127    double q,
128    double radius,
129    double thickness,
130    double alpha,
131    double beta,
132    double sld,
133    double sld_solvent)
134{
135    double form = _integrate_psi(q, radius, thickness, alpha, beta);
136    double contrast = sld - sld_solvent;
137    double volume = M_PI*radius*radius*thickness;
138    return 1.0e-4*form * square(contrast * volume);
139}
Note: See TracBrowser for help on using the repository browser.