package nxm.sys.test.junit.lib;

import static java.lang.Math.acos;
import static java.lang.Math.max;
import static java.lang.Math.min;
import static java.lang.Math.sqrt;
import static java.lang.Math.toDegrees;
import static nxm.sys.inc.Constants.EARTH_EQUATORIAL_RADIUS;
import static nxm.sys.inc.Constants.TWOPI;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

import nxm.sys.lib.GeodeticUtil;
import nxm.sys.lib.Midas;
import nxm.sys.lib.Position;
import nxm.sys.lib.Table;
import nxm.sys.lib.Transform;
import nxm.sys.test.junit.NMTestCase;

/** This test the GeodeticUtil class
    @author  M. Scott Meuleners
 */
public final class GeodeticUtilTest extends NMTestCase {

  // Test Polygons

  private static final double[] TEST_POLYGON_1  // wierd shape that overlaps 180/-180 boundary
          = {  0.0,    85.0,
               30.0,   85.0,
               60.0,   85.0,
               90.0,   85.0,
               120.0,  85.0,
               150.0,  85.0,
               155.0,  85.0,
               160.0,  85.0,
               160.0,  85.0,
               160.0,  75.0,
               160.0,  70.0,
               160.0,  65.0,
               160.0,  60.0,
               160.0,  55.0,
               160.0,  50.0,
               160.0,  45.0,
               160.0,  40.0,
               170.0,  35.0,
               175.0,  30.0,
              -175.0,  30.0,
              -170.0,  35.0,
              -165.0,  40.0,
              -170.0,  45.0,
              -175.0,  50.0,
               175.0,  50.0,
               170.0,  55.0,
               165.0,  60.0,
               170.0,  65.0,
               175.0,  70.0,
              -175.0,  70.0,
              -170.0,  75.0,
              -165.0,  80.0,
              -160.0,  85.0,
              -155.0,  85.0,
              -120.0,  85.0,
              -90.0,   85.0,
              -60.0,   85.0,
              -30.0,   85.0 };

  private static final double[][] TEST_POINTS_1
          = { { 0.0,  90.0   }, // inside (north pole)
              { 0.0, -90.0   }, // outside (south pole)
              { 0.0,    0.0  }, // outside
              { 177.0,  87.0 }, // inside
              {-177.0,  87.0 }, // inside
              { 162.0,  60.0 }, // inside
              { 170.0,  60.0 }, // outside
              { 177.0,  60.0 }, // outside
              {-177.0,  60.0 }, // outside
              {-170.0,  60.0 }, // outside
              { 170.0,  40.0 }, // inside
              { 177.0,  40.0 }, // inside
              {-177.0,  40.0 }, // inside
              {-170.0,  40.0 }, // inside
              { 175.0,  75.0 }};// inside

  private static final boolean[] EXPECTED_1
          = { true, false, false, true, true, true, false, false, false, false, true, true, true, true, true};

  private static final double[] TEST_POLYGON_2 // circle around north pole
          = {  0.0,     10.0,
               30.0,    10.0,
               60.0,    10.0,
               90.0,    10.0,
               120.0,   10.0,
               150.0,   10.0,
               175.0,   10.0,
              -150.0,   10.0,
              -120.0,   10.0,
              -90.0,    10.0,
              -60.0,    10.0,
              -30.0,    10.0,
              -5.0,     10.0};

  private static final double[][] TEST_POINTS_2
          = { { 0.0,    90.0 }, // inside (north pole)
              { 0.0,   -90.0 }, // outside (south pole)
              { 0.0,   -20.0 }, // outside
              { 177.0, -30.0 }, // outside
              { 170.0, -30.0 }, // outside
              {-170.0, -30.0 }, // outside
              {-140.0, -30.0 }, // outside
              { 0.0,    20.0 }, // inside
              { 177.0,  30.0 }, // inside
              { 170.0,  30.0 }, // inside
              {-170.0,  30.0 }, // inside
              {-140.0,  30.0 }, // inside
              { 120.0, -20.0 }};// outside

  private static final boolean[] EXPECTED_2
          = { true, false, false, false, false, false, false, true, true, true, true, true, false};


 private static final double[] TEST_POLYGON_3 // circle around south pole
          = {  0.0,    -10.0,
               30.0,   -10.0,
               60.0,   -10.0,
               90.0,   -10.0,
               120.0,  -10.0,
               150.0,  -10.0,
               175.0,  -10.0,
              -150.0,  -10.0,
              -120.0,  -10.0,
              -90.0,   -10.0,
              -60.0,   -10.0,
              -30.0,   -10.0,
              -5.0,    -10.0};

  private static final double[][] TEST_POINTS_3
          = { { 0.0,    90.0 },// outside (north pole)
              { 0.0,   -90.0 },// inside (south pole)
              { 20.0,  -20.0 },// inside
              { 177.0, -30.0 },// inside
              { 170.0, -30.0 },// inside
              {-170.0, -30.0 },// inside
              {-140.0, -30.0 },// inside
              { 0.0,    20.0 },// outside
              { 177.0,  30.0 },// outside
              { 170.0,  30.0 },// outside
              {-170.0,  30.0 },// outside
              {-140.0,  30.0 },// outside
              { 120.0, -20.0 }};// inside

  private static final boolean[] EXPECTED_3
          = { false, true, true, true, true, true, true, false, false, false, false, false, true};

  private static final double[][] TEST_POINTS_4 // ellipse centered at 10.0,10.0
          = { { 0.0,    90.0 }, // outside (north pole)
              { 0.0,   -90.0 }, // outside (south pole)
              { 10.0,   10.0 }, // inside (center of ellipse
              { 25.0,   25.0 }};// outside

  private static final boolean[] EXPECTED_4
          = { false, false, true, false };

  private static final double[][] TEST_POINTS_5 // Donut centered at 170, 85
          = { { 0.0,    90.0 }, // outside (north pole)
              { 0.0,   -90.0 }, // outside (south pole)
              { 170.0,  85.0 }, // outside (center of donut)
              { 170.0,  70.0 }, // inside
              {-10.0,   80.0 }};// inside

  private static final boolean[] EXPECTED_5
          = { false, false, false, true, true };


  private static final double[] DATE_BOUNDARY_NP_1 = {
      +95, +77,
     +160,  +75,
     +170, +80,
      -81, 85,
      -22, 80,
      +47, 80,
      +95, +77
  };

  private static final double[] DATE_BOUNDARY_NP_2 = {
      +75, +77,
     +140,  +75,
     +150, +80,
     -101, 85,
      -42, 80,
      +27, 80,
      +75, +77
  };

  private static final double[] DATELINE_AT_LAT0_1 = {
     +170, +20,
      -85, +20,
      -85, -20,
     +170, -20,
     +170, +20
  };

  private static final double[] DATELINE_AT_LAT0_2 = {
      +150, +20,
      -105, +20,
      -105, -20,
      +150, -20,
      +150, +20
   };

  /** public constructor to allow running this test case directly with JUnit. */
  public GeodeticUtilTest (String methodName) { this(methodName, null); }

