source: sasmodels/sasmodels/models/ellipsoid.c @ 130d4c7

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 130d4c7 was 130d4c7, checked in by Paul Kienzle <pkienzle@…>, 3 years ago

restore broken ellipsoid 2-D calculation to 4.0.1 equivalent

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