source: sasmodels/sasmodels/models/mono_gauss_coil.c @ a34b811

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since a34b811 was a34b811, checked in by Paul Kienzle <pkienzle@…>, 5 months ago

use radius_effective/radius_effective_mode/radius_effective_modes consistently throughout the code

  • Property mode set to 100644
File size: 1.7 KB
Line 
1static double
2form_volume(double rg)
3{
4    return 1.0;
5}
6
7static double
8radius_effective(int mode, double rg)
9{
10    switch (mode) {
11    default:
12    case 1: // R_g
13        return rg;
14    case 2: // 2R_g
15        return 2.0*rg;
16    case 3: // 3R_g
17        return 3.0*rg;
18    case 4: // (5/3)^0.5*R_g
19        return sqrt(5.0/3.0)*rg;
20    }
21}
22
23static double
24gauss_coil(double qr)
25{
26    const double x = qr*qr;
27
28    // Use series expansion at low q for higher accuracy. We could use
29    // smaller polynomials if we sacrifice some digits of precision or
30    // introduce an additional series expansion around x == 1.
31    // See explore/precision.py, gauss_coil function.
32#if FLOAT_SIZE>4 // DOUBLE_PRECISION
33    // For double precision: use O(5) Pade with 0.5 cutoff (10 mad + 1 divide)
34    if (x < 0.5) {
35        // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
36        const double A1=1./12., A2=2./99., A3=1./2640., A4=1./23760., A5=-1./1995840.;
37        const double B1=5./12., B2=5./66., B3=1./132., B4=1./2376., B5=1./95040.;
38        return (((((A5*x + A4)*x + A3)*x + A2)*x + A1)*x + 1.)
39                / (((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.);
40    }
41#else
42    // For single precision: use O(7) Taylor with 0.8 cutoff (7 mad)
43    if (x < 0.8) {
44        const double C0 = +1.;
45        const double C1 = -1./3.;
46        const double C2 = +1./12.;
47        const double C3 = -1./60.;
48        const double C4 = +1./360.;
49        const double C5 = -1./2520.;
50        const double C6 = +1./20160.;
51        const double C7 = -1./181440.;
52        return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
53    }
54#endif
55
56    return 2.0 * (expm1(-x) + x)/(x*x);
57}
58
59static double
60Iq(double q, double i_zero, double rg)
61{
62    return i_zero * gauss_coil(q*rg);
63}
Note: See TracBrowser for help on using the repository browser.