  /** Creates a new instance of the class.
      @param methodName The name of the method being called.
      @param midas      The midas context testTo use when running the tests.
   */
  public GeodeticUtilTest (String methodName, Midas midas) {
    super(methodName, midas);
  }

  //////////////////////////////////////////////////////////////////////////////////////////////////

  public void testInsidePolygon() {
    // Adding test for new overloaded insidePolygon added in 4.0.0 for NXM-3727
    Position[] boxArray = new Position[] {Position.fromGeo(0, 0, 0),
        Position.fromGeo(0, 4, 0),
        Position.fromGeo(0, 4, 2),
        Position.fromGeo(0, 0, 2),
        Position.fromGeo(0, 0, 0)};
    List<Position> box = Arrays.asList(boxArray);

    for (double lon = -180; lon <= +180; lon++) {
      for (double lat = -90; lat <= +90; lat++) {
        Position p  = Position.fromGeo(0, lat, lon);
        boolean in = (0 <= lat) && (lat <= 4) &&
            (0 <= lon) && (lon <= 2);
        assertEquals("insidePolygon("+p+", box)", GeodeticUtil.insidePolygon(p,box), in);
      }
    }

    // Original insidePolygon testing
    Position p = new Position();

    // Test arbitrary shape, crossing -180/180 longitude and encircling north pole
    for (int i = 0; i<TEST_POINTS_1.length; i++) {
      p.setGeo(0.0, TEST_POINTS_1[i][1], TEST_POINTS_1[i][0]);
      assertEquals("GeodeticUtil.insidePolygon(p, 2, TEST_POLYGON_1) i=" + i, EXPECTED_1[i], GeodeticUtil.insidePolygon(p, 2, TEST_POLYGON_1));
    }

    //Test circle around north pole at 10 degrees latitude
    for (int i = 0; i<TEST_POINTS_2.length; i++) {
      p.setGeo(0.0, TEST_POINTS_2[i][1], TEST_POINTS_2[i][0]);
      assertEquals("GeodeticUtil.insidePolygon(p, 2, TEST_POLYGON_2) i=" + i, EXPECTED_2[i], GeodeticUtil.insidePolygon(p, 2, TEST_POLYGON_2));
    }

    // Test circle around south pole at -10 degrees latitude
    for (int i = 0; i<TEST_POINTS_3.length; i++) {
      p.setGeo(0.0, TEST_POINTS_3[i][1], TEST_POINTS_3[i][0]);
      assertEquals("GeodeticUtil.insidePolygon(p, 2, TEST_POLYGON_3) i=" + i, EXPECTED_3[i], GeodeticUtil.insidePolygon(p, 2, TEST_POLYGON_3));
    }

    // Test ellipse centered at 10.0, 10.0, with 10000 semi-major, 5000 semi-minor and 45 degree tilt
    double ellipse[] = GeodeticUtil.computeEllipseLL(10.0, 10.0, 10000.0, 5000.0, 45, 60);
    for (int i = 0; i<TEST_POINTS_4.length; i++) {
      p.setGeo(0.0, TEST_POINTS_4[i][1], TEST_POINTS_4[i][0]);
      assertEquals("GeodeticUtil.insidePolygon(p, 2, ellipse) i=" + i, EXPECTED_4[i], GeodeticUtil.insidePolygon(p, 2, ellipse));
    }

    // Test donut centered at 170.0, 85.0 with outer diameter 20 deg latitude and inner diameter 10 deg latitude
    double donutInner[] = GeodeticUtil.computeCircleLL(170.0, 85.0, 10, 60);
    double donutOuter[] = GeodeticUtil.computeCircleLL(170.0, 85.0, 20, 60);
    for (int i = 0; i<TEST_POINTS_5.length; i++) {
      p.setGeo(0.0, TEST_POINTS_5[i][1], TEST_POINTS_5[i][0]);
      assertEquals("GeodeticUtil.insidePolygon(p, 2, donutOuter, donutInner) i=" + i, EXPECTED_5[i], GeodeticUtil.insidePolygon(p, 2, donutOuter, donutInner));
    }

    insidePolygonDateline();
    insidePolygonEdges();
  }

  private List<Position> positionListFromDoubleArray (double[] coords, int dim) {
    List<Position> positions = new ArrayList<>();
    boolean haveAlt = dim == 3;
    for (int i = 0; i < coords.length; i += dim) {
      Position p = new Position();
      double alt = (haveAlt) ? coords[i+2] : 0;
      p.setGeo(alt, coords[i+1], coords[i]);
      positions.add(p);
    }
    return positions;
  }

  private void insidePolygonDateline () {
    List<Position> dateBoundaryNP1List = positionListFromDoubleArray(DATE_BOUNDARY_NP_1,2);
    List<Position> dateBoundaryNP2List = positionListFromDoubleArray(DATE_BOUNDARY_NP_2,2);
    List<Position> datelineAtLat01 = positionListFromDoubleArray(DATELINE_AT_LAT0_1,2);
    List<Position> datelineAtLat02 = positionListFromDoubleArray(DATELINE_AT_LAT0_2,2);

    Position NORTH_POLE = Position.fromGeo(0, +90, 0);
    Position DATELINE_EQUATOR = Position.fromGeo(0, 0, +180);
    Position CENTER_EQUATOR = Position.fromGeo(0, 0, 0);
    assertEquals("GeodeticUtil.insidePolygon North Pole 1",true,GeodeticUtil.insidePolygon(NORTH_POLE, 2, DATE_BOUNDARY_NP_1));
    assertEquals("GeodeticUtil.insidePolygon North Pole 2",true,GeodeticUtil.insidePolygon(NORTH_POLE, 2, DATE_BOUNDARY_NP_2));
    assertEquals("GeodeticUtil.insidePolygon Lon=+180 at Equator 1",true,GeodeticUtil.insidePolygon(DATELINE_EQUATOR, 2, DATELINE_AT_LAT0_1));
    assertEquals("GeodeticUtil.insidePolygon Lon=+180 at Equator 2",true,GeodeticUtil.insidePolygon(DATELINE_EQUATOR, 2, DATELINE_AT_LAT0_2));
    assertEquals("GeodeticUtil.insidePolygon Lon=+180 at Equator 1 (0,0,0)",false,GeodeticUtil.insidePolygon(CENTER_EQUATOR, 2, DATELINE_AT_LAT0_1));
    // Jennifer's new version
    assertEquals("GeodeticUtil.insidePolygon IUO North Pole 1",true,GeodeticUtil.insidePolygon(NORTH_POLE, dateBoundaryNP1List));
    assertEquals("GeodeticUtil.insidePolygon IUO North Pole 2",true,GeodeticUtil.insidePolygon(NORTH_POLE, dateBoundaryNP2List));
    assertEquals("GeodeticUtil.insidePolygon IUO Lon=+180 at Equator 1",true,GeodeticUtil.insidePolygon(DATELINE_EQUATOR, datelineAtLat01));
    assertEquals("GeodeticUtil.insidePolygon IUO Lon=+180 at Equator 2",true,GeodeticUtil.insidePolygon(DATELINE_EQUATOR, datelineAtLat02));

    GeodeticUtil.setinsidePolyFix(false);
    assertEquals("GeodeticUtil.insidePolygon North Pole 1 bad code",false,GeodeticUtil.insidePolygon(NORTH_POLE, 2, DATE_BOUNDARY_NP_1));
    assertEquals("GeodeticUtil.insidePolygon Lon=+180 at Equator 1 bad code",false,GeodeticUtil.insidePolygon(DATELINE_EQUATOR, 2, DATELINE_AT_LAT0_1));
    GeodeticUtil.setinsidePolyFix(true);
  }

