Function pj_rectifying_radius

src/latitudes.cpp:415–420

Compute the rectifying radius = quarter meridian / (pi/2 * a). The accuracy of this series is the same as those used for the computation of the auxiliary latitudes.

// of this series is the same as those used for the computation of the
// auxiliary latitudes.
double pj_rectifying_radius(double n) {
    // Expansion of (quarter meridian) / ((a+b)/2 * pi/2) as series in n^2;
    // these coefficients are ( (2*k - 3)!! / (2*k)!! )^2 for k = 0..3
    static const double coeff_rad[] = {1, 1.0 / 4, 1.0 / 64, 1.0 / 256};
    return pj_polyval(n * n, coeff_rad, 3) / (1 + n);
}