source: sasmodels/sasmodels/models/bcc_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@…>, 22 months ago

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

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