  private void insidePolygonEdges () {
    double[] polygonPoints = DATELINE_AT_LAT0_1;
    List<Position> polygonPositions = positionListFromDoubleArray(DATELINE_AT_LAT0_1,2);
    Position horizontalEdgeCheck = Position.fromGeo(0,-20,-140);
    assertEquals("GeodeticUtil.insidePolygon horizontal edge",true,GeodeticUtil.insidePolygon(horizontalEdgeCheck,2,polygonPoints));
    assertEquals("GeodeticUtil.insidePolygon IUO horizontal edge",true,GeodeticUtil.insidePolygon(horizontalEdgeCheck,polygonPositions));
    Position verticalEdgeCheck = Position.fromGeo(0,0,170);
    assertEquals("GeodeticUtil.insidePolygon vertical edge",true,GeodeticUtil.insidePolygon(verticalEdgeCheck,2,polygonPoints));
    assertEquals("GeodeticUtil.insidePolygon IUO vertical edge",true,GeodeticUtil.insidePolygon(verticalEdgeCheck,polygonPositions));
    Position vertexEdgeCheck = Position.fromGeo(0,20,170);
    assertEquals("GeodeticUtil.insidePolygon vertex",true,GeodeticUtil.insidePolygon(vertexEdgeCheck,2,polygonPoints));
    assertEquals("GeodeticUtil.insidePolygon IUO vertex",true,GeodeticUtil.insidePolygon(vertexEdgeCheck,polygonPositions));

    polygonPoints = new double[]{
        0,  0,
        0, 90,
      180,  0,
      -90,  0,
        0,  0
     };
    polygonPositions = positionListFromDoubleArray(polygonPoints, 2);
    Position northPoleEdgeCheck = Position.fromGeo(0,90,0);
    assertEquals("GeodeticUtil.insidePolygon vertex at the North Pole",true,GeodeticUtil.insidePolygon(northPoleEdgeCheck,2,polygonPoints));
    assertEquals("GeodeticUtil.insidePolygon IUO vertex at the North Pole",true,GeodeticUtil.insidePolygon(northPoleEdgeCheck,polygonPositions));
  }

  public void testComputeCircleLL() {

    //                               Lon    Lat   radius  Distance
    double[][] testComputeCircleLL = { {  0.0, 90.0,  30.0 },
                                       {  0.0, 0.0,   30.0 },
                                       {-20.0,-90.0,  40.0 },
                                       {-20.0, 30.0,  15.0 },
                                       {  0.0, 90.0,  90.0 },
                                       { 540.0, 50.0, 20.0 }};// added for normalized lon testing since 3.9.0
    for (int i = 0; i<6; i++) {
      Position p = Position.fromGeo(0.0,testComputeCircleLL[i][1],testComputeCircleLL[i][0]);
      double radius = testComputeCircleLL[i][2];
      double[] circle = GeodeticUtil.computeCircleLL(p, radius, 60);
      for (int j = 0; j<60; j+=10) {
        Position q = Position.fromGeo(0.0,circle[2*j+1],circle[2*j+0]);
        double distDegLat = Transform.greatArc(p, q) / (TWOPI * EARTH_EQUATORIAL_RADIUS) * 360.0;
        assertEquals("circle[" + j + "]", radius, distDegLat, radius * 0.01);
      }
    }

  }

