docs: @Reference(authors = "T. Vincenty", //
        title = "Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations", //
        booktitle = "Survey Review 23:176", //
        url = "https://doi.org/10.1179/sre.1975.23.176.88", //
        bibkey = "doi:10.1179/sre.1975.23.176.88")
    public static double ellipsoidVincentyFormulaRad(double f, double lat1, double lon1, double lat2, double lon2) {
      final double dlon = (lon2 >= lon1) ? (lon2 - lon1) : (lon1 - lon2);
      final double onemf = 1 - f; // = 1 - (a-b)/a = b/a
  
      // Second eccentricity squared
      final double a_b = 1. / onemf; // = a/b
      final double ecc2 = (a_b + 1) * (a_b - 1); // (a^2-b^2)/(b^2)
  
      // Reduced latitudes:
      final double u1 = atan(onemf * tan(lat1));
      final double u2 = atan(onemf * tan(lat2));
      // Trigonometric values
      final DoubleWrapper tmp = new DoubleWrapper(); // To return cosine
      final double su1 = sinAndCos(u1, tmp), cu1 = tmp.value;
      final double su2 = sinAndCos(u2, tmp), cu2 = tmp.value;
  
      // Eqn (13) - initial value
      double lambda = dlon;
  
      for(int i = 0;; i++) {
        final double slon = sinAndCos(lambda, tmp), clon = tmp.value;
  
        // Eqn (14) - \sin \sigma
        final double term1 = cu2 * slon, term2 = cu1 * su2 - su1 * cu2 * clon;
        final double ssig = sqrt(term1 * term1 + term2 * term2);
        // Eqn (15) - \cos \sigma
        final double csig = su1 * su2 + cu1 * cu2 * clon;
        // Two identical points?
        if(!(ssig > 0)) {
          return 0.;
        }
        // Eqn (16) - \sigma from \tan \sigma
        final double sigma = atan2(ssig, csig);
        // Eqn (17) - \sin \alpha, and this way \cos^2 \alpha
        final double salp = cu1 * cu2 * slon / ssig;
        final double c2alp = (1. + salp) * (1. - salp);
        // Eqn (18) - \cos 2 \sigma_m
        final double ctwosigm = (Math.abs(c2alp) > 0) ? csig - 2.0 * su1 * su2 / c2alp : 0.;
        final double c2twosigm = ctwosigm * ctwosigm;
  
        // Eqn (10) - C
        final double cc = f * .0625 * c2alp * (4.0 + f * (4.0 - 3.0 * c2alp));
        // Eqn (11) - new \lambda
        final double prevlambda = lambda;
        lambda = dlon + (1.0 - cc) * f * salp * //
            (sigma + cc * ssig * (ctwosigm + cc * csig * (-1.0 + 2.0 * c2twosigm)));
        // Check for convergence:
        if(Math.abs(prevlambda - lambda) < PRECISION || i >= MAX_ITER) {
          // TODO: what is the proper result to return on MAX_ITER (antipodal
          // points)?
          // Definition of u^2, rewritten to use second eccentricity.
          final double usq = c2alp * ecc2;
          // Eqn (3) - A
          final double aa = 1.0 + usq / 16384.0 * (4096.0 + usq * (-768.0 + usq * (320.0 - 175.0 * usq)));
          // Eqn (4) - B
          final double bb = usq / 1024.0 * (256.0 + usq * (-128.0 + usq * (74.0 - 47.0 * usq)));
          // Eqn (6) - \Delta \sigma
          final double dsig = bb * ssig * (ctwosigm + .25 * bb * (csig * (-1.0 + 2.0 * c2twosigm) //
              - ONE_SIXTH * bb * ctwosigm * (-3.0 + 4.0 * ssig * ssig) * (-3.0 + 4.0 * c2twosigm)));
          // Eqn (19) - s
          return aa * (sigma - dsig);
        }
      }
    }