source: sasmodels/sasmodels/models/parallelepiped.c @ dbf1a60

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

Add comments to c code and clean documentatin

Added comments to c code in both parallelepiped and core shell
parallelepiped noting the change of integration varialbes in the
computation. Cleaned up and final corrections to the core shell
documentation and did some cleaning of parallelipiped. In particular
tried to bring a bit more consistency between the docs.

addresses #896

  • Property mode set to 100644
File size: 2.3 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        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5
41        inner_total *= 0.5;
42
43        const double si = sas_sinx_x(mu * c_scaled * sigma);
44        outer_total += GAUSS_W[i] * inner_total * si * si;
45    }
46    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5
47    outer_total *= 0.5;
48
49    // Multiply by contrast^2 and convert from [1e-12 A-1] to [cm-1]
50    const double V = form_volume(length_a, length_b, length_c);
51    const double drho = (sld-solvent_sld);
52    return 1.0e-4 * square(drho * V) * outer_total;
53}
54
55
56static double
57Iqabc(double qa, double qb, double qc,
58    double sld,
59    double solvent_sld,
60    double length_a,
61    double length_b,
62    double length_c)
63{
64    const double siA = sas_sinx_x(0.5*length_a*qa);
65    const double siB = sas_sinx_x(0.5*length_b*qb);
66    const double siC = sas_sinx_x(0.5*length_c*qc);
67    const double V = form_volume(length_a, length_b, length_c);
68    const double drho = (sld - solvent_sld);
69    const double form = V * drho * siA * siB * siC;
70    // Square and convert from [1e-12 A-1] to [cm-1]
71    return 1.0e-4 * form * form;
72}
Note: See TracBrowser for help on using the repository browser.