package nxm.sys.lib;

import static java.lang.Math.PI;
import static java.lang.Math.acos;
import static java.lang.Math.asin;
import static java.lang.Math.atan;
import static java.lang.Math.atan2;
import static java.lang.Math.cos;
import static java.lang.Math.max;
import static java.lang.Math.min;
import static java.lang.Math.pow;
import static java.lang.Math.sin;
import static java.lang.Math.sqrt;
import static java.lang.Math.tan;
import static java.lang.Math.toDegrees;
import static java.lang.Math.toRadians;
import static java.util.Collections.reverse;
import static java.util.stream.Collectors.toList;
import static nxm.sys.inc.Constants.DEG2RAD;
import static nxm.sys.inc.Constants.EARTH_EQUATORIAL_RADIUS;
import static nxm.sys.inc.Constants.METERS2DEG;
import static nxm.sys.inc.Constants.RAD2DEG;
import static nxm.sys.inc.Constants.SPEED_OF_LIGHT;
import static nxm.sys.inc.Constants.TWOPI;
import static nxm.sys.lib.Position.fromGeo;

import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

import nxm.sys.inc.InternalUseOnly;
import nxm.sys.inc.MidasReference;
import nxm.sys.inc.ProvisionalUseOnly;
import nxm.sys.inc.Units;


/** UNCLASSIFIED - Utilities for geodetic computations. <br>
    <br>
    <b>Naming Convention:</b>
    All methods that end in <code>LL</code> return arrays that are a series of <code>{LON, LAT}</code>
    points where <code>ALT</code> can always be assumed to be <code>0.0</code>. Those that end in
    <code>LLA</code> return arrays that are a series of <code>{LON, LAT, ALT}</code> points. If you need
    to use <code>ALL</code> (<code>{ALT, LAT, LON}</code> )form, just use {@link Transform#all2lla} to
    do the conversion (note that <code>ALL-&gt;LLA</code> and <code>LLA-&gt;ALL</code> are implemented the
    same as both are just swapping the position of the <code>ALT</code> and <code>LON</code> components).<br>
    <br>
    Methods that end in <code>TOP</code> return a series of topocentric <code>{X, Y, X}</code> values
    (see {@link Transform#TOP}).<br>
    <br>
    <b>Units:</b>
    Unless specified otherwise, the following units should be assumed:
    <pre>
      LON      - Degrees longitude
      LAT      - Degrees latitude
      ALT      - Meters above geoid (WGS-84 or "spherical earth").
      AZIM     - Azimuth, relative to observer,  in degrees on a range of 0 to 360 (0=North,90=East,etc).
      HEADING  - Heading, relative to moving object, in degrees on a range of 0 to 360 (0=North,90=East,etc).
    </pre>
    Note that <b>distances on the surface of the Earth may be measured in degrees or meters</b> and
    it is important to check the docs for each method to see which is being used. Typically small
    distances are measured in meters while large ones are in degrees (this is a common convention).<br>
    <br>
    <b>Error Conditions:</b>
    If illegal input values are given, the behavior of the methods in this class is <b>undefined</b>
    and may change without notice (many methods skip error checking in favor of speed). <b>It is the
    responsibility of the caller to make sure that input values are legal.</b>
    @author  J. Czechowski IV
    @author  M. Scott Meuleners
    @classification UNCLASSIFIED
    @since NeXtMidas 2.7.2
 */
public final class GeodeticUtil {
  // JC4: Most of the code in here was taken from various places in the SYS tree (mostly private
  //      methods or functionality within another method). Wherever I did that, I made a note of
  //      it after the method declaration. The intent being that we can put this class in SYS and
  //      replace the old private methods with it.

  /** Time Difference Of Arrival */
  public static final int MODE_TDOA = 1;
  /** Frequency Difference Of Arrival */
  public static final int MODE_FDOA = 2;

  private static ShapeFile countries;
  /** The ordering of the points in an array. */
  public static enum ArrayForm {
      /** array ordered based on LonLat */
      LonLat,
      /** array ordered based on LonLatAlt */
      LonLatAlt,
      /** array ordered based on LatLon */
      LatLon,
      /** array ordered based on LatLonAlt */
      LatLonAlt,
      /** array ordered based on AltLatLon */
      AltLatLon
    }

  /** Fix the case where adjacent latitudes are more than 180 degrees apart - should count as crossing */
  private static boolean insidePolyFix = true;

  private GeodeticUtil () { } // prevent instantiation since all methods are static

  /** Returns whether a point p is inside a given polygon, that may have "holes" in it. For example:
      <pre>
           ###
          #   ####           S
         #        ######################
        #    P                          #
        #               %%%%%%%          #
        #    $$$$$$    %   R   %         #              # = outer
        #   $  Q   $    %%%%%%%         #               $ = inner[0]
        #    $$$$$$                #####                % = inner[1]
         ###       ################
            #######
      </pre>
      In the above example, there is a polygon (<code># = outer</code>) with two "holes"
      (<code>$ = inner[0]</code> and <code>% = inner[1]</code>). Additionally there are four possible points
      (P, Q, R, S). Point P is inside the polygon, point Q is outside the polygon since it is within
      one of the "holes", point R is also outside the polygon since it is also within a "hole", and
      point S is outside the polygon since it is not contained within the outer boundary. <br>
      <br>
      This method uses {@link #insidePolygon(Position,int,double[])} for most of its computations.
      @param p     Point to check for inclusion.
      @param dim   Dimensionality of the outer and inner arrays (should be 2 for LL, or 3 for LLA)
      @param outer Points that make up the outer boundary in LL or LLA format.
      @param inner Array of points that make up the inner "holes" in LL or LLA format. Each inner
                   boundary must be completely contained within the outer boundary and no two inner
                   boundaries may overlap.
      @return      Is point p is inside the given polygon
   */
  public static boolean insidePolygon (Position p, int dim, double[] outer, double[] ... inner) {
    if (inner != null) {
      for (double[] in : inner) {
        if (insidePolygon(p, dim, in)) return false; // within a "hole"
      }
    }
    return insidePolygon(p, dim, outer);
  }

  /** This started as an optimized version of the method in the NeXtMidas GeodeticUtil class
   *  which determines if a point is included in the polygon, but then was updated to fix issues
   *  when the point is on the edge of the polygon, etc.
   *  @param point   The point.
   *  @param polygon The polygon (must be closed).
   *  @return true if the point is in the polygon, false otherwise.
   *  @since NeXtMidas 4.0.0 for NXM-3727
   */
  @ProvisionalUseOnly("limited testing")
  public static boolean insidePolygon (Position point, List<Position> polygon) {
    boolean a        = false;
    boolean b        = false;
    boolean pole     = false;
    double pointLat = point.getLat();
    double pointLon = point.getLon();
    Position p0       = polygon.get(0);
    double priorLat = p0.getLat();
    double priorLon = p0.getLon();
    double maxLat   = max(priorLat, 0);
    double minLat   = min(priorLat, 0);

    for (int n = 1; n < polygon.size(); n++) {
      Position pN      = polygon.get(n);
      double thisLat = pN.getLat();
      double thisLon = pN.getLon();
      double origThisLon = thisLon;
      double lastLat = priorLat;
      double lastLon = priorLon;
      double thisLastDiff = thisLon - lastLon; // NXM-3273 - spans greater than 180 indicate crossing

      // Check to see if there is a -180/+180 crossing here. If so we adjust the
      // numbers so polygon edge crossing detector will work. We also toggle the
      // pole flag since an off number of crossings means it surrounds a pole.
      if (thisLastDiff > 180) {
        pole = !pole;
        // if point to be checked is near -180/180 boundary, adjust numbers so
        // polygon edge crossing detector will work
        if (pointLon > thisLon)
          lastLon += 360.0;
        else
          thisLon -= 360.0;
      }
      else if (thisLastDiff < -180) {
        pole = !pole;
        // if point to be checked is near -180/180 boundary, adjust numbers so
        // polygon edge crossing detector will work
        if (pointLon > lastLon)
          thisLon += 360.0;
        else
          lastLon -= 360.0;
      }

      // Edge crossing detection...
      if ((thisLat == pointLat) && (thisLon == pointLon)) {          // On a Vertex
        return true;
      }
      else if ((thisLon == pointLon) && (thisLon == lastLon)) {      // Vertical line
        if (((thisLat <= pointLat) && (pointLat <= lastLat)) ||
            ((lastLat <= pointLat) && (pointLat <= thisLat))) {
          return true;
        }
      }
      else if ((thisLat == pointLat) && (thisLat == lastLat)) {     // Horizontal line
        if (((thisLon <= pointLon) && (pointLon <= lastLon)) ||
            ((lastLon <= pointLon) && (pointLon <= thisLon))) {
          return true;
        }
      }
      else if ((lastLon >= pointLon) != (thisLon >= pointLon)) {    // Crossing detector
        if (pointLat > (thisLat - lastLat)*(pointLon-lastLon) / (thisLon - lastLon) + lastLat) {
          a = !a; // Line drawn north from the point crosses edge of the polygon
        }
        else {
          b = !b; // Line drawn south from the point crosses edge of the polygon
        }
      }
      maxLat   = max(thisLat, maxLat);
      minLat   = min(thisLat, minLat);
      priorLat = thisLat;
      priorLon = origThisLon;
    }

    // determine which pole is surrounded (if crossing=true)
    if (!pole) {
      return a; // Not a pole (will have a=b so we just return a)
    }
    else if (maxLat > -minLat) {
      return a; // North pole
    }
    else {
      return b; // South pole
    }
  }

