[153f8f6] | 1 | static double |
---|
| 2 | form_volume(double rg) |
---|
[d277229] | 3 | { |
---|
| 4 | return 1.0; |
---|
| 5 | } |
---|
| 6 | |
---|
| 7 | static double |
---|
[a34b811] | 8 | radius_effective(int mode, double rg) |
---|
[d277229] | 9 | { |
---|
[ee60aa7] | 10 | switch (mode) { |
---|
[d42dd4a] | 11 | default: |
---|
[ee60aa7] | 12 | case 1: // R_g |
---|
[d277229] | 13 | return rg; |
---|
[ee60aa7] | 14 | case 2: // 2R_g |
---|
[d277229] | 15 | return 2.0*rg; |
---|
[ee60aa7] | 16 | case 3: // 3R_g |
---|
[d277229] | 17 | return 3.0*rg; |
---|
[ee60aa7] | 18 | case 4: // (5/3)^0.5*R_g |
---|
[d277229] | 19 | return sqrt(5.0/3.0)*rg; |
---|
| 20 | } |
---|
| 21 | } |
---|
| 22 | |
---|
[ee60aa7] | 23 | static double |
---|
| 24 | gauss_coil(double qr) |
---|
[d277229] | 25 | { |
---|
[ee60aa7] | 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; |
---|
[d277229] | 53 | } |
---|
[ee60aa7] | 54 | #endif |
---|
[d277229] | 55 | |
---|
[ee60aa7] | 56 | return 2.0 * (expm1(-x) + x)/(x*x); |
---|
| 57 | } |
---|
| 58 | |
---|
[153f8f6] | 59 | static double |
---|
| 60 | Iq(double q, double i_zero, double rg) |
---|
[ee60aa7] | 61 | { |
---|
| 62 | return i_zero * gauss_coil(q*rg); |
---|
[d277229] | 63 | } |
---|