source: sasmodels/sasmodels/models/sc_paracrystal.c @ 2a0b2b1

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

restructure all 2D models to work with (qa,qb,qc) = rotate(qx,qy) rather than working with angles directly in preparation for revised jitter algorithm

  • Property mode set to 100644
File size: 3.0 KB
Line 
1static double
2_sq_sc(double qa, double qb, double qc, double dnn, double d_factor)
3{
4    // Rewriting equations for efficiency, accuracy and readability, and so
5    // code is reusable between 1D and 2D models.
6    const double a1 = qa;
7    const double a2 = qb;
8    const double a3 = qc;
9
10    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
11
12    // Numerator: (1 - exp(a)^2)^3
13    //         => (-(exp(2a) - 1))^3
14    //         => -expm1(2a)^3
15    // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2)
16    //         => exp(a)^2 - 2 cos(xk) exp(a) + 1
17    //         => (exp(a) - 2 cos(xk)) * exp(a) + 1
18    const double exp_arg = exp(-arg);
19    const double Sq = -cube(expm1(-2.0*arg))
20        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0)
21          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0)
22          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0));
23
24    return Sq;
25}
26
27// occupied volume fraction calculated from lattice symmetry and sphere radius
28static double
29_sc_volume_fraction(double radius, double dnn)
30{
31    return sphere_volume(radius/dnn);
32}
33
34static double
35form_volume(double radius)
36{
37    return sphere_volume(radius);
38}
39
40
41static double Iq(double q, double dnn,
42  double d_factor, double radius,
43  double sld, double solvent_sld)
44{
45    // translate a point in [-1,1] to a point in [0, 2 pi]
46    const double phi_m = M_PI_4;
47    const double phi_b = M_PI_4;
48    // translate a point in [-1,1] to a point in [0, pi]
49    const double theta_m = M_PI_4;
50    const double theta_b = M_PI_4;
51
52
53    double outer_sum = 0.0;
54    for(int i=0; i<150; i++) {
55        double inner_sum = 0.0;
56        const double theta = Gauss150Z[i]*theta_m + theta_b;
57        double sin_theta, cos_theta;
58        SINCOS(theta, sin_theta, cos_theta);
59        const double qc = q*cos_theta;
60        const double qab = q*sin_theta;
61        for(int j=0;j<150;j++) {
62            const double phi = Gauss150Z[j]*phi_m + phi_b;
63            double sin_phi, cos_phi;
64            SINCOS(phi, sin_phi, cos_phi);
65            const double qa = qab*cos_phi;
66            const double qb = qab*sin_phi;
67            const double fq = _sq_sc(qa, qb, qc, dnn, d_factor);
68            inner_sum += Gauss150Wt[j] * fq;
69        }
70        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
71        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
72    }
73    outer_sum *= theta_m;
74    const double Sq = outer_sum/M_PI_2;
75    const double Pq = sphere_form(q, radius, sld, solvent_sld);
76
77    return _sc_volume_fraction(radius, dnn) * Pq * Sq;
78}
79
80
81static double Iqxy(double qx, double qy,
82    double dnn, double d_factor, double radius,
83    double sld, double solvent_sld,
84    double theta, double phi, double psi)
85{
86    double q, zhat, yhat, xhat;
87    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
88    const double qa = q*xhat;
89    const double qb = q*yhat;
90    const double qc = q*zhat;
91
92    q = sqrt(qa*qa + qb*qb + qc*qc);
93    const double Pq = sphere_form(q, radius, sld, solvent_sld);
94    const double Sq = _sq_sc(qa, qb, qc, dnn, d_factor);
95    return _sc_volume_fraction(radius, dnn) * Pq * Sq;
96}
Note: See TracBrowser for help on using the repository browser.