source: sasmodels/sasmodels/models/sc_paracrystal.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@…>, 3 years ago

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

  • Property mode set to 100644
File size: 2.8 KB
Line 
1static double
2sc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
3{
4    // Equations from Matsuoka 9-10-11, multiplied by |q|
5    const double a1 = qa;
6    const double a2 = qb;
7    const double a3 = qc;
8
9    // Matsuoka 13-14-15
10    //     Z_k numerator: 1 - exp(a)^2
11    //     Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a)
12    // Rewriting numerator
13    //         => -(exp(2a) - 1)
14    //         => -expm1(2a)
15    // Rewriting denominator
16    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1)
17    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1
18    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
19    const double exp_arg = exp(arg);
20    const double Zq = -cube(expm1(2.0*arg))
21        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0)
22          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0)
23          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0));
24
25    return Zq;
26}
27
28// occupied volume fraction calculated from lattice symmetry and sphere radius
29static double
30sc_volume_fraction(double radius, double dnn)
31{
32    return sphere_volume(radius/dnn);
33}
34
35static double
36form_volume(double radius)
37{
38    return sphere_volume(radius);
39}
40
41
42static double
43Iq(double q, double dnn,
44    double d_factor, double radius,
45    double sld, double solvent_sld)
46{
47    // translate a point in [-1,1] to a point in [0, 2 pi]
48    const double phi_m = M_PI_4;
49    const double phi_b = M_PI_4;
50    // translate a point in [-1,1] to a point in [0, pi]
51    const double theta_m = M_PI_4;
52    const double theta_b = M_PI_4;
53
54
55    double outer_sum = 0.0;
56    for(int i=0; i<GAUSS_N; i++) {
57        double inner_sum = 0.0;
58        const double theta = GAUSS_Z[i]*theta_m + theta_b;
59        double sin_theta, cos_theta;
60        SINCOS(theta, sin_theta, cos_theta);
61        const double qc = q*cos_theta;
62        const double qab = q*sin_theta;
63        for(int j=0;j<GAUSS_N;j++) {
64            const double phi = GAUSS_Z[j]*phi_m + phi_b;
65            double sin_phi, cos_phi;
66            SINCOS(phi, sin_phi, cos_phi);
67            const double qa = qab*cos_phi;
68            const double qb = qab*sin_phi;
69            const double form = sc_Zq(qa, qb, qc, dnn, d_factor);
70            inner_sum += GAUSS_W[j] * form;
71        }
72        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
73        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
74    }
75    outer_sum *= theta_m;
76    const double Zq = outer_sum/M_PI_2;
77    const double Pq = sphere_form(q, radius, sld, solvent_sld);
78
79    return sc_volume_fraction(radius, dnn) * Pq * Zq;
80}
81
82
83static double
84Iqabc(double qa, double qb, double qc,
85    double dnn, double d_factor, double radius,
86    double sld, double solvent_sld)
87{
88    const double q = sqrt(qa*qa + qb*qb + qc*qc);
89    const double Pq = sphere_form(q, radius, sld, solvent_sld);
90    const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor);
91    return sc_volume_fraction(radius, dnn) * Pq * Zq;
92}
Note: See TracBrowser for help on using the repository browser.