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

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since d42dd4a was d42dd4a, checked in by pkienzle, 6 years ago

fix compiler warnings for CUDA

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