  /** Returns whether a point p is inside a given polygon. <br>
      <br>
      insidePolygon calculations are based on those in PNPOLY - Point Inclusion in Polygon Test from W. Randolph
      Franklin <a href="http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html">
      http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html</a>. Algorithm was modified to check
      intersection on the y axis instead of the x axis, and Checks were added for polygons that cross over the
      -180/180 degree boundary as well as for polygons that surround the north and south pole.
      <br>
      This method assumes that the polygon will fit inside the surface of half of a sphere.
      @param p       Point to check for inclusion.
      @param dim     Dimensionality of the polygon array (should be 2 for LL, or 3 for LLA)
      @param polygon Points that make up the outer boundary of the polygon in LL or LLA format.
                     (Lon,Lat) or (Lon, Lat, Alt)
      @return        Is point p is inside the given polygon
   */
  public static boolean insidePolygon (Position p, int dim, double[] polygon) {
    boolean a=false,b=false,crossing=false;
    double maxy = 0.0;
    double miny = 0.0;
    int i=0,j=0;
    int length = polygon.length / dim;
    double py  = p.getLat();
    double px  = p.getLon();
    // Normalize the longitude  Since 3.8.1 NXM-3253
    px = normalizeLongitude(px);
    // Need to make a copy Since 3.8.1 NXM-3253
    double[] polygonCopy = polygon.clone();
    // Normalize the longitudes in the array Since 3.8.1 NXM-3253
    for (i = 0; i < (polygonCopy.length / dim); i++) {
      double lon = polygonCopy[(i * dim)];
      polygonCopy[(i * dim)] = normalizeLongitude(lon);
    }

    for (i = 0, j = length-1; i < length; j=i++) {
      double pix = polygonCopy[i*dim];
      double pjx = polygonCopy[j*dim];
      double piy = polygonCopy[i*dim+1];
      double pjy = polygonCopy[j*dim+1];
      maxy = max(piy,maxy);
      miny = min(piy,miny);
      double pijdiff = pix - pjx; // NXM-3273 - spans greater than 180 indicate crossing
      // check to see if a -180/180 crossing happens (an odd number of crossings
      // means that the polygon surrounds a pole)
      if ((pix>+90.0 && pjx<-90.0) || (insidePolyFix && (pijdiff > 180))) {
        crossing = !crossing;
        // if point to be checked is near -180/180 boundary, adjust numbers so
        // polygon edge crossing detector will work
        if (px > pix)
          pjx += 360.0;
        else
          pix -= 360.0;
      }
      else if ((pix<-90.0 && pjx>+90.0) || (insidePolyFix && (pijdiff < -180))) {
        crossing = !crossing;
        // if point to be checked is near -180/180 boundary, adjust numbers so
        // polygon edge crossing detector will work
        if (px > pjx)
          pix += 360.0;
        else
          pjx -= 360.0;
      }

      // Edge crossing detection...
      if ((pjy == py) && (pjx == px)) {          // On a Vertex
        return true;
      }
      else if ((pjx == px) && (pjx == pix)) {      // Vertical line
        if (((pjy <= py) && (py <= piy)) ||
            ((piy <= py) && (py <= pjy))) {
          return true;
        }
      }
      else if ((pjy == py) && (pjy == piy)) {     // Horizontal line
        if (((pjx <= px) && (px <= pix)) ||
            ((pix <= px) && (px <= pjx))) {
          return true;
        }
      }
      else if ((pix>=px) != (pjx>=px)) {
        // detect if a line drawn north from target crosses an edge of
        // the polygon
        if (py >  (pjy - piy)*(px-pix) / (pjx - pix) + piy)
          a = !a;
        // detect if a line drawn south from target crosses an edge of
        // the polygon
        if (py <= (pjy - piy)*(px-pix) / (pjx - pix) + piy)
          b = !b;
      }
    }

    // determine which pole is surrounded (if crossing=true)
    boolean north = maxy > (-miny);

    return (a && !crossing) || (a && crossing && north) || (b && crossing && !north);
  }


  /**
   * Helper function to convert non normalized longitude values (i.e. greater than 180 or less than -180)
   * @param lon The longitude to be normalized. If the longitude is good, it will not change the value.
   * @return The normalized value
   */
  public static double normalizeLongitude (double lon) {
    // Need to adjust the points to be normalized Since 3.8.1 NXM-3253
    if ( lon < -180) {
      lon = -(((-lon + 180) % 360) - 180);
    }
    else if (lon > 180) {
        lon = ((lon + 180) % 360) - 180;
      }
    return lon;
  }

  /** Calculates the TDOA or FDOA from a pair of positions to another position
   *
   * @param p1          Position of first observer
   * @param p2          Position of second observer
   * @param t1          Position of target
   * @param mode        Mode ({@link #MODE_TDOA} or {@link #MODE_FDOA})
   * @param soi_freq    Frequency (used only for MODE_FDOA)
   */
  private static double predictPair (Position p1, Position p2, Position t1, int mode, double soi_freq) {

    double xta,yta,zta,xtb,ytb,ztb,rta,rtb,rtai,rtbi,rrta,rrtb;
    double cci;

    // calc distances/rates to target
    xta = p1.getX() - t1.getX();   yta = p1.getY() - t1.getY();   zta=p1.getZ() - t1.getZ();
    xtb = p2.getX() - t1.getX();   ytb = p2.getY() - t1.getY();   ztb=p2.getZ() - t1.getZ();
    rta = Math.sqrt(xta*xta + yta*yta + zta*zta); // calc distance from a to target
    rtai=1/rta; // inverse of rta
    rtb = Math.sqrt(xtb*xtb + ytb*ytb + ztb*ztb); // calc distance from b to target
    rtbi=1/rtb; // inverse of rtb
    rrta = ( xta*p1.getVX() + yta*p1.getVY() + zta*p1.getVZ() ) * rtai; // calculate rate of position change between a and target
    rrtb = ( xtb*p2.getVX() + ytb*p2.getVY() + ztb*p2.getVZ() ) * rtbi; // calculate rate of position change between b and target

    // calculate tdoa/fdoa
    cci = 1/SPEED_OF_LIGHT;

    if (mode == MODE_TDOA) return (cci*(rtb - rta));
    else                   return (-cci*( soi_freq*(rrtb-rrta)));
  }

  /** A test to determine whether a double is between two other doubles.
   *
   * @param test    value to test
   * @param a       first value to test against
   * @param b       second value to test against
   * @return        true if value is between a and b, false if not
   */
  private static boolean isBetween(double test, double a, double b) {
    return (((test > a) && (test < b)) || ((test < a) && (test > b)));
  }

  /**
   * Returns distance in degrees between two positions on the surface of the earth. Note: This does not take into
   * account the earth flattening factor. see {@link GeodeticUtil#computeDistanceMeters} and
   * {@link GeodeticUtil#computeDistanceDeg}
   * @param p1 First Position
   * @param p2 Second Position
   * @return The distance in meters latitude between p1 and p2.
   * @deprecated Since NeXtMidas 3.9.0. Please use the new functionality with using either
   *             {@link GeodeticUtil#computeDistanceMeters(Position, Position)} or
   *             {@link GeodeticUtil#computeDistanceDeg(Position, Position)}
   */
  @Deprecated
  public static double computeDistance(Position p1, Position p2) {
    double lat1 = p1.lat;
    double lat2 = p2.lat;
    double lon1 = p1.lon;
    double lon2 = p2.lon;
    return EARTH_EQUATORIAL_RADIUS * (acos(sin(DEG2RAD*lat1)*sin(DEG2RAD*lat2) + cos(DEG2RAD*lat1)*cos(DEG2RAD*lat2)*cos(DEG2RAD*lon1-DEG2RAD*lon2)));
  }

  /** Returns the great arc distance between two positions in meters,  Note that this method calls
   * {@link Transform#geo2car(Position)} on both positions prior to computing the distance.
   * @param position1 First Position. (Must be in {@link Units#GEODETIC} mode).
   * @param position2 Second Position. (Must be in {@link Units#GEODETIC} mode).
   * @return The distance in meters between p and q
   * @since NeXtMidas 3.8.0
   */
  public static double computeDistanceMeters (Position position1, Position position2) {
    return Transform.greatArc(position1,position2);
  }

  /** Returns the great arc distance between two positions in degrees latitude,  Note that this method calls
   * {@link Transform#geo2car(Position)} on both positions prior to computing the distance.
   * @param p First Position. (Must be in {@link Units#GEODETIC} mode).
   * @param q Second Position. (Must be in {@link Units#GEODETIC} mode).
   * @return The distance in degrees latitude between p and q
   * @since NeXtMidas 3.8.0
   */
  public static double computeDistanceDeg (Position p, Position q) {
    return METERS2DEG * Transform.greatArc(p,q);
  }

