source: sasmodels/sasmodels/models/parallelepiped.c @ 14838a3

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 14838a3 was 14838a3, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

update parallelepiped and core_shell_parallelepiped to use ORIENT macros

  • Property mode set to 100644
File size: 2.6 KB
Line 
1double form_volume(double length_a, double length_b, double length_c);
2double Iq(double q, double sld, double solvent_sld,
3    double length_a, double length_b, double length_c);
4double Iqxy(double qx, double qy, double sld, double solvent_sld,
5    double length_a, double length_b, double length_c,
6    double theta, double phi, double psi);
7
8double form_volume(double length_a, double length_b, double length_c)
9{
10    return length_a * length_b * length_c;
11}
12
13
14double Iq(double q,
15    double sld,
16    double solvent_sld,
17    double length_a,
18    double length_b,
19    double length_c)
20{
21    const double mu = 0.5 * q * length_b;
22   
23    // Scale sides by B
24    const double a_scaled = length_a / length_b;
25    const double c_scaled = length_c / length_b;
26       
27    // outer integral (with gauss points), integration limits = 0, 1
28    double outer_total = 0; //initialize integral
29
30    for( int i=0; i<76; i++) {
31        const double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
32        const double mu_proj = mu * sqrt(1.0-sigma*sigma);
33
34        // inner integral (with gauss points), integration limits = 0, 1
35        // corresponding to angles from 0 to pi/2.
36        double inner_total = 0.0;
37        for(int j=0; j<76; j++) {
38            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
39            double sin_uu, cos_uu;
40            SINCOS(M_PI_2*uu, sin_uu, cos_uu);
41            const double si1 = sinc(mu_proj * sin_uu * a_scaled);
42            const double si2 = sinc(mu_proj * cos_uu);
43            inner_total += Gauss76Wt[j] * square(si1 * si2);
44        }
45        inner_total *= 0.5;
46
47        const double si = sinc(mu * c_scaled * sigma);
48        outer_total += Gauss76Wt[i] * inner_total * si * si;
49    }
50    outer_total *= 0.5;
51
52    // Multiply by contrast^2 and convert from [1e-12 A-1] to [cm-1]
53    const double V = form_volume(length_a, length_b, length_c);
54    const double drho = (sld-solvent_sld);
55    return 1.0e-4 * square(drho * V) * outer_total;
56}
57
58
59double Iqxy(double qx, double qy,
60    double sld,
61    double solvent_sld,
62    double length_a,
63    double length_b,
64    double length_c,
65    double theta,
66    double phi,
67    double psi)
68{
69    double q, cos_val_a, cos_val_b, cos_val_c;
70    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a);
71
72    const double siA = sinc(0.5*q*length_a*cos_val_a);
73    const double siB = sinc(0.5*q*length_b*cos_val_b);
74    const double siC = sinc(0.5*q*length_c*cos_val_c);
75    const double V = form_volume(length_a, length_b, length_c);
76    const double drho = (sld - solvent_sld);
77    const double form = V * drho * siA * siB * siC;
78    // Square and convert from [1e-12 A-1] to [cm-1]
79    return 1.0e-4 * form * form;
80}
Note: See TracBrowser for help on using the repository browser.