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, 5 years ago

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

  • Property mode set to 100644
File size: 4.5 KB
RevLine 
[becded3]1static double
2form_volume(double length_a, double length_b, double length_c)
[c5b7d07]3{
[a807206]4    return length_a * length_b * length_c;
[c5b7d07]5}
6
[d277229]7static double
[99658f6]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
[d277229]39effective_radius(int mode, double length_a, double length_b, double length_c)
40{
[ee60aa7]41    switch (mode) {
[d42dd4a]42    default:
[99658f6]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
[6d5601c]46        return cbrt(length_a*length_b*length_c/M_4PI_3);
[99658f6]47    case 3: // half length_a
[d277229]48        return 0.5 * length_a;
[99658f6]49    case 4: // half length_b
[d277229]50        return 0.5 * length_b;
[99658f6]51    case 5: // half length_c
[d277229]52        return 0.5 * length_c;
[99658f6]53    case 6: // equivalent circular cross-section
[d277229]54        return sqrt(length_a*length_b/M_PI);
[99658f6]55    case 7: // half ab diagonal
[d277229]56        return 0.5*sqrt(length_a*length_a + length_b*length_b);
[99658f6]57    case 8: // half diagonal
[d277229]58        return 0.5*sqrt(length_a*length_a + length_b*length_b + length_c*length_c);
59    }
60}
61
[71b751d]62static void
63Fq(double q,
64    double *F1,
65    double *F2,
[c5b7d07]66    double sld,
67    double solvent_sld,
[a807206]68    double length_a,
69    double length_b,
70    double length_c)
[c5b7d07]71{
[14838a3]72    const double mu = 0.5 * q * length_b;
[2a0b2b1]73
[c5b7d07]74    // Scale sides by B
[14838a3]75    const double a_scaled = length_a / length_b;
76    const double c_scaled = length_c / length_b;
[2a0b2b1]77
[14838a3]78    // outer integral (with gauss points), integration limits = 0, 1
[71b751d]79    double outer_total_F1 = 0.0; //initialize integral
80    double outer_total_F2 = 0.0; //initialize integral
[74768cb]81    for( int i=0; i<GAUSS_N; i++) {
82        const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 );
[14838a3]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.
[71b751d]87        double inner_total_F1 = 0.0;
88        double inner_total_F2 = 0.0;
[74768cb]89        for(int j=0; j<GAUSS_N; j++) {
90            const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 );
[14838a3]91            double sin_uu, cos_uu;
92            SINCOS(M_PI_2*uu, sin_uu, cos_uu);
[1e7b0db0]93            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
94            const double si2 = sas_sinx_x(mu_proj * cos_uu);
[71b751d]95            const double fq = si1 * si2;
96            inner_total_F1 += GAUSS_W[j] * fq;
97            inner_total_F2 += GAUSS_W[j] * fq * fq;
[c5b7d07]98        }
[dbf1a60]99        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5
[71b751d]100        inner_total_F1 *= 0.5;
101        inner_total_F2 *= 0.5;
[c5b7d07]102
[1e7b0db0]103        const double si = sas_sinx_x(mu * c_scaled * sigma);
[71b751d]104        outer_total_F1 += GAUSS_W[i] * inner_total_F1 * si;
105        outer_total_F2 += GAUSS_W[i] * inner_total_F2 * si * si;
[14838a3]106    }
[dbf1a60]107    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5
[71b751d]108    outer_total_F1 *= 0.5;
109    outer_total_F2 *= 0.5;
[14838a3]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);
[71b751d]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;
[c5b7d07]117}
118
[becded3]119static double
[108e70e]120Iqabc(double qa, double qb, double qc,
[c5b7d07]121    double sld,
122    double solvent_sld,
[a807206]123    double length_a,
124    double length_b,
[becded3]125    double length_c)
[c5b7d07]126{
[2a0b2b1]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);
[14838a3]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;
[c5b7d07]135}
Note: See TracBrowser for help on using the repository browser.