source: sasmodels/sasmodels/models/parallelepiped.c @ 108e70e

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 108e70e was 108e70e, checked in by Paul Kienzle <pkienzle@…>, 4 years ago

use Iqac/Iqabc? for the new orientation interface but Iqxy for the old

  • Property mode set to 100644
File size: 2.2 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
7
8static double
9Iq(double q,
10    double sld,
11    double solvent_sld,
12    double length_a,
13    double length_b,
14    double length_c)
15{
16    const double mu = 0.5 * q * length_b;
17
18    // Scale sides by B
19    const double a_scaled = length_a / length_b;
20    const double c_scaled = length_c / length_b;
21
22    // outer integral (with gauss points), integration limits = 0, 1
23    double outer_total = 0; //initialize integral
24
25    for( int i=0; i<GAUSS_N; i++) {
26        const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 );
27        const double mu_proj = mu * sqrt(1.0-sigma*sigma);
28
29        // inner integral (with gauss points), integration limits = 0, 1
30        // corresponding to angles from 0 to pi/2.
31        double inner_total = 0.0;
32        for(int j=0; j<GAUSS_N; j++) {
33            const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 );
34            double sin_uu, cos_uu;
35            SINCOS(M_PI_2*uu, sin_uu, cos_uu);
36            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
37            const double si2 = sas_sinx_x(mu_proj * cos_uu);
38            inner_total += GAUSS_W[j] * square(si1 * si2);
39        }
40        inner_total *= 0.5;
41
42        const double si = sas_sinx_x(mu * c_scaled * sigma);
43        outer_total += GAUSS_W[i] * inner_total * si * si;
44    }
45    outer_total *= 0.5;
46
47    // Multiply by contrast^2 and convert from [1e-12 A-1] to [cm-1]
48    const double V = form_volume(length_a, length_b, length_c);
49    const double drho = (sld-solvent_sld);
50    return 1.0e-4 * square(drho * V) * outer_total;
51}
52
53
54static double
55Iqabc(double qa, double qb, double qc,
56    double sld,
57    double solvent_sld,
58    double length_a,
59    double length_b,
60    double length_c)
61{
62    const double siA = sas_sinx_x(0.5*length_a*qa);
63    const double siB = sas_sinx_x(0.5*length_b*qb);
64    const double siC = sas_sinx_x(0.5*length_c*qc);
65    const double V = form_volume(length_a, length_b, length_c);
66    const double drho = (sld - solvent_sld);
67    const double form = V * drho * siA * siB * siC;
68    // Square and convert from [1e-12 A-1] to [cm-1]
69    return 1.0e-4 * form * form;
70}
Note: See TracBrowser for help on using the repository browser.