source: sasmodels/sasmodels/models/ellipsoid.c @ 1e7b0db0

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 1e7b0db0 was 73e08ae, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

core_shell_ellipsoid, cylinder, ellipsoid: code cleanup

  • Property mode set to 100644
File size: 2.5 KB
Line 
1double form_volume(double radius_polar, double radius_equatorial);
2double Iq(double q, double sld, double sld_solvent, double radius_polar, double radius_equatorial);
3double Iqxy(double qx, double qy, double sld, double sld_solvent,
4    double radius_polar, double radius_equatorial, double theta, double phi);
5
6double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha);
7double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha)
8{
9    double ratio = radius_polar/radius_equatorial;
10    // Given the following under the radical:
11    //     1 + sin^2(T) (v^2 - 1)
12    // we can expand to match the form given in Guinier (1955)
13    //     = (1 - sin^2(T)) + v^2 sin^2(T) = cos^2(T) + sin^2(T)
14    // Instead of using pythagoras we could pass in sin and cos; this may be
15    // slightly better for 2D which has already computed it, but it introduces
16    // an extra sqrt and square for 1-D not required by the current form, so
17    // leave it as is.
18    const double r = radius_equatorial
19                     * sqrt(1.0 + sin_alpha*sin_alpha*(ratio*ratio - 1.0));
20    const double f = sph_j1c(q*r);
21
22    return f*f;
23}
24
25double form_volume(double radius_polar, double radius_equatorial)
26{
27    return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial;
28}
29
30double Iq(double q,
31    double sld,
32    double sld_solvent,
33    double radius_polar,
34    double radius_equatorial)
35{
36    // translate a point in [-1,1] to a point in [0, 1]
37    const double zm = 0.5;
38    const double zb = 0.5;
39    double total = 0.0;
40    for (int i=0;i<76;i++) {
41        //const double sin_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2;
42        const double sin_alpha = Gauss76Z[i]*zm + zb;
43        total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha);
44    }
45    // translate dx in [-1,1] to dx in [lower,upper]
46    const double form = total*zm;
47    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);
48    return 1.0e-4 * s * s * form;
49}
50
51double Iqxy(double qx, double qy,
52    double sld,
53    double sld_solvent,
54    double radius_polar,
55    double radius_equatorial,
56    double theta,
57    double phi)
58{
59    double q, sin_alpha, cos_alpha;
60    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
61    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha);
62    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);
63
64    return 1.0e-4 * form * s * s;
65}
66
Note: See TracBrowser for help on using the repository browser.