  /** Function to generate TDOA / FDOA Lines of Position (LOPs)

     <pre>
          %                       % = First Position
           #                      # = Subsequent Positions
            #xxx                  x = Considered for Next Position
            x#  x
           x  #  x
            x   x
             xxx
     </pre>

     The above example illustrates the method in which the LOPs are calculated.
     The first Position is passed in to the function, <code>t</code>.  The algorithm
     attempts to trace a line where all subsequent points have the same TDOA
     or FDOA as <code>t</code>.  For each subsequent point, TDOA / FDOA values are
     calculated for positions around the current point until the appropriate point
     is found.  This process continues until the line drawn returns to the starting
     <code>t</code>.  Various steps are taken to optimize the process.<br>
     <br>
     Note: Currently limited to  approximately 2450 coordinate pairs
     This method uses <code>predictPair (Position, Position, Position, int, double)</code>
     for most of its computations.

     @param p1        Position of first observer
     @param p2        Position of second observer
     @param t         Position of First point on LOP
     @param soi_freq  Signal of interest frequency (required for FDOA only).
     @param mode      TDOA or FDOA (use MODE_TDOA or MODE_FDOA)
     @param aRes      Maximum number of points to consider for each subsequent point (at
                      least 512 is recommended, actual number of points considered
                      for each Position is usually considerably less.)
     @param vRes      Array of distance resolution in DR form, where D is the distance
                      from the target, and R is the distance resolution for all points
                      less than that distance away from the target.  It is assumed that
                      all points will be pre-sorted by increasing D values.  If a null
                      value is entered, a default array will be used.
     @return          Array of double in LL format.
   */
  public static double[] drawLOPLL(Position p1, Position p2, Position t, double soi_freq, int mode, int aRes, double[] vRes) {
    //int nPoints = 256;
    double tPoint         = predictPair(p1,p2,t,mode,soi_freq); // calculate T/FDOA at position t
    double[] pointsArray  = new double[aRes];
    double div            = 360.0/aRes;
    Position cp           = t.copy();
    double[] buf          = new double[5000];
    if (vRes == null)vRes = new double[] { 100000.0, 500.0, 1000000.0, 10000.0, 4000000.0, 50000.0, Double.MAX_VALUE, 400000.0 };
    double distRes        = vRes[1];

    // record initial position
    buf[0] = cp.lon;
    buf[1] = cp.lat;

    int j = 2;
    int l = 0;
    double distToT = 0;

    while (!(((2*vRes[1] > distToT) && (j>10)) || (j>(buf.length-10)))) {
      double closestValue = Double.MAX_VALUE;
      int closest = 0;
      boolean found = false;
      double angle = 0;

      // variable resolution
      for (int i=0; i<vRes.length; i+=2) {
        if (distToT < vRes[i]) {
          distRes = vRes[i+1];
          break;
        }
      }

      int previousPoint = (l+aRes)%aRes;
      pointsArray[previousPoint] = predictPair(p1,p2,computePosition(cp, div*(previousPoint), distRes),mode,soi_freq);

      for (int i=0; i<aRes; i++) {
        int currentPoint;
        int pPoint;
        if (i%2==0) {
          currentPoint = (l+i/2 + 1)%aRes;
          pPoint = (currentPoint+aRes-1)%aRes;
        }
        else {
          currentPoint = (l+aRes-i/2-1)%aRes;
          pPoint = (currentPoint + 1)%aRes;
        }

        pointsArray[currentPoint] = predictPair(p1,p2,computePosition(cp, div*currentPoint, distRes),mode,soi_freq);

        if (isBetween(tPoint, pointsArray[currentPoint], pointsArray[pPoint])) {
          // compute angle to travel
          if (((pointsArray[currentPoint] >= pointsArray[pPoint]) && (currentPoint > pPoint)) || ((pointsArray[currentPoint] <= pointsArray[pPoint]) && (currentPoint < pPoint)))
            angle = div*currentPoint + div*Math.abs(tPoint - pointsArray[pPoint])/(pointsArray[currentPoint] - pointsArray[pPoint]);
          else
            angle = div*currentPoint + div*Math.abs(pointsArray[pPoint] - tPoint)/(pointsArray[pPoint] - pointsArray[currentPoint]);
            found = true;
            l = currentPoint;
          break;
        }

        if (Math.abs(tPoint - pointsArray[currentPoint]) < closestValue) {
          closestValue = Math.abs(tPoint - pointsArray[currentPoint]);
          closest = currentPoint;
        }
      }

      // if point was not found between two other points, use the one that was
      // closest in value
      if (!found) {
        angle = div*closest;
        // calculate where to start on next loop
        l = closest;
      }

      // calculate next position
      Position np = computePosition(cp, angle, distRes);

      // calculate start direction for next loop
      if (np.lat >= 90) l = aRes/2;
      else if (np.lat <= -90) l = 0;
      // compute course from next position to last position, then go the
      // opposite direction
      else l = ((int)(computeCourse(np,cp) * aRes / 360.0) + aRes/2) % aRes;

      // copy next position to current position
      cp = np;

      // compute the distance to the position t
      distToT = computeDistance(cp,t);

      buf[j++] = cp.lon;
      buf[j++] = cp.lat;
    }

    // Connect first position to last position
    buf[j++] = buf[0];
    buf[j++] = buf[1];

    double[] returnArray = new double[j];
    System.arraycopy(buf, 0, returnArray, 0, j);
    return returnArray;
  }

  /** Computes the initial course (azimuth) from <code>p1</code> to <code>p2</code> in degrees.
      Note: This uses the algorithm for 'great circle', not 'flat Earth'
      @param p1 Current Position
      @param p2 Target Position
      @return initial course in degrees from 0.0 - 360.0
   */
  public static double computeCourse(Position p1, Position p2) {
    double p1lonrad = p1.lon * DEG2RAD;
    double p2lonrad = p2.lon * DEG2RAD;
    double p1latrad = p1.lat * DEG2RAD;
    double p2latrad = p2.lat * DEG2RAD;
    double tc = atan2(sin(p2lonrad-p1lonrad)*cos(p2latrad),cos(p1latrad)*sin(p2latrad)-sin(p1latrad)*cos(p2latrad)*cos(p2lonrad-p1lonrad));
    tc = RAD2DEG * (tc % (2*PI));
    return (tc+360)%360;  // ensure that value is 0.0 <= tc <= 360.0
  }

  /** Returns the elevation angle from Position <code>q</code> to position <code>p</code>.
      @param p   The viewer.
      @param q   The target location (the point the viewer is looking at).
      @param def The default value (returned if viewer is directly over q.
      @return    The elevation angle.
   */
  public static double computeElevation (Position p, Position q, double def) {
    Transform tp = new Transform();
    Transform tq = new Transform();
    tp.setTOP(p.alt, p.lat, p.lon);
    tq.setTOP(q.alt, q.lat, q.lon);
    double[] qpos = new double[3];  // initialized to 0.0, 0.0, 0.0
    Transform.top2top(qpos, Transform.POS, 1, tp, tq);
    double x = qpos[0];
    double y = qpos[1];
    double z = qpos[2];
    if (z == 0) return def;
    double r1 = sqrt(x*x + y*y);
    return toDegrees(atan(z / r1));
  }

  /** Returns position on earths surface distance <code>d</code> meters away from nadir at heading
      <code>h</code>. (The result will always be on the surface of the Earth.) This method assumes a
      spherical earth, of radius {@link nxm.sys.inc.Constants#EARTH_EQUATORIAL_RADIUS}. Calculations are based
      on those in Ed Williams Aviation Formulary V1.44,
      <a href="http://williams.best.vwh.net/avform.htm">http://williams.best.vwh.net/avform.htm</a>
      @param p  The starting position.
      @param h  The heading in deg (0=North,90=East,etc).
      @param d  The distance in meters to target (on surface of earth, must be positive).
      @return   New resulting position.
   */
  public static Position computePosition (Position p, double h, double d) { // new method
    if (d == 0) return Position.fromGeo(0.0, p.getLat(), p.getLon());
    else if ((p.getLat() >= 90.0) || (p.getLat() <= -90.0)) {
      double distDegLat = d / (TWOPI * EARTH_EQUATORIAL_RADIUS) * 360.0;
      // Handle the special case of being directly over the N/S pole.
      return Position.fromGeo(0.0, (p.getLat()>=90)? (90 - distDegLat): (-90 + distDegLat), h);
    }
    else {
      double hRad     = h * (PI / 180.0);
      double dRad     = d / EARTH_EQUATORIAL_RADIUS;
      double lat1Rad  = DEG2RAD * p.getLat();
      double lon1Rad  = DEG2RAD * p.getLon();
      double cosD     = cos(dRad);
      double sinD     = sin(dRad);
      double sinLat1  = sin(lat1Rad);
      double cosLat1  = cos(lat1Rad);
      double cosH     = cos(hRad);
      double sinH     = sin(hRad);
      double latrad   = asin(sinLat1 * cosD + cosLat1 * sinD * cosH);
      double dlonrad  = atan2(sinH * sinD * cosLat1, cosD - sinLat1 * sin(latrad));
      double lonrad   = (lon1Rad + dlonrad + PI) % TWOPI - PI;
      double latdeg   = RAD2DEG * latrad;
      double londeg   = RAD2DEG * lonrad;
      londeg = (londeg > 180.0) ? londeg - 360.0 : (londeg < -180) ? londeg + 360.0 : londeg;
      return Position.fromGeo(0.0, latdeg, londeg);
    }
  }


  /** Returns distance from the nadir point directly below an aircraft/spacecraft at altitude
      <code>alt</code> to a viewer when the elevation <code>elev</code> of the aircraft/spacecraft relative
      the viewer is known.
      @param alt  Altitude of aircraft/spacecraft.
      @param elev Look angle from viewer to aircraft/spacecraft, must be between 0 and 90
                  (0=looking horizontal, 90=looking directly up).
      @return The distance in meters.
   */
  public static double computeDistance (double alt, double elev) { // new method
    double e = DEG2RAD * elev + PI/2;
    double r = EARTH_EQUATORIAL_RADIUS;
    double d = alt + r;
    double a = asin((sin(e) * EARTH_EQUATORIAL_RADIUS) / d);
    double b = PI - a - e;
    return EARTH_EQUATORIAL_RADIUS * b; //<-- dist in meters;
  }

  /** Compute the distance to the horizon of the Earth as seen from a given altitude. This is based
      on the following computation:
      <pre>
        Let P a point external to a circle (with center point O) and T be a point on the circle
        such that the line PT is tangent to the circle. By the definition of a tangent, PT and
        OT are perpendicular, therefore PTO is a right triangle and we can solve for the angle
        w that lies between PT and PO.

            sin(w) = r/d
            w = asin(r/d)

        The other interior angle of the triangle (90-w) is therefore the radius of the horizon
        line as measured in deg lat or lon.
      </pre>
      @param alt The altitude in meters.
      @return The horizon distance as measured in deg latitude.
   */
  public static double computeHorizonDistance (double alt) { // 100% taken from LayerVehicle.java
    double r = EARTH_EQUATORIAL_RADIUS;
    double d = alt+r;
    double w = asin(r/d); // returns angle in radians
    return 90.0 - toDegrees(w);
  }

