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

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since ee60aa7 was ee60aa7, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

clean up effective radius functions; improve mono_gauss_coil accuracy; start moving VR into C

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