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

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since d277229 was d277229, checked in by grethevj, 6 years ago

Models updated to include choices for effective interaction radii

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