  /** Computes a circle on the surface of the earth that approximates the horizon line for a given
      observer. This is merely a combination of calls to {@link #computeHorizonDistance} and
      {@link #computeCircleLL}.
      @param pos    The position of the observer.
      @param div    Number of divisions (line segments) to make a circle (-1=UseDefault).
      @return Array of {lon,lat} points that make up the horizon line. This will be a closed shape
              (i.e. first and last points are identical).
   */
  public static double[] computeHorizonLL (Position pos, int div) { // new method
    double radius = computeHorizonDistance(pos.getAlt());
    return computeCircleLL(pos.getLon(), pos.getLat(), radius, div);
  }

  /** Computes a circle on the surface of the earth that approximates the horizon line for a given
      observer. This is merely a combination of calls to {@link #computeHorizonDistance} and
      {@link #computeCircleLL}.
      @param lon    The longitude of the observer.
      @param lat    The latitude of the observer.
      @param alt    The altitude of the observer.
      @param div    Number of divisions (line segments) to make a circle (-1=UseDefault).
      @return Array of {lon,lat} points that make up the horizon line. This will be a closed shape
              (i.e. first and last points are identical).
   */
  public static double[] computeHorizonLL (double lon, double lat, double alt, int div) { // new method
    double radius = computeHorizonDistance(alt);
    return computeCircleLL(lon, lat, radius, div);
  }

  /** Computes a circle on the surface of the earth with a radius measured in deg latitude. This
      method explicitly takes into account the fact that longitude values "narrow" near the poles
      (i.e. 1deg at Equator = N meters, 1deg at Pole = 0 meters). After computing the line that
      approximates the circle.
      @param pos    The location of the center point.
      @param radius The radius of the circle measured in deg latitude.
      @param div    Number of divisions (line segments) to make a circle (-1=UseDefault).
      @return Array of {lon,lat} points that make up the limit of the circle. This will be a closed
              shape (i.e. first and last points are identical).
   */
  public static double[] computeCircleLL (Position pos, double radius, int div) { // new method
    return computeCircleLL(pos.getLon(), pos.getLat(), radius, div);
  }

  /** Computes a circle on the surface of the earth with a radius measured in deg latitude. This
      method explicitly takes into account the fact that longitude values "narrow" near the poles
      (i.e. 1deg at Equator = N meters, 1deg at Pole = 0 meters). After computing the line that
      approximates the circle.
      @param lon    The longitude of the center point.
      @param lat    The latitude of the center point.
      @param radius The radius of the circle measured in deg latitude.
      @param div    Number of divisions (line segments) to make a circle (-1=UseDefault).
      @return Array of {lon,lat} points that make up the limit of the circle. This will be a closed
              shape (i.e. first and last points are identical).
   */
  public static double[] computeCircleLL (double lon, double lat, double radius, int div) { // 90% taken from MPlot.java
    if (div == -1) { div = 90; } // JC4: We can optimize the default value in the future.

    int      len = (div+1)*2;
    double[] buf = new double[len]; // Points in the line { x1, y1, x2, y2, ... }
    double   ang = 0;               // Current angle to solve for.

/*
    if ((lat >= 90) || (lat <= -90)) {
      // Handle the special case of being directly over the N/S pole.
      double deg = 360 / div; // Use deg rather than radians

      for (int i = 0; (i < buf.length); i+=2, ang+=deg) {
        buf[i+0] = ang;
        buf[i+1] = (lat==90)? (90 - radius): (-90 + radius);
      }
*/
    if (lat >= 90) {
      // Handle the special case of being directly over the N pole.
      for (int ii=0, jj=0; (ii < buf.length); ii+=2, ++jj) {
        buf[ii] = 180.0 - jj*360.0/div;
                                if (buf[ii]<-180.0) buf[ii]+=360.0;
                                if (buf[ii]>+180.0) buf[ii]-=360.0;
        buf[ii+1] = 90 - radius;
      }
    } else if (lat <= -90) {
      // Handle the special case of being directly over the S pole.
      for (int ii=0, jj=0; (ii < buf.length); ii+=2, ++jj) {
        buf[ii] = 0.0 + jj*360.0/div;
                                if (buf[ii]<-180.0) buf[ii]+=360.0;
                                if (buf[ii]>+180.0) buf[ii]-=360.0;
        buf[ii+1] = radius - 90;
      }
    } else {
      // For speed reasons, do all computation in radians and the convert
      // to deg just before adding the value to the buffer.
      double w      = TWOPI / div;
      double latRad = toRadians(lat);
      double sinLat = sin(latRad);
      double cosLat = cos(latRad);
      double d      = toRadians(radius);
      double cosD   = cos(d);
      double sinD   = sin(d);

      for (int i = 0; (i < len); i+=2, ang+=w) {
        // The computation or xRad and yRad is loosely based on equations posted on
        // http://viavirtualearth.com/Wiki/Draw+a+circle.ashx These equations take into
        // account the "narrowing" of the longitudinal distance (in meters) near the poles.
        double yRad = asin(sinLat*cosD + cosLat*sinD*cos(ang));
        double xRad = atan2(sin(ang)*sinD*cosLat, cosD-sinLat*sin(yRad));

        buf[i+0] = normalizeLongitude(toDegrees(xRad) + lon); // added normalized since 3.9.0
        buf[i+1] = toDegrees(yRad);
      }
    }
    buf[len-2] = buf[0];
    buf[len-1] = buf[1];

    return buf;
  }

  /** Draw a conic (scaled according to view).
      This is the same as <code>computeConeLLA(p,q,angle,div,0,1)</code>.
      @param p     The position the cone is originating from (the viewer).
      @param q     The target location (the point the viewer is looking at).
      @param angle The cone angle (in deg).
      @param div   Number of divisions (line segments) to make a cone (-1=UseDefault).
      @return Array of {lon,lat,alt} points that make up the limit of the cone.
   */
  public static double[] computeConeLLA (Position p, Position q, double angle, int div) {
    return computeConeLLA(p,q,angle,div,0,1);
  }

  /** Draw a conic. Draws a cone on the earth by first drawing a circle on a plane perpendicular
      to a line from <code>p</code> (viewer) to <code>q</code> (target), then projecting that circle on
      to the surface of the earth.<br>
      <br>
      Based on nxm.sys.libg MPlot#drawCone, options for <code>cut</code> and <code>mode=-1</code> were removed
      to create a single continuous contour for each function call. Option <code>line</code> was removed as
      it is no longer needed. Also parameters were re-ordered to match other methods in this class.
      @param p     The position the cone is originating from (the viewer).
      @param q     The target location (the point the viewer is looking at - point on/towards earth).
      @param div   Number of divisions (line segments) to make a cone (-1=UseDefault).
      @param angle The half-angle of the cone (in deg).
      @param range The range, set this to 0 to compute range from <code>p</code> and <code>q</code>.
      @param mode  The clipping mode when clipping at the Earth's surface
                   (see {@link Transform#clipAtEarth}):<pre>
                   mode = 1 - Cone extends to horizon.
                   mode = 0 - Clip cone at a plane perpendicular to q.</pre>
      @return Array of {lon,lat,alt} points that make up the limit of the cone.
   */
  public static double[] computeConeLLA (Position p, Position q, double angle, int div,
                                         double range, int mode) {
    return computeConeLLA(p, q, angle, div, range, mode, false);
  }

