source: sasmodels/sasmodels/models/bcc_paracrystal.c @ becded3

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

update oriented models to new interface (which will be in the next commit)

  • Property mode set to 100644
File size: 3.4 KB
Line 
1static double
2bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
3{
4#if 0  // Equations as written in Matsuoka
5    const double a1 = (+qa + qb + qc)/2.0;
6    const double a2 = (-qa - qb + qc)/2.0;
7    const double a3 = (-qa + qb - qc)/2.0;
8#else
9    const double a1 = (+qa + qb - qc)/2.0;
10    const double a2 = (+qa - qb + qc)/2.0;
11    const double a3 = (-qa + qb + qc)/2.0;
12#endif
13
14#if 1
15    // Numerator: (1 - exp(a)^2)^3
16    //         => (-(exp(2a) - 1))^3
17    //         => -expm1(2a)^3
18    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a))
19    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1)
20    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1)
21    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
22    const double exp_arg = exp(arg);
23    const double Zq = -cube(expm1(2.0*arg))
24        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0)
25          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0)
26          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0));
27#else
28    // Alternate form, which perhaps is more approachable
29    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
30    const double sinh_qd = sinh(arg);
31    const double cosh_qd = cosh(arg);
32    const double Zq = sinh_qd/(cosh_qd - cos(dnn*a1))
33                    * sinh_qd/(cosh_qd - cos(dnn*a2))
34                    * sinh_qd/(cosh_qd - cos(dnn*a3));
35#endif
36
37    return Zq;
38}
39
40
41// occupied volume fraction calculated from lattice symmetry and sphere radius
42static double
43bcc_volume_fraction(double radius, double dnn)
44{
45    return 2.0*sphere_volume(sqrt(0.75)*radius/dnn);
46}
47
48static double
49form_volume(double radius)
50{
51    return sphere_volume(radius);
52}
53
54
55static double Iq(double q, double dnn,
56    double d_factor, double radius,
57    double sld, double solvent_sld)
58{
59    // translate a point in [-1,1] to a point in [0, 2 pi]
60    const double phi_m = M_PI;
61    const double phi_b = M_PI;
62    // translate a point in [-1,1] to a point in [0, pi]
63    const double theta_m = M_PI_2;
64    const double theta_b = M_PI_2;
65
66    double outer_sum = 0.0;
67    for(int i=0; i<150; i++) {
68        double inner_sum = 0.0;
69        const double theta = Gauss150Z[i]*theta_m + theta_b;
70        double sin_theta, cos_theta;
71        SINCOS(theta, sin_theta, cos_theta);
72        const double qc = q*cos_theta;
73        const double qab = q*sin_theta;
74        for(int j=0;j<150;j++) {
75            const double phi = Gauss150Z[j]*phi_m + phi_b;
76            double sin_phi, cos_phi;
77            SINCOS(phi, sin_phi, cos_phi);
78            const double qa = qab*cos_phi;
79            const double qb = qab*sin_phi;
80            const double form = bcc_Zq(qa, qb, qc, dnn, d_factor);
81            inner_sum += Gauss150Wt[j] * form;
82        }
83        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
84        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
85    }
86    outer_sum *= theta_m;
87    const double Zq = outer_sum/(4.0*M_PI);
88    const double Pq = sphere_form(q, radius, sld, solvent_sld);
89    return bcc_volume_fraction(radius, dnn) * Pq * Zq;
90}
91
92
93static double Iqxy(double qa, double qb, double qc,
94    double dnn, double d_factor, double radius,
95    double sld, double solvent_sld)
96{
97    const double q = sqrt(qa*qa + qb*qb + qc*qc);
98    const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor);
99    const double Pq = sphere_form(q, radius, sld, solvent_sld);
100    return bcc_volume_fraction(radius, dnn) * Pq * Zq;
101}
Note: See TracBrowser for help on using the repository browser.