1 | static double |
---|
2 | form_volume(double rg) |
---|
3 | { |
---|
4 | return 1.0; |
---|
5 | } |
---|
6 | |
---|
7 | static double |
---|
8 | radius_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 | |
---|
23 | static double |
---|
24 | gauss_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 | |
---|
59 | static double |
---|
60 | Iq(double q, double i_zero, double rg) |
---|
61 | { |
---|
62 | return i_zero * gauss_coil(q*rg); |
---|
63 | } |
---|