  /** Computes a cone on the earth by first drawing a circle on a plane perpendicular
      to a line from <code>p</code> (viewer) to <code>q</code> (target), then projecting that circle on
      to the surface of the earth.<br>
      <br>
      Based on nxm.sys.libg MPlot#drawCone, options for <code>cut</code> and <code>mode=-1</code> were removed
      to create a single continuous contour for each function call. Option <code>line</code> was removed as
      it is no longer needed. Also parameters were re-ordered to match other methods in this class.
      @param p     The position the cone is originating from (the viewer).
      @param q     The target location (the point the viewer is looking at - point on/towards earth).
      @param div   Number of divisions (line segments) to make a cone (-1=UseDefault).
      @param angle The half-angle of the cone (in deg).
      @param range The range, set this to 0 to compute range from <code>p</code> and <code>q</code>.
      @param mode  The clipping mode when clipping at the Earth's surface
                   (see {@link Transform#clipAtEarth}):<pre>
                   mode = 1 - Cone extends to horizon.
                   mode = 0 - Clip cone at a plane perpendicular to q.</pre>
      @param trimToGeosurface trim the cone to the earth's surface
      @return Array of {lon,lat,alt} points that make up the limit of the cone.
*/
  @InternalUseOnly
  public static double[] computeConeLLA (Position p, Position q, double angle, int div,
      double range, int mode, boolean trimToGeosurface) {
    // taken from MPlot.java, minor mods as noted, plus mods for readability
    if (div == -1) { div = 90; } // JC4: We can optimize the default value in the future.

    int       npts  = div+1;               // number of points (+1 since last=first) - MPlot has div+2, but uses elen for alloc array
    double    rclip = EARTH_EQUATORIAL_RADIUS+10000.0;
    double[]  tdata = new double[3];       // target data (internal)
    double[]  vdata = new double[3];       // viewer data (internal)
    double[]  cidata = new double[3];      // cone cross section center data (internal)
    double[]  cdata = new double[npts*3];  // cone data - note this has only cone data, where in MPlot cidata & vdata are at beginning and end
    Transform ctp   = new Transform();     // Transform used when drawing cones (viewer)
    Transform ctq   = new Transform();     // Transform used when drawing cones (target)
    if (p.referenceFrame!=Position.GEODETIC) p.car2geo();
    if (q.referenceFrame!=Position.GEODETIC) q.car2geo();
    ctp.setTOP(p.alt,p.lat,p.lon);  // Set TOPOCENTRIC reference frame relative to viewing position
    ctq.setTOP(q.alt,q.lat,q.lon);  // Set TOPOCENTRIC reference frame relative to target position

    // Transform origin of target into coordinates of viewing frame of reference (tdata is {0,0,0} prior to this)
    Transform.top2top(tdata, Transform.POS, 1, ctq, ctp); // like MPlot, except results stored in data there and tdata here

    // Based on transformed coordinates, calculate r, roll and elev from viewer to target
    double x = tdata[0];
    double y = tdata[1];
    double z = tdata[2];
    double r = sqrt(x*x+y*y+z*z);
    double azim = 0;
    double elev = (x==0 && z==0)? 0 : RAD2DEG*atan2(-x,z);
    double roll = (r == 0)?       0 : RAD2DEG*asin(min(1.0,max(-1.0,-y/r)));

    // set reference frame of viewer so that x-y plane is perpendicular to line
    // from target to viewer, with large z values progressing towards the target
    ctp.setTOP(p.alt,p.lat,p.lon,azim,elev,roll);
    if (range>0) r = range;
    double dtheta = TWOPI/div;
    double cr = r*tan(DEG2RAD*angle); // cone cross section at distance r from viewer
    double cz = (range<0)? 0 : r; // z value of cone cross-section (height above viewer's topocentric x-y plane)
    cidata[2] = cz; ///in MPlot.drawCone this is at beginning of cdata, but then 1st 3 skipped during draw

    // draw circle of radius cr, with z axis cz
    int j = 0; // leave off cidata point from calcCone, where j starts loop as 3
    for (int i=0; i<div; i++) { // in MPlot this is < div, rather than <= div
      double theta = i*dtheta;
      cdata[j++] = cr*cos(theta);
      cdata[j++] = cr*sin(theta);
      cdata[j++] = cz;
    }

    // Transform from topocentric frame of reference to earth centered rotational
    vdata[2] = (range<0)? r : 0; // vdata[0]=vdata[1]=0 (by default)  // OPPOSITE of ciData[2]
    Transform.top2ecr(cidata,Transform.POS,   1, ctp); // matches cdata[0,1,2] in drawCone
    Transform.top2ecr(vdata, Transform.POS,   1, ctp); // same as last cdata points MPlot.drawCone
    Transform.top2ecr(cdata, Transform.POS, div, ctp); // matches cdata[3- in calcCone, but 0,1,2 are not drawn .. i.e. ignore them

    // Extract the transformed position to calculate distance from target to viewer (r) and check
    // to see whether cone cross section intersects earth, if so, clip the cross section at the
    // surface of the earth
    x=cidata[0]; y=cidata[1]; z=cidata[2]; r=sqrt(x*x+y*y+z*z);
    if (r<rclip) Transform.clipAtEarth(cdata, 0, npts, vdata, mode, trimToGeosurface, p, angle);

    // Convert to geodetic LLA for return to caller and copy first point to end so first=last
    Transform.car2geo(cdata, Transform.POS, div); // note: MPlot is div + 3 and this is just div (leaving off first set of points)
    Transform.all2lla(cdata, div); // note: MPlot is div + 3 and this is just div (leaving off first set of points)
    System.arraycopy(cdata, 0, cdata, div*3, 3); // copy first point to last
    return cdata;
  }

  /** Computes an ellipse on a topocentric plane centered at <code>lon</code>, and <code>lat</code>,
      then transformed to LL (ignoring altitude).<br>
      <br>
      Based on nxm.sys.libg MPlot#drawEllipseOnMap(double, double, double, double, double, double, int, nxm.sys.libg Line, java.awt.Graphics),
        point plotting code was removed, leaving only code to compute
      ellipse shape.  Option <code>line</code> was removed as it is no longer needed.  Also, altitude was
      forced to zero on the resulting ellipse.
      @param lon       Longitude in degrees.
      @param lat       Latitude in degrees.
      @param semiMajor Semi-major axis in meters.
      @param semiMinor Semi-minor axis in meters.
      @param tilt      Tilt of major axis in degrees (0=North).
      @param div       Number of divisions (line segments) to make an ellipse (-1=Use Default).
                       Should be a multiple of four.
      @return Array of {lon,lat} points that make up the limit of the ellipse.
      @since NeXtMidas 3.1.0
   */
  public static double[] computeMapEllipse (double lon, double lat, double semiMajor,
                        double semiMinor, double tilt, int div) { // taken from MPlot.java, minor mods as noted, plus mods for readability

    return computeEllipseLL(lon, lat, semiMajor, semiMinor, 90-tilt, div);
  }

  /** Computes an ellipse on a topocentric plane centered at <code>lon</code>, and <code>lat</code>,
      then transformed to LL (ignoring altitude).<br>
      <br>
      Based on nxm.sys.libg MPlot#drawEllipse, point plotting code was removed, leaving only code to compute
      ellipse shape.  Option <code>line</code> was removed as it is no longer needed.  Also, altitude was
      forced to zero on the resulting ellipse.
      @param lon       Longitude in degrees.
      @param lat       Latitude in degrees.
      @param semiMajor Semi-major axis in meters.
      @param semiMinor Semi-minor axis in meters.
      @param tilt      Tilt of major axis in degrees. <b>NOTE: this rotating is done based on Java's method
                       starting at positive X, moving towards positive Y. To convert to 0 degrees = North,
                       do (90-tilt) or use {@link #computeMapEllipse(double, double, double, double, double, int)}</b>
      @param div       Number of divisions (line segments) to make an ellipse (-1=Use Default).
                       Should be a multiple of four.
      @return Array of {lon,lat} points that make up the limit of the ellipse.
  */
  public static double[] computeEllipseLL (double lon, double lat, double semiMajor,
                          double semiMinor, double tilt, int div) { // taken from MPlot.java, minor mods as noted, plus mods for readability
    if (div == -1) { div = 90; } // JC4: We can optimize the default value in the future.
    int elen        = (div+1)*3;
    double[]  edata = new double[elen]; // Data points used when drawing an ellipse.
    Transform et    = new Transform();  // Transform used when drawing ellipse
    int i,j;
    double r,a2i,b2i,theta,dtheta,ctheta,stheta;

    // Compute topocentric frame of given point
    et.setTOP(0.0,lat,lon,tilt,0.0,0.0);

    //compute ellipse
    a2i = 1.0/(semiMajor*semiMajor);
    b2i = 1.0/(semiMinor*semiMinor);
    dtheta = TWOPI/div;
    for (i=j=0; i<div; i++) {
      theta=i*dtheta; ctheta=Math.cos(theta); stheta=Math.sin(theta);
      r = Math.sqrt(1.0/(a2i*ctheta*ctheta+b2i*stheta*stheta));
      edata[j++] = r*ctheta;
      edata[j++] = r*stheta;
      edata[j++] = 0;
    }

    //transform ellipse to LLA
    Transform.top2ecr (edata,1,div,et);
    Transform.car2geo (edata,1,div);
    Transform.all2lla (edata,div);

    // close circle
    for (i=0; i<3; i++) edata[j++]=edata[i];

    // convert to LL, then return
    return fromLLAToLL(edata);
  }

  /** Converts all altitudes less than 0.5mm to 0.0. This method is available since most of the
      methods in this class are prone to small rounding errors (typically in the 8'th decimal
      place). Although this is usually of little consequence, altitudes that are approximately
      0.0m often need to be reported as 0.0 (especially if the rounding error causes the altitude
      value to drop slightly below 0.0). <b>This should never be used in cases where negative
      altitudes are valid or expected.</b>
      <br>
      This method makes the changes in-place but returns input array to permits nesting of calls,
      for example:
      <pre>
        double[] points = fixAltitudeLLA(computeConeLLA(...));
      </pre>
      @param points The array of points to round in LLA order (rounding is done in-place).
      @return <code>points</code>
   */
  public static double[] fixAltitudeLLA (double[] points) {
    for (int i = 2; i < points.length; i+=3) {
      if (points[i] < 0.0005) points[i] = 0.0;
    }
    return points;
  }

  /** Converts all altitudes to zero.
      <br>
      This method makes the changes in-place but returns input array to permits nesting of calls,
      for example:
      <pre>
        double[] points = zeroAltitudeLLA(computeConeLLA(...));
      </pre>
      @param points The array of points in LLA order.
      @return <code>points</code>
   */
  public static double[] zeroAltitudeLLA (double[] points) {
    for (int i = 2; i < points.length; i+=3) {
      points[i] = 0.0;
    }
    return points;
  }

  /** Converts a double array of LLA to an array of LL.
      @param points The array of points in LLA order (must have length divisible by 3).
      @return <code>data</code>
   */
  public static double[] fromLLAToLL (double[] points) {
    double data[] = new double[points.length * 2 / 3];
    for (int i=0,j=0; i < points.length; i+=3) {
      data[j++] = points[i+0];
      data[j++] = points[i+1];
    }
    return data;
  }

  /** Converts a double array of LL {lon,lat} to an array of ALL {alt, lat, lon}.
      @param ll The array of points in LL order (must have length divisible by 2).
      @return <code>all</code>
      @since NeXtMidas 3.5.3 Bug 2874
  */
  public static double[] fromLLToALL (final double[] ll) {
    final int llLength = ll.length;
    final double all[] = new double[llLength * 3 / 2];
    for (int i=0,j=0; i < llLength; i+=2) {
      all[j++] = 0.0;     // alt
      all[j++] = ll[++i]; // lat
      all[j++] = ll[--i]; // lon
    }
    return all;
  }

