source: sasmodels/sasmodels/models/parallelepiped.c @ 99658f6

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 99658f6 was 99658f6, checked in by grethevj, 8 months ago

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

  • Property mode set to 100644
File size: 4.5 KB
Line 
1static double
2form_volume(double length_a, double length_b, double length_c)
3{
4    return length_a * length_b * length_c;
5}
6
7static double
8radius_from_excluded_volume(double length_a, double length_b, double length_c)
9{
10    double r_equiv, length;
11    double lengths[3] = {length_a, length_b, length_c};
12    double lengthmax = fmax(lengths[0],fmax(lengths[1],lengths[2]));
13    double length_1 = lengthmax;
14    double length_2 = lengthmax;
15    double length_3 = lengthmax;
16
17    for(int ilen=0; ilen<3; ilen++) {
18        if (lengths[ilen] < length_1) {
19            length_2 = length_1;
20            length_1 = lengths[ilen];
21            } else {
22                if (lengths[ilen] < length_2) {
23                        length_2 = lengths[ilen];
24                }
25            }
26    }
27    if(length_2-length_1 > length_3-length_2) {
28        r_equiv = sqrt(length_2*length_3/M_PI);
29        length  = length_1;
30    } else  {
31        r_equiv = sqrt(length_1*length_2/M_PI);
32        length  = length_3;
33    }
34
35    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length + (r_equiv + length)*(M_PI*r_equiv + length)));
36}
37
38static double
39effective_radius(int mode, double length_a, double length_b, double length_c)
40{
41    switch (mode) {
42    default:
43    case 1: // equivalent cylinder excluded volume
44        return radius_from_excluded_volume(length_a,length_b,length_c);
45    case 2: // equivalent volume sphere
46        return cbrt(length_a*length_b*length_c/M_4PI_3);
47    case 3: // half length_a
48        return 0.5 * length_a;
49    case 4: // half length_b
50        return 0.5 * length_b;
51    case 5: // half length_c
52        return 0.5 * length_c;
53    case 6: // equivalent circular cross-section
54        return sqrt(length_a*length_b/M_PI);
55    case 7: // half ab diagonal
56        return 0.5*sqrt(length_a*length_a + length_b*length_b);
57    case 8: // half diagonal
58        return 0.5*sqrt(length_a*length_a + length_b*length_b + length_c*length_c);
59    }
60}
61
62static void
63Fq(double q,
64    double *F1,
65    double *F2,
66    double sld,
67    double solvent_sld,
68    double length_a,
69    double length_b,
70    double length_c)
71{
72    const double mu = 0.5 * q * length_b;
73
74    // Scale sides by B
75    const double a_scaled = length_a / length_b;
76    const double c_scaled = length_c / length_b;
77
78    // outer integral (with gauss points), integration limits = 0, 1
79    double outer_total_F1 = 0.0; //initialize integral
80    double outer_total_F2 = 0.0; //initialize integral
81    for( int i=0; i<GAUSS_N; i++) {
82        const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 );
83        const double mu_proj = mu * sqrt(1.0-sigma*sigma);
84
85        // inner integral (with gauss points), integration limits = 0, 1
86        // corresponding to angles from 0 to pi/2.
87        double inner_total_F1 = 0.0;
88        double inner_total_F2 = 0.0;
89        for(int j=0; j<GAUSS_N; j++) {
90            const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 );
91            double sin_uu, cos_uu;
92            SINCOS(M_PI_2*uu, sin_uu, cos_uu);
93            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
94            const double si2 = sas_sinx_x(mu_proj * cos_uu);
95            const double fq = si1 * si2;
96            inner_total_F1 += GAUSS_W[j] * fq;
97            inner_total_F2 += GAUSS_W[j] * fq * fq;
98        }
99        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5
100        inner_total_F1 *= 0.5;
101        inner_total_F2 *= 0.5;
102
103        const double si = sas_sinx_x(mu * c_scaled * sigma);
104        outer_total_F1 += GAUSS_W[i] * inner_total_F1 * si;
105        outer_total_F2 += GAUSS_W[i] * inner_total_F2 * si * si;
106    }
107    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5
108    outer_total_F1 *= 0.5;
109    outer_total_F2 *= 0.5;
110
111    // Multiply by contrast^2 and convert from [1e-12 A-1] to [cm-1]
112    const double V = form_volume(length_a, length_b, length_c);
113    const double contrast = (sld-solvent_sld);
114    const double s = contrast * V;
115    *F1 = 1.0e-2 * s * outer_total_F1;
116    *F2 = 1.0e-4 * s * s * outer_total_F2;
117}
118
119static double
120Iqabc(double qa, double qb, double qc,
121    double sld,
122    double solvent_sld,
123    double length_a,
124    double length_b,
125    double length_c)
126{
127    const double siA = sas_sinx_x(0.5*length_a*qa);
128    const double siB = sas_sinx_x(0.5*length_b*qb);
129    const double siC = sas_sinx_x(0.5*length_c*qc);
130    const double V = form_volume(length_a, length_b, length_c);
131    const double drho = (sld - solvent_sld);
132    const double form = V * drho * siA * siB * siC;
133    // Square and convert from [1e-12 A-1] to [cm-1]
134    return 1.0e-4 * form * form;
135}
Note: See TracBrowser for help on using the repository browser.