  public void testComputeConeLLA() {
    //                               lonP latP altP lonQ latQ altQ angle
    double[][] testComputeConeLLA = {{  0,   0, 4E7,   0,  80,   0,    5},
                                     {  0,   0, 4E7,   0,  80,   0,   10},
                                     {  0,   0, 4E7,   0,  80,   0,   15},
                                     {  0,  90, 4E7,   0,  15,   0,   15},
                                     { 40,  40, 4E7,   5,   5,   0,   25}};

    Transform.setUseLegacyConeCode(true);
    for (int i = 0; i < testComputeConeLLA.length; i++) {
      Position p = Position.fromGeo(testComputeConeLLA[i][2],testComputeConeLLA[i][1],testComputeConeLLA[i][0]);
      Position q = Position.fromGeo(testComputeConeLLA[i][5],testComputeConeLLA[i][4],testComputeConeLLA[i][3]);
      double[] cone = GeodeticUtil.computeConeLLA(p, q, testComputeConeLLA[i][6], 60, 0, 1);
      for (int j = 0; j<60; j+=4) {
        Position r = Position.fromGeo(cone[3*j+2],cone[3*j+1],cone[3*j+0]);
        p.geo2car();
        q.geo2car();
        r.geo2car();
        // check angle between line pq and line pr, should be equal to angle above
        double pqx = (p.x-q.x);
        double pqy = (p.y-q.y);
        double pqz = (p.z-q.z);
        double pqr = sqrt((pqx)*(pqx) + (pqy)*(pqy) + (pqz)*(pqz));
        double prx = (p.x-r.x);
        double pry = (p.y-r.y);
        double prz = (p.z-r.z);
        double prr = sqrt((prx)*(prx) + (pry)*(pry) + (prz)*(prz));
        double rqx = (r.x-q.x);
        double rqy = (r.y-q.y);
        double rqz = (r.z-q.z);
        double rqr = sqrt((rqx)*(rqx) + (rqy)*(rqy) + (rqz)*(rqz));
        double angle = toDegrees(acos(-(rqr*rqr-pqr*pqr-prr*prr) / (2 * pqr * prr)));  // law of cosines for triangles

        assertEquals("angle j=" + j, testComputeConeLLA[i][6], angle, 0.1);
      }
    }

    /* Test new Transform.clipAtEarth code against original Transform.clipAtEarth code
     * Scenario #1: bad cone (legacy code) does not intersect with the Earth
     * Scenario #2: bad cone (legacy code) partially intersects with the Earth
     * NOTE: not testing case where bad cone (legacy code) completely intersects with the Earth since bad cone (legacy code)
     *       shows no discrepancies when all points in the cone intersect with the Earth
     */
    //                                   lonP  latP        altP lonQ  latQ    altQ   angle lonTestPt latTestPt altTestPt
    double[][] testComputeConeLLAClip = {{ 20,   15,      3.7E7,  20,   50,      0,  14.75,      120,     61.5,        0},  // Scenario #1
                                         {  0,   90,        3E7,   0,   90,      0,     14,        0,    10.09,        0},
                                         {  0,  -90,        3E7,   0,  -90,      0,     14,        0,   -10.09,        0},
                                         {  0,    0,        3E7,   0,    0,      0,     14,    69.36,       60,        0},
                                         {  0,    0,        3E7,  60,    0,      0,     22,     69.4,       60,        0},
                                         {  0,    0,        3E7,  45,   45,      0,     22,     69.4,       60,        0},
                                         { 20,  -20,        3E7,   0,   30,      0,     22,        0,     58.2,        0},
                                         {  0,   45,        3E6,   0,   90,      0,      4,   133.45,       87,        0},  // Scenario #2
                                         {  0,   45,        3E6,   0,   90,      0,     14,      -90,    75.08,        0},
                                         {  0,   60,        3E7,   0,   90,      0,     14,     -150,       37,        0},
                                         {  0,  -45,        3E6,   0,  -90,      0,      4,   133.45,      -87,        0},
                                         {  0,  -45,        3E6,   0,  -90,      0,     14,   133.45,   -86.93,        0},
                                         {  0,  -60,        3E7,   0,  -90,      0,     14,     -150,      -37,        0},
                                         {  0,    0,        3E7,  60,    0,      0,     10,       75,    47.36,        0},
                                         {  0,    0,        3E7,  45,   45,      0,     10,    69.45,       60,        0},
                                         { 20,  -20,        3E7,   0,   30,      0,     10,        0,     58.2,        0},
                                         {  0,  -30,      3.7E7,   0,    0,      0,     13,    124.1,      -60,        0},
                                         { 16, -3.9, 3.481502E7,  67, 34.3, 0.1997, 0.6332,       90,    39.44,        0}};

    // Whether or not the good cone (default code) should be inside the bad cone (legacy code)
    boolean[] defaultInsideLegacy = {true,true,true,false,false,false,false, // Scenario #1
        true,true,true,true,true,true,false,false,false,true,false};         // Scenario #2

    boolean runMetrics = false; // Toggle metrics on or off

    /* Based on whether or not the good cone (default code) should be inside the bad cone (legacy code),
     * Tests old code and new code by checking whether or not a given test point for each scenario is either:
     *   1. Outside the resulting good cone (default code) and inside the resulting bad cone (legacy code)
     *   2. Outside the resulting bad cone (legacy code) and inside the resulting good cone (default code)
     */
    for (int i = 0; i < testComputeConeLLAClip.length; i++) {
      if (runMetrics) System.out.println(i);
      Position p = Position.fromGeo(testComputeConeLLAClip[i][2],testComputeConeLLAClip[i][1],testComputeConeLLAClip[i][0]); // POS
      Position q = Position.fromGeo(testComputeConeLLAClip[i][5],testComputeConeLLAClip[i][4],testComputeConeLLAClip[i][3]); // POSP

      // Test point to check in between cones
      Position testPoint = Position.fromGeo(testComputeConeLLAClip[i][9],testComputeConeLLAClip[i][8],testComputeConeLLAClip[i][7]);

      /* If the good cone (default code) is inside the bad cone (legacy code):
       *   check that the test point is outside the good cone (default code) for the current scenario
       * If the bad cone (legacy code) is inside the good cone (default code):
       *   check that the test point is outside the bad cone (legacy code) for the current scenario
       */
      Transform.setUseLegacyConeCode(!defaultInsideLegacy[i]);
      double[] cone = GeodeticUtil.computeConeLLA(p, q, testComputeConeLLAClip[i][6], 60, 0, 1);
      assertFalse(GeodeticUtil.insidePolygon(testPoint,3,cone));
      int metricsLoopCount = 10000;
      if (runMetrics) runMetrics(() -> GeodeticUtil.computeConeLLA(p, q, 14.75, 60, 0, 1), "GeodeticUtil.computeConeLLA new code", metricsLoopCount);

      /* If the good cone (default code) is inside the bad cone (legacy code):
       *   check that the test point is inside the bad cone (legacy code) for the current scenario
       * If the bad cone (legacy code) is inside the good cone (default code):
       *   check that the test point is inside the good cone (default code) for the current scenario
       */
      Transform.setUseLegacyConeCode(defaultInsideLegacy[i]);
      cone = GeodeticUtil.computeConeLLA(p, q, testComputeConeLLAClip[i][6], 60, 0, 1);
      assertTrue(GeodeticUtil.insidePolygon(testPoint,3,cone));
      if (runMetrics) runMetrics(() -> GeodeticUtil.computeConeLLA(p, q, 14.75, 60, 0, 1), "GeodeticUtil.computeConeLLA original code", metricsLoopCount);
    }
    Transform.setUseLegacyConeCode(false);
  }

  /**
   * Note: Expected from Don Cross calculations (independent source, not necessarily more accurate)
   * @since 3.6.0 for Bug 2874
   */
  public void testComputeCourse() {
                                //  lon1 lat1   alt1        lon2  lat2   alt2        azimuth
    double[][] testPositionPairs = {{ 0.0,  0.0, 40000000.0, 40.0,  40.0, 40000000.0,  37.4537},
                                    { 0.0,  0.0, 40000000.0, 80.0,   0.0, 40000000.0,  90.0},
                                    {90.0,  0.0, 40000000.0, 15.0,   0.0, 40000000.0, 270.0},
                                    {40.0, 40.0, 40000000.0,  5.0,   5.0, 40000000.0, 231.30}};
    for (int i=0; i < testPositionPairs.length; i++) {
      Position p1 = Position.fromGeo(testPositionPairs[i][2],testPositionPairs[i][1],testPositionPairs[i][0]);
      Position p2 = Position.fromGeo(testPositionPairs[i][5],testPositionPairs[i][4],testPositionPairs[i][3]);
      double azimuth = GeodeticUtil.computeCourse(p1, p2);
      assertEquals("azimuth from position pair " + i, testPositionPairs[i][6], azimuth, 0.0001);
    }
  }

  /**
   * Note: Expected from Don Cross calculations (independent source, not necessarily more accurate)
   * @since 3.6.0 for Bug 2874
   */
  public void testComputeElevation() {
                                 //  lon1 lat1   alt1        lon2  lat2   alt2        elevation angle
    double[][] testPositionPairs = {{ 0.0,  0.0, 40000000.0, 40.0,  40.0, 40000000.0, -27.0368},
                                    { 0.0,  0.0, 40000000.0, 80.0,   0.0, 40000000.0, -40.0},
                                    {90.0,  0.0, 40000000.0, 15.0,   0.0, 40000000.0, -37.5},
                                    {40.0, 40.0, 40000000.0,  5.0,   5.0, 40000000.0, -23.5312}};
    for (int i=0; i < testPositionPairs.length; i++) {
      Position p1 = Position.fromGeo(testPositionPairs[i][2],testPositionPairs[i][1],testPositionPairs[i][0]);
      Position p2 = Position.fromGeo(testPositionPairs[i][5],testPositionPairs[i][4],testPositionPairs[i][3]);
      double elevationAngle = GeodeticUtil.computeElevation(p1, p2, 0);
      assertEquals("elevation angle from position pair " + i, testPositionPairs[i][6], elevationAngle, 0.01);
    }
  }

