Java tutorial
/* Contributed in the public domain. * Licensed to CS Systmes d'Information (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 org.apache.commons.math3.geometry.euclidean.threed.Line; import org.apache.commons.math3.geometry.euclidean.threed.Rotation; import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; import org.apache.commons.math3.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.errors.OrekitException; 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; import java.util.Arrays; 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; /** * Unit tests for {@link Geoid}. * * @author Evan Ward */ public class GeoidTest { /** maximum degree and order used in testing {@link Geoid}. */ @SuppressWarnings("javadoc") private static final int maxOrder = 360, 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 = potential.getReferenceDate(); } /** {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. * * @throws OrekitException on error * @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() throws OrekitException { /* * 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. * * @throws OrekitException on error */ @Test public void testGetIntersectionPoint() throws OrekitException { // 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 um. assertThat(message, actualReversed, geodeticPointCloseTo(gp, 1e-6)); assertThat(message, actual, geodeticPointCloseTo(gp, 1e-6)); } } /** * check {@link Geoid#getIntersectionPoint(Line, Vector3D, Frame, * AbsoluteDate)} handles frame transformations correctly * * @throws OrekitException on error */ @Test public void testGetIntersectionPointFrame() throws OrekitException { // 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 * * @throws OrekitException on error */ @Test public void testGetIntersectionPointNoIntersection() throws OrekitException { 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. * * @throws OrekitException on error */ @Test public void testTransformVector3DFrameAbsoluteDate() throws OrekitException { // 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). * * @throws OrekitException on error */ @Test public void testTransformGeodeticPoint() throws OrekitException { // 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)} * * @throws OrekitException on error */ @Test public void testProjectToGround() throws OrekitException { //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, 1e-9)); } }