source: sasmodels/sasmodels/models/cylinder.c @ a34b811

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since a34b811 was a34b811, checked in by Paul Kienzle <pkienzle@…>, 7 months ago

use radius_effective/radius_effective_mode/radius_effective_modes consistently throughout the code

  • Property mode set to 100644
File size: 2.5 KB
Line 
1#define INVALID(v) (v.radius<0 || v.length<0)
2
3static double
4form_volume(double radius, double length)
5{
6    return M_PI*radius*radius*length;
7}
8
9static double
10_fq(double qab, double qc, double radius, double length)
11{
12    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length);
13}
14
15static double
16radius_from_excluded_volume(double radius, double length)
17{
18    return 0.5*cbrt(0.75*radius*(2.0*radius*length + (radius + length)*(M_PI*radius + length)));
19}
20
21static double
22radius_from_volume(double radius, double length)
23{
24    return cbrt(M_PI*radius*radius*length/M_4PI_3);
25}
26
27static double
28radius_from_diagonal(double radius, double length)
29{
30    return sqrt(radius*radius + 0.25*length*length);
31}
32
33static double
34radius_effective(int mode, double radius, double length)
35{
36    switch (mode) {
37    default:
38    case 1:
39        return radius_from_excluded_volume(radius, length);
40    case 2:
41        return radius_from_volume(radius, length);
42    case 3:
43        return radius;
44    case 4:
45        return 0.5*length;
46    case 5:
47        return (radius < 0.5*length ? radius : 0.5*length);
48    case 6:
49        return (radius > 0.5*length ? radius : 0.5*length);
50    case 7:
51        return radius_from_diagonal(radius,length);
52    }
53}
54
55static void
56Fq(double q,
57    double *F1,
58    double *F2,
59    double sld,
60    double solvent_sld,
61    double radius,
62    double length)
63{
64    // translate a point in [-1,1] to a point in [0, pi/2]
65    const double zm = M_PI_4;
66    const double zb = M_PI_4;
67
68    double total_F1 = 0.0;
69    double total_F2 = 0.0;
70    for (int i=0; i<GAUSS_N ;i++) {
71        const double theta = GAUSS_Z[i]*zm + zb;
72        double sin_theta, cos_theta; // slots to hold sincos function output
73        // theta (theta,phi) the projection of the cylinder on the detector plane
74        SINCOS(theta , sin_theta, cos_theta);
75        const double form = _fq(q*sin_theta, q*cos_theta, radius, length);
76        total_F1 += GAUSS_W[i] * form * sin_theta;
77        total_F2 += GAUSS_W[i] * form * form * sin_theta;
78    }
79    // translate dx in [-1,1] to dx in [lower,upper]
80    total_F1 *= zm;
81    total_F2 *= zm;
82    const double s = (sld - solvent_sld) * form_volume(radius, length);
83    *F1 = 1e-2 * s * total_F1;
84    *F2 = 1e-4 * s * s * total_F2;
85}
86
87
88
89static double
90Iqac(double qab, double qc,
91    double sld,
92    double solvent_sld,
93    double radius,
94    double length)
95{
96    const double form = _fq(qab, qc, radius, length);
97    const double s = (sld-solvent_sld) * form_volume(radius, length);
98    return 1.0e-4 * square(s * form);
99}
100
Note: See TracBrowser for help on using the repository browser.