  public void testComputeEllipseLL() {
//                                    lon  lat   semi-major  semi-minor tilt
    double[][] testComputeConeLLA = {{0.0, 0.0,  2000000.0,    1000000.0,   45.0},
                                     {0.0, 0.0,  20000.0,      10000.0,     0.0},
                                     {0.0, 0.0,  2000.0,       1000.0,      90.0},
                                     {90.0,90.0, 2000000.0,    1000000.0,   45.0},
                                     {90.0,90.0, 20000.0,      10000.0,     0.0},
                                     {90.0,90.0, 2000.0,       1000.0,      90.0}};


    for (int i = 0; i<testComputeConeLLA.length; i++) {
      double[] ellipse = GeodeticUtil.computeEllipseLL(testComputeConeLLA[i][0], testComputeConeLLA[i][1],
                                                       testComputeConeLLA[i][2], testComputeConeLLA[i][3],
                                                       testComputeConeLLA[i][4], 60);


      double maxAngle = 0.0;
      double minAngle = 180.0;
      Position p = Position.fromGeo(testComputeConeLLA[i][2],testComputeConeLLA[i][1],testComputeConeLLA[i][0]);
      for (int j = 0; j<60; j++) {
        Position q = Position.fromGeo(0.0,ellipse[2*j+1],ellipse[2*j+0]);

        // check angle between line pq and line pr, should be equal to angle above
        double pqx = (p.x-q.x);
        double pqy = (p.y-q.y);
        double pqz = (p.z-q.z);
        double pqr = sqrt((pqx)*(pqx) + (pqy)*(pqy) + (pqz)*(pqz));
        double prx = (p.x-0.0);
        double pry = (p.y-0.0);
        double prz = (p.z-0.0);
        double prr = sqrt((prx)*(prx) + (pry)*(pry) + (prz)*(prz));
        double rqx = (0.0-q.x);
        double rqy = (0.0-q.y);
        double rqz = (0.0-q.z);
        double rqr = sqrt((rqx)*(rqx) + (rqy)*(rqy) + (rqz)*(rqz));
        double angle = toDegrees(acos(-(pqr*pqr-rqr*rqr-prr*prr) / (2 * rqr * prr)));
        maxAngle = max(angle,maxAngle);
        minAngle = min(angle,minAngle);
      }
      double max = testComputeConeLLA[i][2] / (TWOPI * EARTH_EQUATORIAL_RADIUS) * 360.0;
      double min = testComputeConeLLA[i][3] / (TWOPI * EARTH_EQUATORIAL_RADIUS) * 360.0;

      //System.out.println("max: " + max + " maxAngle: " + maxAngle);
      //System.out.println("min: " + min + " minAngle: " + minAngle);
      assertEquals("maxAngle i=" + i, max, maxAngle, 1);
      assertEquals("minAngle i=" + i, min, minAngle, 1);
    }
  }

  public void testComputeHorizonDistance() {
    double[] testComputeHorizonDistance = { 10.0, 500.0, 10000.0, 40000.0  };
    double[] expected                   = {  0.10,  0.72,    3.21,    6.40 };

    for (int i=0; i<4; i++) {
      double dist = GeodeticUtil.computeHorizonDistance(testComputeHorizonDistance[i]);

      assertEquals("dist i=" + i, expected[i], dist, 0.01);
    }

  }

  public void testComputeHorizonLL() {} // covered in testComputeCircleLL() / testComputeHorizonDistance()

  /**
   * computeMapEllipse() is the same as computeEllipseLL() except tilt is adjusted before
   * calling ComputeMapEllipse() with "90-tilt"
   * @since 3.6.0 for Bug 2874
   */
  public void testComputeMapEllipse() {
                                  //  lon  lat   semi-major  semi-minor tilt
    double[][] testComputeConeLLA = {{0.0, 0.0,  2000000.0,    1000000.0,   45.0},
                                     {0.0, 0.0,  20000.0,      10000.0,     90.0},
                                     {0.0, 0.0,  2000.0,       1000.0,       0.0},
                                     {90.0,90.0, 2000000.0,    1000000.0,   45.0},
                                     {90.0,90.0, 20000.0,      10000.0,     90.0},
                                     {90.0,90.0, 2000.0,       1000.0,      45.0}};


    for (int i = 0; i<testComputeConeLLA.length; i++) {
      double[] ellipse = GeodeticUtil.computeMapEllipse(testComputeConeLLA[i][0], testComputeConeLLA[i][1],
                         testComputeConeLLA[i][2], testComputeConeLLA[i][3],
                         testComputeConeLLA[i][4], 60);


      double maxAngle = 0.0;
      double minAngle = 180.0;
      Position p = Position.fromGeo(testComputeConeLLA[i][2],testComputeConeLLA[i][1],testComputeConeLLA[i][0]);
      for (int j = 0; j<60; j++) {
        Position q = Position.fromGeo(0.0,ellipse[2*j+1],ellipse[2*j+0]);

        // check angle between line pq and line pr, should be equal to angle above
        double pqx = (p.x-q.x);
        double pqy = (p.y-q.y);
        double pqz = (p.z-q.z);
        double pqr = sqrt((pqx)*(pqx) + (pqy)*(pqy) + (pqz)*(pqz));
        double prx = (p.x-0.0);
        double pry = (p.y-0.0);
        double prz = (p.z-0.0);
        double prr = sqrt((prx)*(prx) + (pry)*(pry) + (prz)*(prz));
        double rqx = (0.0-q.x);
        double rqy = (0.0-q.y);
        double rqz = (0.0-q.z);
        double rqr = sqrt((rqx)*(rqx) + (rqy)*(rqy) + (rqz)*(rqz));
        double angle = toDegrees(acos(-(pqr*pqr-rqr*rqr-prr*prr) / (2 * rqr * prr)));
        maxAngle = max(angle,maxAngle);
        minAngle = min(angle,minAngle);
      }
      double max = testComputeConeLLA[i][2] / (TWOPI * EARTH_EQUATORIAL_RADIUS) * 360.0;
      double min = testComputeConeLLA[i][3] / (TWOPI * EARTH_EQUATORIAL_RADIUS) * 360.0;

      //System.out.println("max: " + max + " maxAngle: " + maxAngle);
      //System.out.println("min: " + min + " minAngle: " + minAngle);
      assertEquals("maxAngle i=" + i, max, maxAngle, 1);
      assertEquals("minAngle i=" + i, min, minAngle, 1);
    }
  }

