docs: @Reference(authors = "Erich Schubert, Arthur Zimek, Hans-Peter Kriegel", //
        title = "Geodetic Distance Queries on R-Trees for Indexing Geographic Data", //
        booktitle = "Int. Symp. Advances in Spatial and Temporal Databases (SSTD'2013)", //
        url = "https://doi.org/10.1007/978-3-642-40235-7_9", //
        bibkey = "DBLP:conf/ssd/SchubertZK13")
    public static double latlngMinDistRadFull(double plat, double plng, double rminlat, double rminlng, double rmaxlat, double rmaxlng) {
      // FIXME: add support for rectangles crossing the +-180 deg boundary!
  
      // Degenerate rectangles:
      if((rminlat >= rmaxlat) && (rminlng >= rmaxlng)) {
        return haversineFormulaRad(rminlat, rminlng, plat, plng);
      }
  
      // The simplest case is when the query point is in the same "slice":
      if(rminlng <= plng && plng <= rmaxlng) {
        return (rminlat <= plat && plat <= rmaxlat) ? 0 // Inside
            : (plat < rminlat) ? rminlat - plat : plat - rmaxlat; // S, N
      }
  
      // Determine whether going east or west is shorter.
      double lngE = rminlng - plng;
      lngE = lngE >= 0 ? lngE : lngE + TWOPI;
      double lngW = plng - rmaxlng; // we keep this negative!
      lngW = lngW >= 0 ? lngW : lngW + TWOPI;
  
      // Compute sine and cosine values we will certainly need below:
      final DoubleWrapper tmp = new DoubleWrapper(); // To return cosine
      final double slatQ = sinAndCos(plat, tmp), clatQ = tmp.value;
      final double slatN = sinAndCos(rmaxlat, tmp), clatN = tmp.value;
      final double slatS = sinAndCos(rminlat, tmp), clatS = tmp.value;
  
      // Head east, to min edge:
      if(lngE <= lngW) {
        final double slngD = sinAndCos(lngE, tmp), clngD = tmp.value;
  
        // Bearing to south
        // atan2(slngD * clatS, clatQ * slatS - slatQ * clatS * clngD);
        // Bearing from south
        final double bs = atan2(slngD * clatQ, clatS * slatQ - slatS * clatQ * clngD);
        // Bearing to north
        // atan2(slngD * clatN, clatQ * slatN - slatQ * clatN * clngD);
        // Bearing from north
        final double bn = atan2(slngD * clatQ, clatN * slatQ - slatN * clatQ * clngD);
        if(bs < HALFPI && bn > HALFPI) {
          // Radians from south pole = abs(ATD)
          final double radFromS = -HALFPI - plat;
  
          // Cross-track-distance to longitude line.
          return asin(sin(radFromS) * -slngD);
        }
        if(bs - HALFPI < HALFPI - bn) {
          // Haversine to north corner.
          final double slatN2 = sin((plat - rmaxlat) * .5);
          final double slon = sin(lngE * .5);
          final double aN = slatN2 * slatN2 + slon * slon * clatQ * clatN;
          return 2 * atan2(sqrt(aN), sqrt(1 - aN));
        }
        else {
          // Haversine to south corner.
          final double slatS2 = sin((plat - rminlat) * .5);
          final double slon = sin(lngE * .5);
          final double aS = slatS2 * slatS2 + slon * slon * clatQ * clatS;
          return 2 * atan2(sqrt(aS), sqrt(1 - aS));
        }
      }
      else { // Head west, to max edge
        final double slngD = -sinAndCos(lngW, tmp), clngD = tmp.value;
  
        // Bearing to south
        // atan2(slngD * clatS, clatQ * slatS - slatQ * clatS * clngD);
        // Bearing from south
        final double bs = atan2(slngD * clatQ, clatS * slatQ - slatS * clatQ * clngD);
        // Bearing to north
        // atan2(slngD * clatN, clatQ * slatN - slatQ * clatN * clngD);
        // Bearing from north
        final double bn = atan2(slngD * clatQ, clatN * slatQ - slatN * clatQ * clngD);
        if(bs > -HALFPI && bn < -HALFPI) {
          // Radians from south = abs(ATD) = distance from pole
          final double radFromS = -HALFPI - plat;
          // Cross-track-distance to longitude line.
          return asin(sin(radFromS) * slngD);
        }
        if(-HALFPI - bs < bn + HALFPI) {
          // Haversine to north corner.
          final double slatN2 = sin((plat - rmaxlat) * .5);
          final double slon = sin(lngW * .5);
          final double aN = slatN2 * slatN2 + slon * slon * clatQ * clatN;
          return 2 * atan2(sqrt(aN), sqrt(1 - aN));
        }
        else {
          // Haversine to south corner.
          final double slatS2 = sin((plat - rminlat) * .5);
          final double slon = sin(lngW * .5);
          final double aS = slatS2 * slatS2 + slon * slon * clatQ * clatS;
          return 2 * atan2(sqrt(aS), sqrt(1 - aS));
        }
      }
    }