  /** Gets the center of a bounding box given in ECEF/Geodetic coordinates
      @param perimeterECEFGeo bounding box perimeter
      @return coordinates of center of bounding box
      @since NeXtMidas 3.5.3  Bug 2874
   */
  public static double[] getCenterECEFGeo (final double[] perimeterECEFGeo) {
    // transform the AoI to approximately (0.0N 0.0E)
    final double[] crudeCenterEG = getCenterOfBoundingBoxECEFGeo(perimeterECEFGeo);
    TMatrix tm = new TMatrix();
    tm.rotate( 90.0-crudeCenterEG[2],TMatrix.AZIM);
    tm.rotate(     -crudeCenterEG[1],TMatrix.ROLL);
    tm.rotate(-90.0                 ,TMatrix.AZIM);
    final double[] perimeterEGO = Transform.ecefGeo2ecefGeo(tm, perimeterECEFGeo);

    // in this origin-centric reference frame, calculate a more-precise center estimate
    final double[] fineCenterEGO = getCenterOfBoundingBoxECEFGeo(perimeterEGO);

    // undo the previous transform and return the result
    tm.transpose();
    return Transform.ecefGeo2ecefGeo(tm, fineCenterEGO);
  }

  /** Alias for getCenterECEFGeo(double[])}
      @param perimeterECRGeo see linked method for docs
      @return  see linked method for docs
      @see #getCenterECEFGeo(double[])
      @since NeXtMidas 3.5.3  Bug 2874
   */
  public static double[] getCenterECRGeo (final double[] perimeterECRGeo) {
    return getCenterECEFGeo(perimeterECRGeo);
  }

  /** Gets a circle representing a circle of latitude of the Earth in ECR/Geodetic coordinates
      @param lat latitude of circle (-90.0 to +90.0)
      @param numSidesForShape number of sides for the shape
      @return points describing the equator
      @since NeXtMidas 3.5.3  Bug 2874
   */
  public static double[] getCircleOfLatitudeECEFGeo (final double lat, int numSidesForShape) {
    final double[] circle = new double[(numSidesForShape+1)*3];
    double lon = 0.0;
    final double lonInc = 360.0 / numSidesForShape;
    int ii = 0;
    for (int jj=0; jj<numSidesForShape; ++jj) {
      circle[ii++] = 0.0; // altitude
      circle[ii++] = lat; // latitude
      circle[ii++] = lon; // longitude
      lon += lonInc;
      if (lon>180.0) {
        lon -= 360.0;
      }
    }

    // repeat the first point
    circle[ii++] = 0.0; // altitude
    circle[ii++] = lat; // latitude
    circle[ii]   = 0.0; // longitude
    return circle;
  }

  /** Alias for getCircleOfLatitudeECEFGeo(double, int)
      @param lat latitude of circle (-90.0 to +90.0)
      @param numSidesForShape see linked method for docs
      @return see linked method for docs
      @see #getCircleOfLatitudeECEFGeo(double, int)
      @since NeXtMidas 3.5.3  Bug 2874
   */
  public static double[] getCircleOfLatitudeECRGeo (final double lat, int numSidesForShape) {
    return getCircleOfLatitudeECEFGeo(lat, numSidesForShape);
  }

  /** Gets crude center of a bounding box given in ECEF/Geodetic coordinates
      @param perimeterECEFGeo
      @return crude center of a bounding box
      @since NeXtMidas 3.5.3  Bug 2874
   */
  private static double[] getCenterOfBoundingBoxECEFGeo (final double[] perimeterECEFGeo) {
    final int LENGTH = (perimeterECEFGeo.length/3) * 3; // defensively round down if not given a multiple of 3
    // convert to ECR-CARTESIAN coordinates and find the center of a bounding box
    double xMax = -Double.MAX_VALUE;
    double xMin =  Double.MAX_VALUE;
    double yMax = -Double.MAX_VALUE;
    double yMin =  Double.MAX_VALUE;
    double zMax = -Double.MAX_VALUE;
    double zMin =  Double.MAX_VALUE;
    for (int ii=0; ii<LENGTH;) {
      // analogous to Transform.geo2car()
      final double alt  = perimeterECEFGeo[ii++];
      final double lat  = Math.toRadians(perimeterECEFGeo[ii++]); // lat
      final double lon  = Math.toRadians(perimeterECEFGeo[ii++]); // lon
      final double clat = Math.cos(lat);
      final double x    = Math.cos(lon)*clat;
      final double y    = Math.sin(lon)*clat;
      final double z    = Math.sin(lat);

      xMax = Math.max(x, xMax);
      xMin = Math.min(x, xMin);
      yMax = Math.max(y, yMax);
      yMin = Math.min(y, yMin);
      zMax = Math.max(z, zMax);
      zMin = Math.min(z, zMin);
    }
    final double xCenter = (xMax+xMin)/2.0;
    final double yCenter = (yMax+yMin)/2.0;
    final double zCenter = (zMax+zMin)/2.0;

    // convert the center back to ECR-GEODETIC coordinates and force the altitude to zero
    // analogous to Transform.car2geo()
    final double[] centerECRGeo = new double[3];
    centerECRGeo[0] = 0.0;
    if (zCenter>=1.0) { // preclude asin() domain errors
      centerECRGeo[1] = 90.0; // lat
      centerECRGeo[2] = 0.0; // lon
    } else if (zCenter<=-1.0) { // preclude asin() domain errors
      centerECRGeo[1] = -90.0; // lat
      centerECRGeo[2] = 0.0; // lon
    } else if (xCenter==0.0 && yCenter==0.0) { // preclude atan2() domain errors
      centerECRGeo[1] = (zCenter>0.0) ? 90.0 : -90.0; // lat
      centerECRGeo[2] = 0.0; // lon
    } else {
      centerECRGeo[1] = Math.toDegrees(Math.asin(zCenter)); // lat
      centerECRGeo[2] = Math.toDegrees(Math.atan2(yCenter,xCenter)); // lon
    }
    return centerECRGeo;
  }


  /** Converts an array of points to a list of positions.
   *  @param format The format (not null).
   *  @param points The raw points values (not null).
   *  @return The positions.
   *  @since NeXtMidas 3.9.0
   */
  @ProvisionalUseOnly("Not adequately tested")
  public static List<Position> toPositions (ArrayForm format, double[] points) {
    List<Position> pts;

    switch (format) {
      case LonLatAlt:
        pts = new ArrayList<>((points.length / 3) + 1);
        for (int i = 0; i < points.length; i+=3) {
          pts.add(fromGeo(points[i+2], points[i+1], points[i]));
        }
        return pts;
      case LonLat:
        pts = new ArrayList<>((points.length / 2) + 1);
        for (int i = 0; i < points.length; i+=2) {
          pts.add(fromGeo(0.0, points[i+1], points[i]));
        }
        return pts;
      case LatLonAlt:
        pts = new ArrayList<>((points.length / 3) + 1);
        for (int i = 0; i < points.length; i+=3) {
          pts.add(fromGeo(points[i+2], points[i], points[i+1]));
        }
        return pts;
      case LatLon:
        pts = new ArrayList<>((points.length / 2) + 1);
        for (int i = 0; i < points.length; i+=2) {
          pts.add(fromGeo(0.0, points[i], points[i+1]));
        }
        return pts;
      case AltLatLon:
        pts = new ArrayList<>((points.length / 3) + 1);
        for (int i = 0; i < points.length; i+=3) {
          pts.add(fromGeo(points[i], points[i+1], points[i+2]));
        }
        return pts;
      default:
        throw new AssertionError("Invalid format: "+format);
    }
  }

  /** Converts a list of positions to an array of points.
   *  @param format The format (not null).
   *  @param points The raw points values (not null).
   *  @return The positions.
   *  @since NeXtMidas 3.9.0
   */
  @ProvisionalUseOnly("Not adequately tested")
  public static double[] toPoints (ArrayForm format, List<Position> points) {
    double[] pts;
    int j = 0;

    switch (format) {
      case LonLatAlt:
        pts = new double[points.size()*3];
        for (Position p : points) {
          pts[j++] = p.getLon();
          pts[j++] = p.getLat();
          pts[j++] = p.getAlt();
        }
        return pts;
      case LonLat:
        pts = new double[points.size()*2];
        for (Position p : points) {
          pts[j++] = p.getLon();
          pts[j++] = p.getLat();
        }
        return pts;
      case LatLonAlt:
        pts = new double[points.size()*3];
        for (Position p : points) {
          pts[j++] = p.getLat();
          pts[j++] = p.getLon();
          pts[j++] = p.getAlt();
        }
        return pts;
      case LatLon:
        pts = new double[points.size()*2];
        for (Position p : points) {
          pts[j++] = p.getLat();
          pts[j++] = p.getLon();
        }
        return pts;
      case AltLatLon:
        pts = new double[points.size()*3];
        for (Position p : points) {
          pts[j++] = p.getAlt();
          pts[j++] = p.getLat();
          pts[j++] = p.getLon();
        }
        return pts;
      default:
        throw new AssertionError("Invalid format: "+format);
    }
  }

  /** Closes the polygon if it is not already closed.
   *  @param pts The points.
   *  @return {@code pts} (this modified the parameter passed in).
   *  @throws NullPointerException     If the list is null.
   *  @throws IllegalArgumentException If the list has fewer than 3 distinct points.
   *  @since NeXtMidas 3.9.0
   */
  @ProvisionalUseOnly("Not adequately tested")
  public static List<Position> closeShape (List<Position> pts) {
    if (pts.size() < 3) {
      throw new IllegalArgumentException("Can not create polygon with less than 3 distinct points");
    }

    Position first = pts.get(0);
    Position last  = pts.get(pts.size()-1);

    if ((first.getLat() != last.getLat()) ||
        (first.getLon() != last.getLon()) ||
        (first.getAlt() != last.getAlt())) {
      // close the polygon
      pts.add(first);
    }

    if (pts.size() < 4) {
      throw new IllegalArgumentException("Can not create polygon with less than 3 distinct points");
    }
    return pts;
  }