  public void testComputePosition() {
    Position p = new Position();
    Position q = new Position();

    //                               Lon    Lat   Heading  Distance
    double[][] testComputePosition = { {-20.0, 30.0,   0.0,    2500000.0 },
                                       {-20.0,  0.0,  45.0,    2500000.0 },
                                       {-20.0, 30.0,  90.0,          0.0 },
                                       {-20.0, 30.0, 235.0,    2500000.0 },
                                       {  0.0, 90.0,  45.0,    2500000.0 },
                                       { 30.0, 90.0,   0.0,    2500000.0 },
                                       { 30.0, 85.0,   0.0,   10000000.0 }};

    double[][] expected            = {{-20.0,  52.46 },
                                     {  -3.71, 15.67 },
                                     { -20.0,  30.0  },
                                     { -38.98, 15.80 },
                                     {  45.00, 67.54 },
                                     {   0.0,  67.54 },
                                     {-150.0,   5.17 }};


    for (int i = 0; i<7; i++) {
      p.setGeo(0.0, testComputePosition[i][1], testComputePosition[i][0]);
      q = GeodeticUtil.computePosition(p, testComputePosition[i][2], testComputePosition[i][3]);

      assertEquals("q.getLat() i=" + i, expected[i][1], q.getLat(), 0.01);
      assertEquals("q.getLon() i=" + i, expected[i][0], q.getLon(), 0.01);
    }

  }

  public void testFixAltitudeLLA() {

    double[] testFixAltitudeLLA = {  12.0, 12.0,      0.0005,
                                    -12.0, 12.0,      0.005,
                                    -12.0,-12.0,      0.000499,
                                     12.0,-12.0,      0.00005,
                                     12.0, 12.0,      0.5,
                                     12.0, 12.0,     50.0,
                                     12.0, 12.0,  20000.0,
                                     12.0, 12.0,     -0.0005,
                                    -12.0, 12.0,     -0.005,
                                    -12.0,-12.0,     -0.000499,
                                     12.0,-12.0,     -0.00005,
                                     12.0, 12.0,     -0.5,
                                     12.0, 12.0,    -50.0,
                                     12.0, 12.0, -20000.0       };

    double[]   expected =  { 0.0005, 0.005, 0.0, 0.0, 0.5, 50.0, 20000.0,
                             0.0,    0.0,   0.0, 0.0, 0.0,  0.0,     0.0 };

    double[] testFixAltitudeLLAReturn = GeodeticUtil.fixAltitudeLLA(testFixAltitudeLLA);

    // test both returned array and array passed in (changes are made in-place and returned).
    for (int i = 0; i < expected.length; i++) {
      assertEquals("testFixAltitudeLLA[" + (i*3+2) + "]", expected[i], testFixAltitudeLLA[i*3+2]);
    }

    for (int i = 0; i < expected.length; i++) {
      assertEquals("testFixAltitudeLLAReturn[" + (i*3+2) + "]", expected[i], testFixAltitudeLLAReturn[i*3+2]);
    }
  }

  public void testFromLLAToLL() {
    double[] testFromLLAToLL = { 12.0, 12.0, 25.0,
                                -12.0, 12.0, 25.0,
                                -12.0,-12.0, 25.0,
                                 12.0,-12.0, 25.0,
                                 12.0, 12.0,-25.0 };

    double[] expected        = { 12.0, 12.0,
                                -12.0, 12.0,
                                -12.0,-12.0,
                                 12.0,-12.0,
                                 12.0, 12.0 };

    int oldLength = testFromLLAToLL.length;
    double[] testLL   = GeodeticUtil.fromLLAToLL(testFromLLAToLL);

    // test length of array, should be 2/3 that of original
    assertEquals("testLL.length", 2*oldLength/3, testLL.length);

    for (int i = 0; i<expected.length; i++) {
      assertEquals("testLLAToLL[" + i + "]", expected[i], testLL[i]);
    }

  }

  public void testFromLLToALL() {
    double[] testFromLLToALL = { 12.0, 12.0,
                                -12.0, 12.0,
                                -12.0,-12.0,
                                 12.0,-12.0,
                                 12.0, 12.0 };

    double[] expected        = { 0.0, 12.0, 12.0,
                                 0.0, 12.0,-12.0,
                                 0.0,-12.0,-12.0,
                                 0.0,-12.0, 12.0,
                                 0.0, 12.0, 12.0, };

    int oldLength = testFromLLToALL.length;
    double[] testALL   = GeodeticUtil.fromLLToALL(testFromLLToALL);

    // test length of array, should be 3/2 that of original
    assertEquals("testALL.length", 3*oldLength/2, testALL.length);

    assertArrayEquals("testLLToALL[]", expected, testALL);
  }

  public void testGetCenterECEFGeo()               { } // tested in testgetCenterECRGeo()

  public void testGetCenterECRGeo() {
    subtest_getCenterEG( -67.89,  12.345,  8.3, 90);
    subtest_getCenterEG(-177.89,  12.345,  8.3, 90); // cross the IDL
    subtest_getCenterEG( 177.89, -12.345, 16.3, 90); // cross the IDL and the equator
    subtest_getCenterEG( 177.89, -44.345,  8.3, 90); // cross the IDL near 45 degrees lat
    subtest_getCenterEG(   7.89,  84.345,  8.3, 90); // encompass the north pole
    subtest_getCenterEG(  -7.89, -84.345,  8.3, 90); // encompass the south pole
  }
  private void subtest_getCenterEG(double lon, double lat, double radius, int div) {
    double[] centerEG =
      GeodeticUtil.getCenterECRGeo(
        GeodeticUtil.fromLLToALL(
          GeodeticUtil.computeCircleLL(lon, lat, radius, div)));
    assertEquals("centerEG[0]", 0.0, centerEG[0]);
    assertEquals("centerEG[1]", lat, centerEG[1], 0.05);
    assertEquals("centerEG[2]", lon, centerEG[2], 0.00000000000005);
  }

  public void testGetCircleOfLatitudeECEFGeo()     { } // tested in GetCircleOfLatitudeECEFGeo()

  public void testGetCircleOfLatitudeECRGeo() {
    double[] expected        = { 0.0, 0.0,    0.0,
                                 0.0, 0.0,   40.0,
                                 0.0, 0.0,   80.0,
                                 0.0, 0.0,  120.0,
                                 0.0, 0.0,  160.0,
                                 0.0, 0.0, -160.0,
                                 0.0, 0.0, -120.0,
                                 0.0, 0.0,  -80.0,
                                 0.0, 0.0,  -40.0,
                                 0.0, 0.0,    0.0, };

    double[] equatorEG = GeodeticUtil.getCircleOfLatitudeECRGeo(0.0, 9);

    // test length of array, should be 3/2 that of original
    assertEquals("equatorEG.length", expected.length, equatorEG.length);

    assertArrayEquals("equatorEG[]", expected, equatorEG);
  }

