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@…>, 6 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
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
[deb7ee0]7
[becded3]8static double
9Iq(double q,
[c5b7d07]10    double sld,
11    double solvent_sld,
[a807206]12    double length_a,
13    double length_b,
14    double length_c)
[c5b7d07]15{
[14838a3]16    const double mu = 0.5 * q * length_b;
[2a0b2b1]17
[c5b7d07]18    // Scale sides by B
[14838a3]19    const double a_scaled = length_a / length_b;
20    const double c_scaled = length_c / length_b;
[2a0b2b1]21
[14838a3]22    // outer integral (with gauss points), integration limits = 0, 1
23    double outer_total = 0; //initialize integral
24
[74768cb]25    for( int i=0; i<GAUSS_N; i++) {
26        const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 );
[14838a3]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;
[74768cb]32        for(int j=0; j<GAUSS_N; j++) {
33            const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 );
[14838a3]34            double sin_uu, cos_uu;
35            SINCOS(M_PI_2*uu, sin_uu, cos_uu);
[1e7b0db0]36            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
37            const double si2 = sas_sinx_x(mu_proj * cos_uu);
[74768cb]38            inner_total += GAUSS_W[j] * square(si1 * si2);
[c5b7d07]39        }
[14838a3]40        inner_total *= 0.5;
[c5b7d07]41
[1e7b0db0]42        const double si = sas_sinx_x(mu * c_scaled * sigma);
[74768cb]43        outer_total += GAUSS_W[i] * inner_total * si * si;
[14838a3]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;
[c5b7d07]51}
52
53
[becded3]54static double
[108e70e]55Iqabc(double qa, double qb, double qc,
[c5b7d07]56    double sld,
57    double solvent_sld,
[a807206]58    double length_a,
59    double length_b,
[becded3]60    double length_c)
[c5b7d07]61{
[2a0b2b1]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);
[14838a3]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;
[c5b7d07]70}
Note: See TracBrowser for help on using the repository browser.