  /** Computes an ellipse on a topocentric plane. <br>
   *  <br>
   *  WARNING: This uses tilt measured from <b>true north</b> which matches
   *           {@code computeMapEllipse(..)} but differs from {@code computeEllipseLL(..)}.
   *  @param lon       Longitude in degrees.
   *  @param lat       Latitude in degrees.
   *  @param semiMajor Semi-major axis in meters.
   *  @param semiMinor Semi-minor axis in meters.
   *  @param tilt      Tilt of major axis in degrees (0=North).
   *  @param div       Number of divisions (line segments) to make an ellipse (-1=Use Default).
   *                   Should be a multiple of four.
   *  @return The positions ordered clockwise as a closed shape (first and last points match).
   *  @since NeXtMidas 3.9.0
   */
  @ProvisionalUseOnly("Not adequately tested")
  public static List<Position> computeEllipse (double lon, double lat, double semiMajor,
                                               double semiMinor, double tilt, int div) {
    double[]       points    = computeEllipseLL(lon, lat, semiMajor, semiMinor, 90-tilt, div);
    List<Position> positions = toPositions(ArrayForm.LonLat, points);
    reverse(positions);
    return closeShape(positions);
  }

  /** Computes a circle on the surface of the earth with a radius measured in meters.  <br>
   *  <br>
   *  WARNING: This uses <b>meters</b> rather than degrees which matches
   *           {@code computeEllipse(..)} but differs from {@code computeCircleLL(..)}.
   *  @param pos    The location of the center point.
   *  @param radius The radius of the circle in meters.
   *  @param div    Number of divisions (line segments) to make a circle (-1=UseDefault).
   *                Should be a multiple of four.
   *  @return The positions ordered clockwise as a closed shape (first and last points match).
   *  @since NeXtMidas 3.9.0
   */
  @ProvisionalUseOnly("Not adequately tested")
  public static List<Position> computeCircle (Position pos, double radius, int div) {
    return computeEllipse(pos.getLon(), pos.getLat(), radius, radius, 0, div);
  }

  /** Computes a circle on the surface of the earth with a radius measured in meters.  <br>
   *  <br>
   *  WARNING: This uses <b>meters</b> rather than degrees which matches
   *           {@code computeEllipse(..)} but differs from {@code computeCircleLL(..)}.
   *  @param lon    The longitude of the center point.
   *  @param lat    The latitude of the center point.
   *  @param radius The radius of the circle in meters.
   *  @param div    Number of divisions (line segments) to make a circle (-1=UseDefault).
   *                Should be a multiple of four.
   *  @return The positions ordered clockwise (first and last points match).
   *  @since NeXtMidas 3.9.0
   */
  @ProvisionalUseOnly("Not adequately tested")
  public static List<Position> computeCircle (double lon, double lat, double radius, int div) {
    return computeEllipse(lon, lat, radius, radius, 0, div);
  }

  /** Draw a conic.
   *  This is the same as <code>computeConeLLA(p,q,angle,div,0,1)</code>.
   *  @param p     The position the cone is originating from (the viewer).
   *  @param q     The target location (the point the viewer is looking at).
   *  @param angle The cone angle (in deg).
   *  @param div   Number of divisions (line segments) to make a cone (-1=UseDefault).
   *  @return The positions ordered clockwise as a closed shape (first and last points match).
   *  @since NeXtMidas 3.9.0
   */
  @ProvisionalUseOnly("Not adequately tested") /// Jennifer calling this
  public static List<Position> computeCone (Position p, Position q, double angle, int div) {
    return computeCone(p,q,angle,div,0,1);
  }

  /** Draw a conic. Draws a cone on the earth by first drawing a circle on a plane perpendicular
   *  to a line from <code>p</code> (viewer) to <code>q</code> (target), then projecting that circle on
   *  to the surface of the earth.<br>
   *  <br>
   *  Based on nxm.sys.libg MPlot#drawCone, options for <code>cut</code> and <code>mode=-1</code> were removed
   *  to create a single continuous contour for each function call. Option <code>line</code> was removed as
   *  it is no longer needed. Also parameters were re-ordered to match other methods in this class.
   *  @param p     The position the cone is originating from (the viewer).
   *  @param q     The target location (the point the viewer is looking at).
   *  @param div   Number of divisions (line segments) to make a cone (-1=UseDefault).
   *  @param angle The half-angle of the cone (in deg).
   *  @param range The range, set this to 0 to compute range from <code>p</code> and <code>q</code>.
   *  @param mode  The clipping mode when clipping at the Earth's surface
   *               (see {@link Transform#clipAtEarth}):<pre>
   *               mode = 1 - Cone extends to horizon.
   *               mode = 0 - Clip cone at a plane perpendicular to q.</pre>
   *  @return The positions ordered clockwise as a closed shape (first and last points match).
   *  @since NeXtMidas 3.9.0
   */
  @ProvisionalUseOnly("Not adequately tested")
  public static List<Position> computeCone (Position p, Position q, double angle, int div,
                                            double range, int mode) {
    double[]       points    = computeConeLLA(p, q, angle, div, range, mode,true);
    List<Position> positions = toPositions(ArrayForm.LonLatAlt, points);
    return closeShape(positions);
  }

  /** Get the approximate 3 dB beamwidth for a parabolic dish antenna. This is here since it
   *  complements {@code computeCone(..)}. Note that when the value must be divided by two to
   *  get the half-angle before passing to {@code computeCone(..)}. Example:
   *  <pre>
   *    computeCone(p, q, getBeamWidth(d,rf)/2, div)
   *  </pre>
   *  @param d   Dish diameter in meters.
   *  @param rf  Signal frequency in Hz.
   *  @return Beamwidth in deg.
   *  @since NeXtMidas 3.9.0
   */
  @ProvisionalUseOnly("Not adequately tested")
  public static double getBeamWidth (double d, double rf) {
    return (58 * 1.11 * SPEED_OF_LIGHT) / rf / d;
  }
  /** Gets the countries shape file being used.
   *  @return The shape file.
   *  @since NeXtMidas 3.9.0
   */
  @ProvisionalUseOnly("Not adequately tested")
  public static synchronized ShapeFile getCountryShapeFile () {
    if (countries == null) {
      countries = new ShapeFile(null, "classpath:nxm/sys/dat/world.shp");
      countries.open(ShapeFile.INPUT);
    }
    return countries;
  }

  /** Sets the countries shape file being used. Note that this will not work well
   *  with files other than {@code classpath:nxm/sys/dat/world.shp} or (the more
   *  modern)  {@code classpath:nxm/map/dat/world.shp}.
   *  @param ref      The Midas reference to use.
   *  @param fileName The name of the shape file.
   *  @since NeXtMidas 3.9.0
   */
  @ProvisionalUseOnly("Not adequately tested")
  public static void setCountryShapeFile (MidasReference ref, Object fileName) {
    setCountryShapeFile(new ShapeFile(ref, fileName));
  }

  /** Sets the countries shape file being used. Note that this will not work well
   *  with files other than {@code classpath:nxm/sys/dat/world.shp} or (the more
   *  modern)  {@code classpath:nxm/map/dat/world.shp}.
   *  @param f The shape file.
   *  @since NeXtMidas 3.9.0
   */
  @ProvisionalUseOnly("Not adequately tested")
  public static synchronized void setCountryShapeFile (ShapeFile f) {
    if ((f != null) && !f.isOpen()) {
      f.open(ShapeFile.INPUT);
    }
    countries = f;
  }

  /** Gets information about a country that is at the given location.
   *  @param pos    The position.
   *  @param radius The *approximate* "nearness" radius in meters. If 0.0 only inclusiveness will
   *                be checked. If -1.0 a reasonable default will be used.
   *  @return The information from the DBF file accompanying {@link #getCountryShapeFile()} for
   *          each potential match. Priority is given for inclusion over nearness. This may be
   *          empty but will not be null.
   *  @since NeXtMidas 3.9.0
   */
  private static List<ShapeFile.Shape> getCountryShape (Position pos, double radius) {
    if (pos.getAlt() >= 100_000) {
      // Outer space...
      return Collections.emptyList();
    }

    List<ShapeFile.Shape> shapes = getCountryShapeFile().getShapesContaining(pos);

    if (radius < 0) {
      // Use nominal distance for "territorial waters" as a default
      radius = 22_200; // 22.2 km = 12 nmi
    }

    if (shapes.isEmpty() && (radius > 0)) {
      shapes = getCountryShapeFile().getShapesNear(pos, radius);
    }

    return shapes;
  }

  /** Gets information about a country that is at the given location.
   *  @param alt    The altitude.
   *  @param lat    The latitude.
   *  @param lon    The longitude.
   *  @return The information from the DBF file accompanying {@link #getCountryShapeFile()} for
   *          each potential match. Priority is given for inclusion over nearness. This may be
   *          empty but will not be null.
   *  @since NeXtMidas 3.9.0
   */
  public static List<Table> getCountryInfo (double alt, double lat, double lon) {
    return getCountryInfo(fromGeo(alt, lat, lon), -1);
  }

  /** Gets information about a country that is at the given location.
   *  @param pos    The position.
   *  @return The information from the DBF file accompanying {@link #getCountryShapeFile()} for
   *          each potential match. Priority is given for inclusion over nearness. This may be
   *          empty but will not be null.
   *  @since NeXtMidas 3.9.0
   */
  public static List<Table> getCountryInfo (Position pos) {
    return getCountryInfo(pos, -1);
  }