  public void testZeroAltitudeLLA() {
    double[] testZeroAlt = { 12.0, 12.0, 25.0,
                            -12.0, 12.0, 25.0,
                            -12.0,-12.0, 25.0,
                             12.0,-12.0, 25.0,
                             12.0, 12.0,-25.0 };
    double   expected = 0.0;

    double[] testZeroAltReturn = GeodeticUtil.zeroAltitudeLLA(testZeroAlt);

    // test both returned array and array passed in (changes are made in-place and returned).
    for (int i = 0; i < 5; i++) {
      assertEquals("testZeroAlt[" + (i*3+2) + "]", expected, testZeroAlt[i*3+2]);
    }

    for (int i = 0; i < 5; i++) {
      assertEquals("testZeroAltReturn[" + (i*3+2) + "]", expected, testZeroAltReturn[i*3+2]);
    }
  }

  /**
   * Test method that will test {@link GeodeticUtil#computeDistanceMeters(Position, Position)} Note: The expected
   * results of the tests are based off what they produced when created. The results have been compared to the results
   * on the Internet and have been found to be similar. If the tests fail, it is due to change in code making the
   * calculation change.
   * @since NeXtMidas 3.8.0
   */
  public void testComputeDistanceMeters () {
    Position position1 = new Position();
    Position position2 = new Position();
    Position position3 = new Position();
    position3.setLat(48.85837);
    position3.setLon(2.2944813);

    // close spots
    position1.setLat(38.889931);
    position1.setLon(-77.009003);
    position2.setLat(38.889484);
    position2.setLon(-77.035278);
    testComputeDistanceMetersHelper(position1, position2, 2283.180178298971);

    //Distance to Eiffel Tower from US Capitol
    testComputeDistanceMetersHelper(position1, position3,6187781.52671396);

    // same spot
    position1.setLat(90);
    position1.setLon(90);
    position2.setLat(90);
    position2.setLon(90);
    testComputeDistanceMetersHelper(position1, position2, 0);

    // far spots
    position1.setLat(90);
    position1.setLon(90);
    position2.setLat(-90);
    position2.setLon(-90);
    testComputeDistanceMetersHelper(position1, position2, 20037508.342789244);

  }

  /**
   * This is a helper function for testComputeDistanceMeters. It does the calculations and asserts for each test case.
   * @param position1 - first position
   * @param position2 - second position
   * @param expected - expected results
   * @since NeXtMidas 3.8.0
   */
  private void testComputeDistanceMetersHelper (Position position1, Position position2, double expected) {
    double newToMeters = GeodeticUtil.computeDistanceMeters(position1, position2);
    assertEquals("The GeodeticUtil.computeDistanceMeters should calculate correctly", newToMeters, expected);
  }

  /**
   * This will test {@link GeodeticUtil#computeDistance(Position, Position)} Note: The expected results of the tests are
   * based off what they produced when created. The results have been compared to the results on the Internet and have
   * been found to be similar. If the tests fail, it is due to change in code making the calculation change.
   * @since NeXtMidas 3.8.0
   */
  public void testComputeDistanceOriginal () {
    Position position1 = new Position();
    Position position2 = new Position();
    Position position3 = new Position();
    position3.setLat(48.85837);
    position3.setLon(2.2944813);

    // close spots
    position1.setLat(38.889931);
    position1.setLon(-77.009003);
    position2.setLat(38.889484);
    position2.setLon(-77.035278);
    testComputeDistanceOriginalHelper(position1, position2, 2277.1722903465484);

    //Distance to Eiffel Tower from US Capitol
    testComputeDistanceOriginalHelper(position1, position3,6167008.697729923);

    // same spot
    position1.setLat(90);
    position1.setLon(90);
    position2.setLat(90);
    position2.setLon(90);
    testComputeDistanceOriginalHelper(position1, position2, 0);

    // far spots
    position1.setLat(90);
    position1.setLon(90);
    position2.setLat(-90);
    position2.setLon(-90);
    testComputeDistanceOriginalHelper(position1, position2, 20037508.342789244);
  }

  /**
   * This is a helper function for {@link GeodeticUtilTest#testComputeDistanceOriginal()} It does the calculations and
   * asserts for each test case.
   * @param position1 - first position
   * @param position2 - second position
   * @param expected - expected results
   * @since NeXtMidas 3.8.0
   */
  @SuppressWarnings("deprecation") // Since 3.9.0
  private void testComputeDistanceOriginalHelper (Position position1, Position position2, double expected) {

    double distance = GeodeticUtil.computeDistance(position1, position2);
    assertEquals("The GeodeticUtil.computeDistance should calculate correctly", distance, expected);
  }

  /**
   * This will test {@link GeodeticUtil#computeDistanceDeg(Position, Position)} Note: The expected results of the tests
   * are based off what they produced when created. The results have been compared to the results on the Internet and
   * have been found to be similar. If the tests fail, it is due to change in code making the calculation change.
   * @since NeXtMidas 3.8.0
   */
  public void testComputeDistanceDeg () {
    Position position1 = new Position();
    Position position2 = new Position();
    Position position3 = new Position();
    position3.setLat(48.85837);
    position3.setLon(2.2944813);

    // close spots
    position1.setLat(38.889931);
    position1.setLon(-77.009003);
    position2.setLat(38.889484);
    position2.setLon(-77.035278);
    computeDistDegHelper(position1, position2, 0.020510156505646998);

    //Distance to Eiffel Tower from US Capitol
    computeDistDegHelper(position1, position3, 55.58578720239577);

    // same spot
    position1.setLat(90);
    position1.setLon(90);
    position2.setLat(90);
    position2.setLon(90);
    computeDistDegHelper(position1, position2, 0);

    // far spots
    position1.setLat(90);
    position1.setLon(90);
    position2.setLat(-90);
    position2.setLon(-90);
    computeDistDegHelper(position1, position2, 180.0);

  }

  /**
   * This is a helper function for {@link GeodeticUtilTest#testComputeDistanceDeg()} It does the calculations and
   * asserts for each test case.
   * @param position1 - first position
   * @param position2 - second position
   * @param expected - expected results
   * @since NeXtMidas 3.8.0
   */
  private void computeDistDegHelper (Position position1, Position position2, double expected) {
    double degree = GeodeticUtil.computeDistanceDeg(position1, position2);
    assertEquals("The GeodeticUtil.computeDistanceDeg should calculate correctly", degree, expected);
  }

  /**
   *  Test created for using non normalized values - i.e longitude of 185 - normalized is -175
   */
  public void testInsidePolygonNormalize () {
    double DATELINE_LONGITUDE = 180.0;
    double EQUATOR_LATITUDE = 0.0;
    double RADIUS_DEGREES = 10.0;
    int NUMBER_OF_SEGMENTS = 4;
    double TEST_POINT_1_LATITUDE = 0.0;
    double TEST_POINT_1_LONGITUDE = -175.0;
    double TEST_POINT_2_LATITUDE = 0.0;
    double TEST_POINT_2_LONGITUDE = 185.0;
    double TEST_POINT_3_LATITUDE = 0.0;
    double TEST_POINT_3_LONGITUDE = -108175; // -175 + (-360 * 300)
    // Generate a circle that crosses the International Date Line:
    double[] circleContour = GeodeticUtil.computeCircleLL(
        DATELINE_LONGITUDE, // Center the circle on the date line
        EQUATOR_LATITUDE,
        RADIUS_DEGREES,
        NUMBER_OF_SEGMENTS);
    // Confirm that those points, which are not normalized, are considered inside the polygon.
    assertTrue(isPointIncluded(TEST_POINT_1_LATITUDE, TEST_POINT_1_LONGITUDE, circleContour));
    assertTrue(isPointIncluded(TEST_POINT_2_LATITUDE, TEST_POINT_2_LONGITUDE, circleContour));
    assertTrue(isPointIncluded(TEST_POINT_3_LATITUDE, TEST_POINT_3_LONGITUDE, circleContour));
  }

