Code: Select all
double g(const double s) { return s / std::tanh(s) - 1.0; }
double R(const double t, const double s) {
double s2 = s * s;
double s4 = s2 * s2;
if (s < 0.03) {
return (3072.0 + t * (384.0 + t * (24.0 + t))) / 3072.0 - t * (2688.0 + t * (80.0 + 21.0 * t)) / 322560.0 * s2 +
t * (2816.0 - t * (88.0 + 63.0 * t)) / 3548160.0 * s4;
}
double s6 = s2 * s4;
double t2 = t * t;
double t3 = t2 * t;
double gval = g(s);
return 1.0 + 3.0 * t * gval / (8.0 * s2) - (5.0 * t2 * (-8.0 * s2 + gval * (24.0 + 3.0 * gval))) / (128.0 * s4) +
(35.0 * t3 * (-40.0 * s2 + gval * (120 + gval * (24.0 + 3.0 * gval)))) / (1024.0 * s6);
}