  /** Gets information about a country that is at the given location.
   *  @param alt    The altitude.
   *  @param lat    The latitude.
   *  @param lon    The longitude.
   *  @param radius The *approximate* "nearness" radius in meters. If 0.0 only inclusiveness will
   *                be checked. If -1.0 a reasonable default will be used.
   *  @return The information from the DBF file accompanying {@link #getCountryShapeFile()} for
   *          each potential match. Priority is given for inclusion over nearness. This may be
   *          empty but will not be null.
   *  @since NeXtMidas 3.9.0
   */
  public static List<Table> getCountryInfo (double alt, double lat, double lon, double radius) {
    return getCountryInfo(fromGeo(alt, lat, lon), radius);
  }

  /** Gets information about a country that is at the given location.
   *  @param pos    The position.
   *  @param radius The *approximate* "nearness" radius in meters. If 0.0 only inclusiveness will
   *                be checked. If -1.0 a reasonable default will be used.
   *  @return The information from the DBF file accompanying {@link #getCountryShapeFile()} for
   *          each potential match. Priority is given for inclusion over nearness. This may be
   *          empty but will not be null.
   *  @since NeXtMidas 3.9.0
   */
  public static List<Table> getCountryInfo (Position pos, double radius) {
    List<ShapeFile.Shape> shapes = getCountryShape(pos,radius);
    List<Table>           info   = new ArrayList<>(8);

    if (shapes != null) {
      shapes.forEach((s) -> {
        if (s != null) {
          info.add(s.getDataTable());
        }
      });
    }
    return info;
  }

  /** Gets information about a country that is at the given location. This can be used for
   *  country NAME or CODE lookup as follows:
   *  <pre>
   *    getCountryInfo("NAME", ...);
   *    getCountryInfo("CODE", ...);
   *  </pre>
   *  @param field  The name of the field (not null).
   *  @param alt    The altitude.
   *  @param lat    The latitude.
   *  @param lon    The longitude.
   *  @return The country code information from the DBF file accompanying
   *          {@link #getCountryShapeFile()} for each potential match. Priority is given
   *          for inclusion over nearness. This may be empty but will not be null.
   *  @since NeXtMidas 3.9.0
   */
  public static List<String> getCountryInfo (String field, double alt, double lat, double lon) {
    return getCountryInfo(field, fromGeo(alt, lat, lon), -1);
  }

  /** Gets information about a country that is at the given location. This can be used for
   *  country NAME or CODE lookup as follows:
   *  <pre>
   *    getCountryInfo("NAME", ...);
   *    getCountryInfo("CODE", ...);
   *  </pre>
   *  @param field  The name of the field (not null).
   *  @param pos    The position.
   *  @return The country code information from the DBF file accompanying
   *          {@link #getCountryShapeFile()} for each potential match. Priority is given
   *          for inclusion over nearness. This may be empty but will not be null.
   *  @since NeXtMidas 3.9.0
   */
  public static List<String> getCountryInfo (String field, Position pos) {
    return getCountryInfo(field, pos, -1);
  }

  /** Gets information about a country that is at the given location. This can be used for
   *  country NAME or CODE lookup as follows:
   *  <pre>
   *    getCountryInfo("NAME", ...);
   *    getCountryInfo("CODE", ...);
   *  </pre>
   *  @param field  The name of the field (not null).
   *  @param alt    The altitude.
   *  @param lat    The latitude.
   *  @param lon    The longitude.
   *  @param radius The *approximate* "nearness" radius in meters. If 0.0 only inclusiveness will
   *                be checked. If -1.0 a reasonable default will be used.
   *  @return The country code information from the DBF file accompanying
   *          {@link #getCountryShapeFile()} for each potential match. Priority is given
   *          for inclusion over nearness. This may be empty but will not be null.
   *  @since NeXtMidas 3.9.0
   */
  public static List<String> getCountryInfo (String field, double alt, double lat, double lon, double radius) {
    return getCountryInfo(field, fromGeo(alt, lat, lon), radius);
  }

  /** Gets information about a country that is at the given location. This can be used for
   *  country NAME or CODE lookup as follows:
   *  <pre>
   *    getCountryInfo("NAME", ...);
   *    getCountryInfo("CODE", ...);
   *  </pre>
   *  @param field  The name of the field (not null).
   *  @param pos    The position.
   *  @param radius The *approximate* "nearness" radius in meters. If 0.0 only inclusiveness will
   *                be checked. If -1.0 a reasonable default will be used.
   *  @return The country code information from the DBF file accompanying
   *          {@link #getCountryShapeFile()} for each potential match. Priority is given
   *          for inclusion over nearness. This may be empty but will not be null.
   *  @since NeXtMidas 3.9.0
   */
  public static List<String> getCountryInfo (String field, Position pos, double radius) {
    return getCountryShape(pos,radius).stream().map((s) -> s.getDataTable().getS(field)).collect(toList());
  }

  /** Calculates radius of CEP95 circle from EEP95 ellipse parameters. The computation is not
   *  specific to the units used, so the units output will match the units passed in (e.g.
   *  meters in yields meters out).
   *  @param sig_smaj_95 EEP semi-major axis.
   *  @param sig_smin_95 EEP semi-minor axis.
   *  @return The CEP.
   *  @since NeXtMidas 3.9.0
   */
  @ProvisionalUseOnly("Not adequately tested")
  public static double eep2cep (double sig_smaj_95, double sig_smin_95) {
    double sig_smaj_245 = sig_smaj_95 / 2.4477468306808166; // 2.45
    double sig_smin_245 = sig_smin_95 / 2.4477468306808166; // 2.45
    double m            = pow(sig_smaj_245, 2) + pow(sig_smin_245, 2);
    double mm9          = 9*m*m; // = 9 * m**2
    double v            = 2*(pow(sig_smaj_245, 4) + pow(sig_smin_245, 4));
    double cep          = sqrt(m * pow(1.6447*sqrt(v/mm9) + (1 - v/mm9), 3));

    return (Double.isNaN(cep))? 0.0 : cep;
  }

  /** Create 'wire cage' by making cross sections down each axis
   *
   *  @param semiMajor   ellipsoid semi along X axis
   *  @param semiMinor   ellipsoid semi along Y axis
   *  @param semiMedian  ellipsoid semi along Z axis
   *  @param slices      number of slices (i.e. ellipse cross sections) to make along each axis
   *  @return the distance along the discretized axis and the semis for all the wire cage ellipses
   */
  @InternalUseOnly
  public static double[][][] computeXYZEllipsoidCrossSections (double semiMajor,
      double semiMinor, double semiMedian, int slices) {
    double sliceJump = 0.4;
    int axes = 3;
    int crossSecDirections = 3; // X, Y, Z
    int x=0, y=1, z=2;
    double[] semi = {semiMajor, semiMinor, semiMedian};
    double[][][] xyzCrossSections = new double[crossSecDirections][slices][axes];
    // first slice each direction is full size
    xyzCrossSections[0][0] = new double[] {0,semiMinor,semiMedian};
    xyzCrossSections[1][0] = new double[] {semiMajor,0,semiMedian};
    xyzCrossSections[2][0] = new double[] {semiMajor,semiMinor,0};
    for (int direction = 0; direction < crossSecDirections; direction++) {
      for (int slice = 1, index=1; slice < (slices/2+1); slice++,index+=2) {
        xyzCrossSections[direction][index] = new double[axes];
        xyzCrossSections[direction][index+1] = new double[axes];
        double ratio = slice*sliceJump;
        double crossLngth = semi[direction] * ratio;
        for (int axis = 0; axis < axes; axis++) {
          if (axis == direction) {
            xyzCrossSections[direction][index][axis] = crossLngth;
            xyzCrossSections[direction][index+1][axis] = -crossLngth;
          }
          else {
            xyzCrossSections[direction][index][axis] =
                semi[axis] * Math.sqrt((1-Math.pow(ratio, 2)));
            xyzCrossSections[direction][index+1][axis] =
                semi[axis] * Math.sqrt((1-Math.pow(ratio, 2)));
          }
        }
      }
    }
    return xyzCrossSections;
  }

  /** The rotation matrix, given the yaw, ptich, and roll
   *
   *  @param zRotation  rotation about Z axis in radians
   *  @param yRotation  rotation about Y axis in radians
   *  @param xRotation  rotation about X axis in radians
   *  @return rotation matrix
   */
  @InternalUseOnly
  public static double[][] getZYXRotationTMatrix (double zRotation, double yRotation, double xRotation) {
    double[][] zyxRotationMatrix = new double[3][3];
    // first row
    zyxRotationMatrix[0][0] = Math.cos(zRotation) * Math.cos(yRotation);
    zyxRotationMatrix[0][1] = Math.cos(zRotation) * Math.sin(yRotation) * Math.sin(xRotation) -
                      Math.sin(zRotation) * Math.cos(xRotation);
    zyxRotationMatrix[0][2] = Math.cos(zRotation) * Math.sin(yRotation) * Math.cos(xRotation) +
                      Math.sin(zRotation) * Math.sin(xRotation);
    // second row
    zyxRotationMatrix[1][0] = Math.sin(zRotation) * Math.cos(yRotation);
    zyxRotationMatrix[1][1] = Math.sin(zRotation) * Math.sin(yRotation) * Math.sin(xRotation) +
                      Math.cos(zRotation) * Math.cos(xRotation);
    zyxRotationMatrix[1][2] = Math.sin(zRotation) * Math.sin(yRotation) * Math.cos(xRotation) -
                      Math.cos(zRotation) * Math.sin(xRotation);
    // third row
    zyxRotationMatrix[2][0] = -Math.sin(yRotation);
    zyxRotationMatrix[2][1] = Math.cos(yRotation) * Math.sin(xRotation);
    zyxRotationMatrix[2][2] = Math.cos(yRotation) * Math.cos(xRotation);
    return zyxRotationMatrix;
  }

  /**  Use the insidePolygon fix that allows it to correctly detect pole crossings
   *
   *   @param useFix false to turn off the fix
   */
  public static void setinsidePolyFix (boolean useFix) {
    GeodeticUtil.insidePolyFix = useFix;
  }

}