  private boolean isPointIncluded (double lat, double lon, double[] contour) {
    Position position = Position.fromGeo(0, lat, lon);
    if (GeodeticUtil.insidePolygon(position, 2, contour)) {
      return true;
    }
    else {
      return false;
    }
  }

  /** This will test all versions of getCountryInfo. Countries are returned based on inclusion over nearness.
   *  This mean that if a point is inside a country, only that country will be returned regardless of the radius.
   *  If a point is in the water, all countries within the radius will be returned.
   *  @since NeXtMidas 3.9.4, modified in NeXtMidas 4.1.2
   */
  public void testGetCountry () {
    // In Slovenia (zero radius)
    List<Table> countryInfo = GeodeticUtil.getCountryInfo(0, 46.52615215310537, 13.718178651141322, 0);
    assertEquals("GeodeticUtil.getCountryInfo Slovenia:","Slovenia",namesOfCountries(countryInfo));

    // In Slovenia (no radius provided, default radius reaches both Slovenia and surrounding countries)
    countryInfo = GeodeticUtil.getCountryInfo(0, 46.52615215310537, 13.718178651141322);
    assertEquals("GeodeticUtil.getCountryInfo Slovenia:","Slovenia",namesOfCountries(countryInfo));

    // In water, outside Taiwan (no radius provided, default radius does not reach Taiwan, using Position)
    countryInfo = GeodeticUtil.getCountryInfo(Position.fromGeo(0,24,122));
    assertEquals("GeodeticUtil.getCountryInfo []:","",namesOfCountries(countryInfo));

    // In water, outside Taiwan (radius provided reaches Taiwan)
    countryInfo = GeodeticUtil.getCountryInfo(0, 24, 122, 1e5);
    assertEquals("GeodeticUtil.getCountryInfo Taiwan:","Taiwan",namesOfCountries(countryInfo));

    // In water, outside Mexico and the United States (radius provided reaches both Mexico and the United States)
    countryInfo = GeodeticUtil.getCountryInfo(0, 26.048, -97.15, 1e4);
    assertEquals("GeodeticUtil.getCountryInfo Mexico,United States:","Mexico,United States",namesOfCountries(countryInfo));

    // In water, outside Mexico and the United States (radius provided does not reach Mexico nor the United States, non-existent field)
    List<String> countryInfoString = GeodeticUtil.getCountryInfo("NONEXISTENT",0, 26.048, -97.15, 1);
    assertEquals("GeodeticUtil.getCountryInfo [] (non-existent field):","[]",countryInfoString.toString());

    // In water, outside Mexico and the United States (no radius provided, default radius reaches both Mexico and the United States, non-existent field)
    countryInfoString = GeodeticUtil.getCountryInfo("NONEXISTENT",0,26.048,-97.15);
    assertEquals("GeodeticUtil.getCountryInfo Mexico,United States (non-existent field):","[null, null]",countryInfoString.toString());

    // In water, outside Mexico and the United States (no radius provided, default radius reaches both Mexico and the United States, NAME field, using Position)
    countryInfoString = GeodeticUtil.getCountryInfo("NAME",Position.fromGeo(0,26.048,-97.15));
    assertEquals("GeodeticUtil.getCountryInfo Mexico,United States (NAME field):","[Mexico, United States]",countryInfoString.toString());
  }

  private String namesOfCountries (List<Table> onesNear) {
    String countryNames = "";
    boolean waitingOnFirst = true;
    for (Table table : onesNear) {
      String thisCountry = table.get("NAME").toString();
      if (waitingOnFirst) {
        countryNames = thisCountry;
        waitingOnFirst = false;
      }
      else {
        countryNames += "," + thisCountry;
      }
    }
    return countryNames;
  }

  public void testComputeXYZEllipsoidCrossSections () {
    double semiAxisOrig = 10;
    int numSlices = 3;
    double[][][] circleCross = GeodeticUtil.computeXYZEllipsoidCrossSections(semiAxisOrig,semiAxisOrig,semiAxisOrig,numSlices);

    double ratio = 0.4;
    double semiDirectionNew = semiAxisOrig*ratio;
    double semiOthersNew = semiAxisOrig*Math.sqrt(1-Math.pow(ratio, 2));
    double[][][] expected = {{{0,semiAxisOrig,semiAxisOrig},{semiDirectionNew,semiOthersNew,semiOthersNew},{-semiDirectionNew,semiOthersNew,semiOthersNew}},
        {{semiAxisOrig,0,semiAxisOrig},{semiOthersNew,semiDirectionNew,semiOthersNew},{semiOthersNew,-semiDirectionNew,semiOthersNew}},
        {{semiAxisOrig,semiAxisOrig,0},{semiOthersNew,semiOthersNew,semiDirectionNew},{semiOthersNew,semiOthersNew,-semiDirectionNew}}};
    String[] axes = {"X","Y","Z"};
    for (int i = 0; i < circleCross.length ; i++) {
      for (int slice=0; slice < circleCross[0].length; slice++) {
        assertEquals("Discretize "+axes[i]+", new semi-major:",circleCross[i][slice][0],expected[i][slice][0]);
        assertEquals("Discretize "+axes[i]+", new semi-minor:",circleCross[i][slice][1],expected[i][slice][1]);
        assertEquals("Discretize "+axes[i]+", new semi-median:",circleCross[i][slice][2],expected[i][slice][2]);
      }
    }
  }

  public void testEep2cep () {
    // The eep2cep(..) function is a port from Python. These test cases are
    // just comparing expected input/output values to what we get when we
    // use the Python version which has been around for some time (i.e. we
    // just assume the Python version is correct).
    double tol = 5e-9;
    testEep2cep(0.0,     0.0,     0.0,           tol);
    testEep2cep(1e3,     0.0,     790.73559226,  tol);
    testEep2cep(1e3,     1e3,     995.380339255, tol);
    testEep2cep(2e3,     1e3,     1658.53654023, tol);
    testEep2cep(75.6154, 42.0468, 63.5011585632, tol);
  }

  private void testEep2cep (double semiMaj, double semiMin, double cep, double tol) {
    assertEquals("eep2cep("+semiMaj+", "+semiMin+") -> "+cep, GeodeticUtil.eep2cep(semiMaj, semiMin), cep, tol);
    assertEquals("eep2cep("+semiMin+", "+semiMaj+") -> "+cep, GeodeticUtil.eep2cep(semiMin, semiMaj), cep, tol);
  }
}
