/* Contributed in the public domain.
    * Licensed to CS GROUP (CS) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * CS licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *   http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package org.orekit.models.earth;
  
  import static org.hamcrest.CoreMatchers.is;
  import static org.hamcrest.CoreMatchers.nullValue;
  import static org.hamcrest.CoreMatchers.sameInstance;
  import static org.hamcrest.MatcherAssert.assertThat;
  import static org.junit.Assert.assertEquals;
  import static org.orekit.OrekitMatchers.closeTo;
  import static org.orekit.OrekitMatchers.geodeticPointCloseTo;
  import static org.orekit.OrekitMatchers.vectorCloseTo;
  
  import java.util.Arrays;
  
  import org.hipparchus.geometry.euclidean.threed.Line;
  import org.hipparchus.geometry.euclidean.threed.Rotation;
  import org.hipparchus.geometry.euclidean.threed.Vector3D;
  import org.hipparchus.util.FastMath;
  import org.junit.Assert;
  import org.junit.Before;
  import org.junit.BeforeClass;
  import org.junit.Test;
  import org.orekit.Utils;
  import org.orekit.bodies.GeodeticPoint;
  import org.orekit.forces.gravity.potential.EGMFormatReader;
  import org.orekit.forces.gravity.potential.GravityFieldFactory;
  import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
  import org.orekit.frames.Frame;
  import org.orekit.frames.FramesFactory;
  import org.orekit.frames.Transform;
  import org.orekit.time.AbsoluteDate;
  
  /**
   * Unit tests for {@link Geoid}.
   *
   * @author Evan Ward
   */
  public class GeoidTest {
  
      /** maximum degree used in testing {@link Geoid}. */
      private static final int maxOrder = 360;
      /** maximum order used in testing {@link Geoid}. */
      private static final int maxDegree = 360;
      /** The WGS84 reference ellipsoid. */
      private static ReferenceEllipsoid WGS84 = new ReferenceEllipsoid(
              6378137.00, 1 / 298.257223563, FramesFactory.getGCRF(),
              3.986004418e14, 7292115e-11);
      /**
       * The potential to use in {@link #getComponent()}. Set in {@link
       * #setUpBefore()}.
       */
      private static NormalizedSphericalHarmonicsProvider potential;
      /** date to use in test cases */
      private static AbsoluteDate date;
  
      /**
       * load orekit data and gravity field.
       *
       * @throws Exception on error.
       */
      @BeforeClass
      public static void setUpBefore() throws Exception {
          Utils.setDataRoot("geoid:regular-data");
          GravityFieldFactory.clearPotentialCoefficientsReaders();
          GravityFieldFactory.addPotentialCoefficientsReader(
                  new EGMFormatReader("egm96", false));
          potential = GravityFieldFactory.getConstantNormalizedProvider(
                  maxDegree, maxOrder);
          date = null;
      }
  
      /** {lat, lon, expectedValue} points to evaluate the undulation */
      private double[][] reference;
  
      /** create the array of reference points */
      @Before
      public void setUp() {
          reference = new double[][]{
                  {0, 75, -100.3168},
                  {-30, 60, 14.9214},
                  {0, 130, 76.6053},
                  {60, -30, 63.7979},
                  {40, -75, -34.0402},
                  {28, 92, -30.4056},// the Himalayas
                 {45, 250, -7.4825},// the rockies
                 // this section is taken from the NGA's test file
                 {38.6281550, 269.7791550, -31.628},
                 {-14.6212170, 305.0211140, -2.969},
                 {46.8743190, 102.4487290, -43.575},
                 {-23.6174460, 133.8747120, 15.871},
                 {38.6254730, 359.9995000, 50.066},
                 {-.4667440, .0023000, 17.329}};
     }
 
     /**
      * Gets a new instance of {@link Geoid} to test with. It is given the EGM96
      * potential and the WGS84 ellipsoid.
      *
      * @return a new {@link Geoid}
      */
     private Geoid getComponent() {
         return new Geoid(potential, WGS84);
     }
 
     /** Test constructor and simple getters. */
     @Test
     public void testGeoid() {
         Geoid geoid = getComponent();
         // reference ellipse is the same
         assertEquals(WGS84, geoid.getEllipsoid());
         // geoid and reference ellipse are in the same frame
         assertEquals(WGS84.getBodyFrame(), geoid.getBodyFrame());
     }
 
     /** throws on null */
     @Test(expected = NullPointerException.class)
     public void testGeoidNullPotential() {
         new Geoid(null, WGS84);
     }
 
     /** throws on null */
     @Test(expected = NullPointerException.class)
     public void testGeoidNullEllipsoid() {
         new Geoid(potential, null);
     }
 
     /**
      * Test several pre-computed points from the Online Geoid Height Evaluation
      * tool, which takes into account terrain.
      *
      * @see <a href="http://geographiclib.sourceforge.net/cgi-bin/GeoidEval">Online
      * Geoid Height Evaluation</a>
      * @see <a href="http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html">Geoid
      * height for WGS84 and EGM96</a>
      */
     @Test
     public void testGetUndulation() {
         /*
          * allow 3 meter of error, which is what the approximations would
          * suggest, see the comment for Geoid.
          */
         final double maxError = 3;
 
         // run the test
         Geoid geoid = getComponent();
         for (double[] row : reference) {
             double lat = row[0];
             double lon = row[1];
             double undulation = geoid.getUndulation(FastMath.toRadians(lat),
                     FastMath.toRadians(lon), date);
             double expected = row[2];
             // System.out.format("%10g %10g %10g %10g%n", lat, lon, expected,
             // undulation - expected);
             Assert.assertEquals(String.format("lat: %5g, lon: %5g", lat, lon),
                     undulation, expected, maxError);
         }
     }
 
     /**
      * check {@link Geoid#getIntersectionPoint(Line, Vector3D, Frame,
      * AbsoluteDate)} with several points.
      */
     @Test
     public void testGetIntersectionPoint() {
         // setup
         Geoid geoid = getComponent();
         Frame frame = geoid.getBodyFrame();
 
         for (double[] point : reference) {
             GeodeticPoint gp = new GeodeticPoint(FastMath.toRadians(point[0]),
                     FastMath.toRadians(point[1]), 0);
             Vector3D expected = geoid.transform(gp);
             // glancing line: 10% vertical and 90% north (~6 deg elevation)
             Vector3D slope = gp.getZenith().scalarMultiply(0.1)
                     .add(gp.getNorth().scalarMultiply(0.9));
             Vector3D close = expected.add(slope.scalarMultiply(100e3));
             Vector3D pointOnLine = expected.add(slope);
             Line line = new Line(close, pointOnLine, 0);
             // line directed the other way
             Line otherDirection = new Line(pointOnLine, close, 0);
 
             // action
             GeodeticPoint actual = geoid.getIntersectionPoint(line, close,
                     frame, date);
             // other direction
             GeodeticPoint actualReversed = geoid.getIntersectionPoint(
                     otherDirection, close, frame, date);
 
             // verify
             String message = String.format("point: %s%n",
                     Arrays.toString(point));
             // position accuracy on Earth's surface to 1.3 um.
             assertThat(message, actualReversed, geodeticPointCloseTo(gp, 1.3e-6));
             assertThat(message, actual, geodeticPointCloseTo(gp, 1.3e-6));
         }
     }
 
     /**
      * check {@link Geoid#getIntersectionPoint(Line, Vector3D, Frame,
      * AbsoluteDate)} handles frame transformations correctly
      */
     @Test
     public void testGetIntersectionPointFrame() {
         // setup
         Geoid geoid = getComponent();
         Frame frame = new Frame(
                 geoid.getBodyFrame(),
                 new Transform(
                         date,
                         new Transform(
                                 date,
                                 new Vector3D(-1, 2, -3),
                                 new Vector3D(4, -5, 6)),
                         new Transform(
                                 date,
                                 new Rotation(-7, 8, -9, 10, true),
                                 new Vector3D(-11, 12, -13))),
                 "test frame");
         GeodeticPoint gp = new GeodeticPoint(FastMath.toRadians(46.8743190),
                 FastMath.toRadians(102.4487290), 0);
         Vector3D expected = geoid.transform(gp);
         // glancing line: 10% vertical and 90% north (~6 deg elevation)
         Vector3D slope = gp.getZenith().scalarMultiply(0.1)
                 .add(gp.getNorth().scalarMultiply(0.9));
         Vector3D close = expected.add(slope.scalarMultiply(100));
         Line line = new Line(expected.add(slope), close, 0);
         Transform xform = geoid.getBodyFrame().getTransformTo(frame, date);
         // transform to test frame
         close = xform.transformPosition(close);
         line = xform.transformLine(line);
 
         // action
         GeodeticPoint actual = geoid.getIntersectionPoint(line, close, frame,
                 date);
 
         // verify, 1 um position accuracy at Earth's surface
         assertThat(actual, geodeticPointCloseTo(gp, 1e-6));
     }
 
     /**
      * check {@link Geoid#getIntersectionPoint(Line, Vector3D, Frame,
      * AbsoluteDate)} returns null when there is no intersection
      */
     @Test
     public void testGetIntersectionPointNoIntersection() {
         Geoid geoid = getComponent();
         Vector3D closeMiss = new Vector3D(geoid.getEllipsoid()
                 .getEquatorialRadius() + 18, 0, 0);
         Line line = new Line(closeMiss, closeMiss.add(Vector3D.PLUS_J), 0);
 
         // action
         final GeodeticPoint actual = geoid.getIntersectionPoint(line,
                 closeMiss, geoid.getBodyFrame(), date);
 
         // verify
         assertThat(actual, nullValue());
     }
 
     /**
      * check altitude is referenced to the geoid. h<sub>ellipse</sub> =
      * h<sub>geoid</sub> + N. Where N is the undulation of the geoid.
      */
     @Test
     public void testTransformVector3DFrameAbsoluteDate()
             {
         // frame and date are the same
         Frame frame = FramesFactory.getGCRF();
         AbsoluteDate date = AbsoluteDate.CCSDS_EPOCH;
 
         Geoid geoid = getComponent();
         // test point at 0,0,0
         Vector3D point = new Vector3D(WGS84.getEquatorialRadius(), 0, 0);
         double undulation = geoid.getUndulation(0, 0, date);
 
         // check ellipsoidal points and geoidal points differ by undulation
         GeodeticPoint ellipsoidal = geoid.getEllipsoid().transform(
                 point, frame, date);
         GeodeticPoint geoidal = geoid.transform(point, frame, date);
         assertThat(ellipsoidal.getAltitude() - geoidal.getAltitude(),
                 is(undulation));
 
         // check it is the reverse of transform(GeodeticPoint)
         point = new Vector3D(0.5, 0.4, 0.31).scalarMultiply(WGS84
                 .getEquatorialRadius());
         Vector3D expected = geoid
                 .transform(geoid.transform(point, frame, date));
         // allow 2 upls of error
         assertThat(point, vectorCloseTo(expected, 2));
 
     }
 
     /**
      * check that the altitude is referenced to the geoid (includes
      * undulation).
      */
     @Test
     public void testTransformGeodeticPoint() {
         // geoid
         Geoid geoid = getComponent();
         // ellipsoid
         ReferenceEllipsoid ellipsoid = geoid.getEllipsoid();
         // point to test with orthometric height
         GeodeticPoint orthometric = new GeodeticPoint(0, 75, 5);
         // undulation at point
         double undulation = geoid.getUndulation(orthometric.getLatitude(),
                 orthometric.getLongitude(), date);
         // same point with height referenced to ellipsoid
         GeodeticPoint point = new GeodeticPoint(orthometric.getLatitude(),
                 orthometric.getLongitude(), orthometric.getAltitude()
                 + undulation);
 
         // test they are the same
         Vector3D expected = ellipsoid.transform(point);
         Vector3D actual = geoid.transform(orthometric);
         assertThat(actual, is(expected));
 
         // test the point 0,0,0
         expected = new Vector3D(WGS84.getEquatorialRadius()
                 + geoid.getUndulation(0, 0, date), 0, 0);
         actual = geoid.transform(new GeodeticPoint(0, 0, 0));
         assertThat(actual, vectorCloseTo(expected, 0));
     }
 
     /** check {@link Geoid#getEllipsoid()} */
     @Test
     public void testGetEllipsoid() {
         //setup
         Geoid geoid = new Geoid(potential, WGS84);
 
         //action + verify
         assertThat(geoid.getEllipsoid(), sameInstance(WGS84));
     }
 
     /**
      * check {@link Geoid#projectToGround(Vector3D, AbsoluteDate, Frame)}
      */
     @Test
     public void testProjectToGround() {
         //setup
         Vector3D p = new Vector3D(7e8, 1e3, 200);
         Geoid geoid = new Geoid(potential, WGS84);
 
         //action
         Vector3D actual = geoid.projectToGround(p, date, FramesFactory.getGCRF());
 
         //verify
         assertThat(
                 geoid.transform(actual, geoid.getBodyFrame(), date).getAltitude(),
                 closeTo(0.0, 1.1e-9));
     }
 
 }