rho function

double rho(
  1. double phi
)

Distance from Earth center to point on ellipsoid at latitude phi.

Result is fraction of equatorial radius.

Implementation

double rho(double phi) {
  return 0.9983271 + 0.0016764 * math.cos(2 * phi) -
      0.0000035 * math.cos(4 * phi);
}