Java tutorial
/* Copyright 2011-2012 Space Applications Services * Licensed to CS Communication & Systmes (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 java.io.BufferedReader; import java.io.File; import java.io.FileInputStream; import java.io.FileNotFoundException; import java.io.InputStream; import java.io.InputStreamReader; import java.util.Collection; import java.util.StringTokenizer; import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; import org.apache.commons.math3.util.FastMath; import org.junit.Assert; import org.junit.BeforeClass; import org.junit.Test; import org.orekit.Utils; import org.orekit.data.DataProvidersManager; 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.FramesFactory; import org.orekit.models.earth.GeoMagneticFieldFactory.FieldModel; import org.orekit.time.AbsoluteDate; import org.orekit.time.TimeScalesFactory; public class GeoMagneticFieldTest { /** 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; // test results for test values provided as part of the WMM2015 Report private final double[][] wmmTestValues = { // Date Alt Lat Lon X Y Z H F I D // km deg deg nT nT nT nT nT deg deg { 2015.0, 0, 80, 0, 6627.1, -445.9, 54432.3, 6642.1, 54836.0, 83.04, -3.85 }, { 2015.0, 0, 0, 120, 39518.2, 392.9, -11252.4, 39520.2, 41090.9, -15.89, 0.57 }, { 2015.0, 0, -80, 240, 5797.3, 15761.1, -52919.1, 16793.5, 55519.8, -72.39, 69.81 }, { 2015.0, 100, 80, 0, 6314.3, -471.6, 52269.8, 6331.9, 52652.0, 83.09, -4.27 }, { 2015.0, 100, 0, 120, 37535.6, 364.4, -10773.4, 37537.3, 39052.7, -16.01, 0.56 }, { 2015.0, 100, -80, 240, 5613.1, 14791.5, -50378.6, 15820.7, 52804.4, -72.57, 69.22 }, { 2017.5, 0, 80, 0, 6599.4, -317.1, 54459.2, 6607.0, 54858.5, 83.08, -2.75 }, { 2017.5, 0, 0, 120, 39571.4, 222.5, -11030.1, 39572.0, 41080.5, -15.57, 0.32 }, { 2017.5, 0, -80, 240, 5873.8, 15781.4, -52687.9, 16839.1, 55313.4, -72.28, 69.58 }, { 2017.5, 100, 80, 0, 6290.5, -348.5, 52292.7, 6300.1, 52670.9, 83.13, -3.17 }, { 2017.5, 100, 0, 120, 37585.5, 209.5, -10564.2, 37586.1, 39042.5, -15.70, 0.32 }, { 2017.5, 100, -80, 240, 5683.5, 14808.8, -50163.0, 15862.0, 52611.1, -72.45, 69.00 } }; // test results for test values provided as part of the WMM2015 Report // the results for the IGRF12 model have been obtained from the NOAA // online calculator: http://www.ngdc.noaa.gov/geomag-web/#igrfwmm private final double[][] igrfTestValues = { // Date Alt Lat Lon X Y Z H F I D // km deg deg nT nT nT nT nT deg deg { 2015.0, 0, 80, 0, 6630.9, -447.2, 54434.5, 6645.9, 54838.7, 83.039, -3.858 }, { 2015.0, 0, 0, 120, 39519.3, 388.6, -11251.7, 39521.3, 41091.7, -15.891, 0.563 }, { 2015.0, 0, -80, 240, 5808.8, 15754.8, -52945.5, 16791.5, 55544.4, -72.403, 69.761 }, { 2015.0, 100, 80, 0, 6317.2, -472.6, 52272.0, 6334.9, 52654.5, 83.090, -4.278 }, { 2015.0, 100, 0, 120, 37536.9, 361.2, -10773.1, 37538.6, 39053.9, -16.012, 0.551 }, { 2015.0, 100, -80, 240, 5622.8, 14786.8, -50401.4, 15819.8, 52825.8, -72.574, 69.180 }, { 2017.5, 0, 80, 0, 6601.0, -316.4, 54455.5, 6608.5, 54855.0, 83.080, -2.744 }, { 2017.5, 0, 0, 120, 39568.1, 225.0, -11041.4, 39568.7, 41080.3, -15.591, 0.325 }, { 2017.5, 0, -80, 240, 5894.7, 15768.1, -52696.8, 16833.9, 55320.2, -72.283, 69.502 }, { 2017.5, 100, 80, 0, 6291.6, -347.2, 52289.9, 6301.2, 52668.2, 83.128, -3.158 }, { 2017.5, 100, 0, 120, 37583.0, 212.3, -10575.1, 37583.6, 39043.0, -15.715, 0.323 }, { 2017.5, 100, -80, 240, 5702.0, 14797.8, -50170.0, 15858.3, 52616.7, -72.458, 68.927 } }; /** * load orekit data and gravity field. * * @throws Exception on error. */ @BeforeClass public static void setUpBefore() throws Exception { Utils.setDataRoot("earth:geoid:regular-data"); GravityFieldFactory.clearPotentialCoefficientsReaders(); GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("egm96", false)); potential = GravityFieldFactory.getConstantNormalizedProvider(maxDegree, maxOrder); } @Test public void testInterpolationYYY5() throws OrekitException { double decimalYear = GeoMagneticField.getDecimalYear(1, 1, 2005); GeoMagneticField field = GeoMagneticFieldFactory.getIGRF(decimalYear); GeoMagneticElements e = field.calculateField(1.2, 0.7, -2.5); Assert.assertEquals(-6.0032, e.getDeclination(), 1.0e-4); decimalYear = GeoMagneticField.getDecimalYear(2, 1, 2005); field = GeoMagneticFieldFactory.getIGRF(decimalYear); e = field.calculateField(1.2, 0.7, -2.5); Assert.assertEquals(-6.0029, e.getDeclination(), 1.0e-4); } @Test public void testInterpolationAtEndOfEpoch() throws OrekitException { double decimalYear = GeoMagneticField.getDecimalYear(31, 12, 2009); GeoMagneticField field1 = GeoMagneticFieldFactory.getIGRF(decimalYear); GeoMagneticField field2 = GeoMagneticFieldFactory.getIGRF(2010.0); Assert.assertNotEquals(field1.getEpoch(), field2.getEpoch()); GeoMagneticElements e1 = field1.calculateField(0, 0, 0); Assert.assertEquals(-6.1068, e1.getDeclination(), 1.0e-4); GeoMagneticElements e2 = field2.calculateField(0, 0, 0); Assert.assertEquals(-6.1064, e2.getDeclination(), 1.0e-4); } @Test public void testInterpolationAtEndOfValidity() throws OrekitException { double decimalYear = GeoMagneticField.getDecimalYear(1, 1, 2020); GeoMagneticField field = GeoMagneticFieldFactory.getIGRF(decimalYear); GeoMagneticElements e = field.calculateField(0, 0, 0); Assert.assertEquals(-4.7446, e.getDeclination(), 1.0e-4); } @Test(expected = OrekitException.class) public void testTransformationOutsideValidityPeriod() throws OrekitException { double decimalYear = GeoMagneticField.getDecimalYear(10, 1, 2020); @SuppressWarnings("unused") GeoMagneticField field = GeoMagneticFieldFactory.getIGRF(decimalYear); } @Test public void testWMM() throws Exception { // test values from sample coordinate file // provided as part of the geomag 7.0 distribution available at // http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html // modification: the julian day calculation of geomag is slightly different // to the one from the WMM code, we use the WMM convention thus the outputs // have been adapted. runSampleFile(FieldModel.WMM, "sample_coords.txt", "sample_out_WMM2015.txt"); final double eps = 1e-1; final double degreeEps = 1e-2; for (int i = 0; i < wmmTestValues.length; i++) { final GeoMagneticField model = GeoMagneticFieldFactory.getWMM(wmmTestValues[i][0]); final GeoMagneticElements result = model.calculateField(wmmTestValues[i][2], wmmTestValues[i][3], wmmTestValues[i][1]); // X Assert.assertEquals(wmmTestValues[i][4], result.getFieldVector().getX(), eps); // Y Assert.assertEquals(wmmTestValues[i][5], result.getFieldVector().getY(), eps); // Z Assert.assertEquals(wmmTestValues[i][6], result.getFieldVector().getZ(), eps); // H Assert.assertEquals(wmmTestValues[i][7], result.getHorizontalIntensity(), eps); // F Assert.assertEquals(wmmTestValues[i][8], result.getTotalIntensity(), eps); // inclination Assert.assertEquals(wmmTestValues[i][9], result.getInclination(), degreeEps); // declination Assert.assertEquals(wmmTestValues[i][10], result.getDeclination(), degreeEps); } } @Test public void testWMMWithHeightAboveMSL() throws Exception { // test results for test values provided as part of the WMM2015 Report // using height above MSL instead of height above ellipsoid // the results have been obtained from the NOAA online calculator: // http://www.ngdc.noaa.gov/geomag-web/#igrfwmm final double[][] testValues = { // Date Alt Lat Lon X Y Z H F I D // km deg deg nT nT nT nT nT deg deg { 2015.0, 100, 80, 0, 6314.2, -471.6, 52269.1, 6331.8, 52651.2, 83.093, -4.271 }, { 2015.0, 100, 0, 120, 37534.4, 364.3, -10773.1, 37536.2, 39051.6, -16.013, 0.556 }, { 2015.0, 100, -80, 240, 5613.2, 14791.9, -50379.6, 15821.1, 52805.4, -72.565, 69.219 } }; final Geoid geoid = new Geoid(potential, WGS84); final double eps = 1e-1; final double degreeEps = 1e-2; for (int i = 0; i < testValues.length; i++) { final AbsoluteDate date = new AbsoluteDate(2015, 1, 1, TimeScalesFactory.getUTC()); final GeoMagneticField model = GeoMagneticFieldFactory.getWMM(testValues[i][0]); final double undulation = geoid.getUndulation(FastMath.toRadians(testValues[i][2]), FastMath.toRadians(testValues[i][3]), date); final GeoMagneticElements result = model.calculateField(testValues[i][2], testValues[i][3], testValues[i][1] + undulation / 1000d); // X Assert.assertEquals(testValues[i][4], result.getFieldVector().getX(), eps); // Y Assert.assertEquals(testValues[i][5], result.getFieldVector().getY(), eps); // Z Assert.assertEquals(testValues[i][6], result.getFieldVector().getZ(), eps); // H Assert.assertEquals(testValues[i][7], result.getHorizontalIntensity(), eps); // F Assert.assertEquals(testValues[i][8], result.getTotalIntensity(), eps); // inclination Assert.assertEquals(testValues[i][9], result.getInclination(), degreeEps); // declination Assert.assertEquals(testValues[i][10], result.getDeclination(), degreeEps); } } @Test public void testIGRF() throws Exception { // test values from sample coordinate file // provided as part of the geomag 7.0 distribution available at // http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html // modification: the julian day calculation of geomag is slightly different // to the one from the WMM code, we use the WMM convention thus the outputs // have been adapted. runSampleFile(FieldModel.IGRF, "sample_coords.txt", "sample_out_IGRF12.txt"); final double eps = 1e-1; final double degreeEps = 1e-2; for (int i = 0; i < igrfTestValues.length; i++) { final GeoMagneticField model = GeoMagneticFieldFactory.getIGRF(igrfTestValues[i][0]); final GeoMagneticElements result = model.calculateField(igrfTestValues[i][2], igrfTestValues[i][3], igrfTestValues[i][1]); final Vector3D b = result.getFieldVector(); // X Assert.assertEquals(igrfTestValues[i][4], b.getX(), eps); // Y Assert.assertEquals(igrfTestValues[i][5], b.getY(), eps); // Z Assert.assertEquals(igrfTestValues[i][6], b.getZ(), eps); // H Assert.assertEquals(igrfTestValues[i][7], result.getHorizontalIntensity(), eps); // F Assert.assertEquals(igrfTestValues[i][8], result.getTotalIntensity(), eps); // inclination Assert.assertEquals(igrfTestValues[i][9], result.getInclination(), degreeEps); // declination Assert.assertEquals(igrfTestValues[i][10], result.getDeclination(), degreeEps); } } @Test(expected = OrekitException.class) public void testUnsupportedTransform() throws Exception { final GeoMagneticField model = GeoMagneticFieldFactory.getIGRF(1910); // the IGRF model of 1910 does not have secular variation, thus time transformation is not supported model.transformModel(1950); } @Test(expected = OrekitException.class) public void testOutsideValidityTransform() throws Exception { final GeoMagneticField model1 = GeoMagneticFieldFactory.getIGRF(2005); final GeoMagneticField model2 = GeoMagneticFieldFactory.getIGRF(2010); // the interpolation transformation is only allowed between 2005 and 2010 model1.transformModel(model2, 2012); } @Test public void testValidTransform() throws Exception { final GeoMagneticField model = GeoMagneticFieldFactory.getWMM(2015); Assert.assertTrue(model.supportsTimeTransform()); final GeoMagneticField transformedModel = model.transformModel(2017); Assert.assertEquals(2015, transformedModel.validFrom(), 1e0); Assert.assertEquals(2020, transformedModel.validTo(), 1e0); Assert.assertEquals(2017, transformedModel.getEpoch(), 1e0); } @Test public void testLoadOriginalWMMModel() throws Exception { GeoMagneticModelLoader loader = new GeoMagneticModelLoader(); InputStream input = getResource("WMM2015.COF"); loader.loadData(input, "WMM2015.COF"); Collection<GeoMagneticField> models = loader.getModels(); Assert.assertNotNull(models); Assert.assertEquals(1, models.size()); GeoMagneticField wmmModel = models.iterator().next(); Assert.assertEquals("WMM-2015", wmmModel.getModelName()); Assert.assertEquals(2015, wmmModel.getEpoch(), 1e-9); final double eps = 1e-1; final double degreeEps = 1e-2; for (int i = 0; i < wmmTestValues.length; i++) { if (wmmTestValues[i][0] != wmmModel.getEpoch()) { continue; } final GeoMagneticElements result = wmmModel.calculateField(wmmTestValues[i][2], wmmTestValues[i][3], wmmTestValues[i][1]); // X Assert.assertEquals(wmmTestValues[i][4], result.getFieldVector().getX(), eps); // Y Assert.assertEquals(wmmTestValues[i][5], result.getFieldVector().getY(), eps); // Z Assert.assertEquals(wmmTestValues[i][6], result.getFieldVector().getZ(), eps); // H Assert.assertEquals(wmmTestValues[i][7], result.getHorizontalIntensity(), eps); // F Assert.assertEquals(wmmTestValues[i][8], result.getTotalIntensity(), eps); // inclination Assert.assertEquals(wmmTestValues[i][9], result.getInclination(), degreeEps); // declination Assert.assertEquals(wmmTestValues[i][10], result.getDeclination(), degreeEps); } } public void runSampleFile(final FieldModel type, final String inputFile, final String outputFile) throws Exception { final BufferedReader inReader = new BufferedReader(new InputStreamReader(getResource(inputFile))); final BufferedReader outReader = new BufferedReader(new InputStreamReader(getResource(outputFile))); // read header line outReader.readLine(); String line = null; while ((line = inReader.readLine()) != null) { if (line.trim().length() == 0) { break; } final StringTokenizer st = new StringTokenizer(line); final double year = getYear(st.nextToken()); final String heightType = st.nextToken(); final String heightStr = st.nextToken(); final double lat = getLatLon(st.nextToken()); final double lon = getLatLon(st.nextToken()); final GeoMagneticField field = GeoMagneticFieldFactory.getField(type, year); double height = Double.valueOf(heightStr.substring(1)); if (heightStr.startsWith("M")) { // convert from m to km height /= 1000d; } else if (heightStr.startsWith("F")) { // convert from feet to km height *= 3.048e-4; } final GeoMagneticElements ge = field.calculateField(lat, lon, height); final String validateLine = outReader.readLine(); // geocentric altitude is not yet supported, ignore by now if (!heightType.startsWith("C")) { validate(ge, validateLine); } String geString = ge.toString(); Assert.assertNotNull(geString); Assert.assertFalse(geString.isEmpty()); } } private double getYear(final String yearStr) { if (yearStr.contains(",")) { final StringTokenizer st = new StringTokenizer(yearStr, ","); final String y = st.nextToken(); final String m = st.nextToken(); final String d = st.nextToken(); return GeoMagneticField.getDecimalYear(Integer.valueOf(d), Integer.valueOf(m), Integer.valueOf(y)); } else { return Double.valueOf(yearStr); } } private double getLatLon(final String str) { if (str.contains(",")) { final StringTokenizer st = new StringTokenizer(str, ","); final int d = Integer.valueOf(st.nextToken()); final int m = Integer.valueOf(st.nextToken()); int s = 0; if (st.hasMoreTokens()) { s = Integer.valueOf(st.nextToken()); } double deg = Math.abs(d) + m / 60d + s / 3600d; if (d < 0) { deg = -deg; } return deg; } else { return Double.valueOf(str); } } private double getDegree(final String degree, final String minute) { double result = Double.valueOf(degree.substring(0, degree.length() - 1)); final double min = Double.valueOf(minute.substring(0, minute.length() - 1)) / 60d; result += (result < 0) ? -min : min; return result; } @SuppressWarnings("unused") private void validate(final GeoMagneticElements ge, final String outputLine) throws Exception { final StringTokenizer st = new StringTokenizer(outputLine); final double year = getYear(st.nextToken()); final String coord = st.nextToken(); final String heightStr = st.nextToken(); final String latStr = st.nextToken(); final String lonStr = st.nextToken(); final double dec = getDegree(st.nextToken(), st.nextToken()); final double inc = getDegree(st.nextToken(), st.nextToken()); final double h = Double.valueOf(st.nextToken()); final double x = Double.valueOf(st.nextToken()); final double y = Double.valueOf(st.nextToken()); final double z = Double.valueOf(st.nextToken()); final double f = Double.valueOf(st.nextToken()); final double eps = 1e-1; Assert.assertEquals(h, ge.getHorizontalIntensity(), eps); Assert.assertEquals(f, ge.getTotalIntensity(), eps); Assert.assertEquals(x, ge.getFieldVector().getX(), eps); Assert.assertEquals(y, ge.getFieldVector().getY(), eps); Assert.assertEquals(z, ge.getFieldVector().getZ(), eps); Assert.assertEquals(dec, ge.getDeclination(), eps); Assert.assertEquals(inc, ge.getInclination(), eps); } private InputStream getResource(final String name) throws FileNotFoundException { // the data path has multiple components, the resources are in the first one final String separator = System.getProperty("path.separator"); final String dataPath = System.getProperty(DataProvidersManager.OREKIT_DATA_PATH).split(separator)[0]; return new FileInputStream(new File(dataPath, name)); } }