Java tutorial
/* Copyright 2002-2015 CS Systmes d'Information * 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.forces.drag; import java.io.FileNotFoundException; import java.text.ParseException; 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.Test; import org.orekit.SolarInputs97to05; import org.orekit.Utils; import org.orekit.bodies.CelestialBodyFactory; import org.orekit.bodies.GeodeticPoint; import org.orekit.bodies.OneAxisEllipsoid; import org.orekit.errors.OrekitException; import org.orekit.frames.Frame; import org.orekit.frames.FramesFactory; import org.orekit.time.AbsoluteDate; import org.orekit.time.DateComponents; import org.orekit.time.TimeComponents; import org.orekit.time.TimeScalesFactory; import org.orekit.utils.Constants; import org.orekit.utils.IERSConventions; import org.orekit.utils.PVCoordinatesProvider; public class JB2006Test { @Test public void testWithOriginalTestsCases() throws OrekitException, ParseException { Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true); PVCoordinatesProvider sun = CelestialBodyFactory.getSun(); OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.460, 1.0 / 298.257222101, itrf); SolarInputs97to05 in = SolarInputs97to05.getInstance(); earth.setAngularThreshold(1e-10); JB2006 atm = new JB2006(in, sun, earth); double myRo; double PI = 3.1415927; // SET SOLAR INDICES USE 1 DAY LAG FOR EUV AND F10 INFLUENCE double S10 = 140; double S10B = 100; double F10 = 135; double F10B = 95; // USE 5 DAY LAG FOR MG FUV INFLUENCE double XM10 = 130; double XM10B = 95; // USE 6.7 HR LAG FOR ap INFLUENCE double AP = 30; // SET TIME OF INTEREST double IYR = 01; double IDAY = 200; if (IYR < 50) IYR = IYR + 100; double IYY = (IYR - 50) * 365 + ((IYR - 1) / 4 - 12); double ID1950 = IYY + IDAY; double D1950 = ID1950; double AMJD = D1950 + 33281; // COMPUTE DENSITY KG/M3 RHO // alt = 400 myRo = atm.getDensity(AMJD, 90. * PI / 180., 20. * PI / 180., 90. * PI / 180., 45. * PI / 180., 1000 * 400, F10, F10B, AP, S10, S10B, XM10, XM10B); double roTestCase = 0.4066e-11; double tzTestCase = 1137.7; double tinfTestCase = 1145.8; Assert.assertEquals(roTestCase * 1e12, FastMath.round(myRo * 1e15) / 1e3, 0); Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp() * 10) / 10.0, 0); Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp() * 10) / 10.0, 0); // alt = 90 myRo = atm.getDensity(AMJD, 90. * PI / 180., 20. * PI / 180., 90. * PI / 180., 45. * PI / 180., 1000 * 90, F10, F10B, AP, S10, S10B, XM10, XM10B); roTestCase = 0.3285e-05; tzTestCase = 183.0; tinfTestCase = 1142.8; Assert.assertEquals(roTestCase * 1e05, FastMath.round(myRo * 1e09) / 1e4, 0); Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp() * 10) / 10.0, 0); Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp() * 10) / 10.0, 0); // alt = 110 myRo = atm.getDensity(AMJD, 90. * PI / 180., 20. * PI / 180., 90. * PI / 180., 45. * PI / 180., 1000 * 110, F10, F10B, AP, S10, S10B, XM10, XM10B); roTestCase = 0.7587e-07; tzTestCase = 257.4; tinfTestCase = 1142.8; Assert.assertEquals(roTestCase * 1e07, FastMath.round(myRo * 1e11) / 1e4, 0); Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp() * 10) / 10.0, 0); Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp() * 10) / 10.0, 0); // alt = 180 myRo = atm.getDensity(AMJD, 90. * PI / 180., 20. * PI / 180., 90. * PI / 180., 45. * PI / 180., 1000 * 180, F10, F10B, AP, S10, S10B, XM10, XM10B); roTestCase = 0.5439; // *1e-9 tzTestCase = 915.0; tinfTestCase = 1130.9; Assert.assertEquals(roTestCase, FastMath.round(myRo * 1e13) / 1e4, 0); Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp() * 10) / 10.0, 0); Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp() * 10) / 10.0, 0); // alt = 230 myRo = atm.getDensity(AMJD, 90. * PI / 180., 20. * PI / 180., 90. * PI / 180., 45. * PI / 180., 1000 * 230, F10, F10B, AP, S10, S10B, XM10, XM10B); roTestCase = 0.1250e-09; tzTestCase = 1047.5; tinfTestCase = 1137.4; Assert.assertEquals(roTestCase * 1e09, FastMath.round(myRo * 1e13) / 1e4, 0); Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp() * 10) / 10.0, 0); Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp() * 10) / 10.0, 0); // alt = 270 myRo = atm.getDensity(AMJD, 90. * PI / 180., 20. * PI / 180., 90. * PI / 180., 45. * PI / 180., 1000 * 270, F10, F10B, AP, S10, S10B, XM10, XM10B); roTestCase = 0.4818e-10; tzTestCase = 1095.6; tinfTestCase = 1142.5; Assert.assertEquals(roTestCase * 1e10, FastMath.round(myRo * 1e14) / 1e4, 0); Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp() * 10) / 10.0, 0); Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp() * 10) / 10.0, 0); // alt = 660 myRo = atm.getDensity(AMJD, 90. * PI / 180., 20. * PI / 180., 90. * PI / 180., 45. * PI / 180., 1000 * 660, F10, F10B, AP, S10, S10B, XM10, XM10B); roTestCase = 0.9451e-13; tzTestCase = 1149.0; tinfTestCase = 1149.9; Assert.assertEquals(roTestCase * 1e13, FastMath.round(myRo * 1e17) / 1e4, 0); Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp() * 10) / 10.0, 0); Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp() * 10) / 10.0, 0); // alt = 890 myRo = atm.getDensity(AMJD, 90. * PI / 180., 20. * PI / 180., 90. * PI / 180., 45. * PI / 180., 1000 * 890, F10, F10B, AP, S10, S10B, XM10, XM10B); roTestCase = 0.8305e-14; tzTestCase = 1142.5; tinfTestCase = 1142.8; Assert.assertEquals(roTestCase * 1e14, FastMath.round(myRo * 1e18) / 1e4, 0); Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp() * 10) / 10.0, 0); Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp() * 10) / 10.0, 0); // alt = 1320 myRo = atm.getDensity(AMJD, 90. * PI / 180., 20. * PI / 180., 90. * PI / 180., 45. * PI / 180., 1000 * 1320, F10, F10B, AP, S10, S10B, XM10, XM10B); roTestCase = 0.2004e-14; tzTestCase = 1142.7; tinfTestCase = 1142.8; Assert.assertEquals(roTestCase * 1e14, FastMath.round(myRo * 1e18) / 1e4, 0); Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp() * 10) / 10.0, 0); Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp() * 10) / 10.0, 0); // alt = 1600 myRo = atm.getDensity(AMJD, 90. * PI / 180., 20. * PI / 180., 90. * PI / 180., 45. * PI / 180., 1000 * 1600, F10, F10B, AP, S10, S10B, XM10, XM10B); roTestCase = 0.1159e-14; tzTestCase = 1142.8; tinfTestCase = 1142.8; Assert.assertEquals(roTestCase * 1e14, FastMath.round(myRo * 1e18) / 1e4, 0); Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp() * 10) / 10.0, 0); Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp() * 10) / 10.0, 0); // OTHER entries AMJD += 50; myRo = atm.getDensity(AMJD, 45. * PI / 180., 10. * PI / 180., 45. * PI / 180., -10. * PI / 180., 400 * 1000, F10, F10B, AP, S10, S10B, XM10, XM10B); roTestCase = 0.4838e-11; tzTestCase = 1137.4; tinfTestCase = 1145.4; Assert.assertEquals(roTestCase * 1e11, FastMath.round(myRo * 1e15) / 1e4, 0); Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp() * 10) / 10.0, 0); Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp() * 10) / 10.0, 0); } public void testComparisonWithDTM2000() throws OrekitException, ParseException, FileNotFoundException { AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 01, 01), TimeComponents.H00, TimeScalesFactory.getUTC()); Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true); PVCoordinatesProvider sun = CelestialBodyFactory.getSun(); OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.460, 1.0 / 298.257222101, itrf); SolarInputs97to05 in = SolarInputs97to05.getInstance(); earth.setAngularThreshold(1e-10); JB2006 jb = new JB2006(in, sun, earth); DTM2000 dtm = new DTM2000(in, sun, earth); // Positions Vector3D pos = new Vector3D(6500000.0, -1234567.0, 4000000.0); // COMPUTE DENSITY KG/M3 RHO // alt = 400 double roJb = jb.getDensity(date, pos, FramesFactory.getEME2000()); double roDtm = dtm.getDensity(date, pos, FramesFactory.getEME2000()); pos = new Vector3D(3011109.360780633, -5889822.669411588, 4002170.0385907636); // COMPUTE DENSITY KG/M3 RHO // alt = 400 roJb = jb.getDensity(date, pos, FramesFactory.getEME2000()); roDtm = dtm.getDensity(date, pos, FramesFactory.getEME2000()); pos = new Vector3D(-1033.4793830 * 1000, 7901.2952754 * 1000, 6380.3565958 * 1000); // COMPUTE DENSITY KG/M3 RHO // alt = 400 roJb = jb.getDensity(date, pos, FramesFactory.getEME2000()); roDtm = dtm.getDensity(date, pos, FramesFactory.getEME2000()); GeodeticPoint point; for (int i = 0; i < 367; i++) { date = date.shiftedBy(Constants.JULIAN_DAY); point = new GeodeticPoint(FastMath.toRadians(40), 0, 300 * 1000); pos = earth.transform(point); roJb = jb.getDensity(date, pos, FramesFactory.getEME2000()); roDtm = dtm.getDensity(date, pos, FramesFactory.getEME2000()); Assert.assertEquals(roDtm, roJb, roJb); } } public void testSolarInputs() throws OrekitException, ParseException { AbsoluteDate date = new AbsoluteDate(new DateComponents(2001, 01, 14), TimeComponents.H00, TimeScalesFactory.getUTC()); SolarInputs97to05 in = SolarInputs97to05.getInstance(); // 2001 14 2451924.0 176.3 164.4 180.0 180.4 163.4 169.2 // 14 176 164 9 12 9 6 4 4 9 7 Assert.assertEquals(176.3, in.getF10(date), 0); Assert.assertEquals(164.4, in.getF10B(date), 0); Assert.assertEquals(180.0, in.getS10(date), 0); Assert.assertEquals(180.4, in.getS10B(date), 0); Assert.assertEquals(163.4, in.getXM10(date), 0); Assert.assertEquals(169.2, in.getXM10B(date), 0); Assert.assertEquals(9, in.getAp(date), 0); date = date.shiftedBy(11 * 3600); Assert.assertEquals(6, in.getAp(date), 0); date = new AbsoluteDate(new DateComponents(1998, 02, 02), new TimeComponents(18, 00, 00), TimeScalesFactory.getUTC()); // 1998 33 2450847.0 89.1 95.1 95.8 97.9 81.3 92.0 1 // 33 89 95 4 5 4 4 2 0 0 3 98 Assert.assertEquals(89.1, in.getF10(date), 0); Assert.assertEquals(95.1, in.getF10B(date), 0); Assert.assertEquals(95.8, in.getS10(date), 0); Assert.assertEquals(97.9, in.getS10B(date), 0); Assert.assertEquals(81.3, in.getXM10(date), 0); Assert.assertEquals(92.0, in.getXM10B(date), 0); Assert.assertEquals(0, in.getAp(date), 0); date = date.shiftedBy(6 * 3600 - 1); Assert.assertEquals(3, in.getAp(date), 0); } @Before public void setUp() { Utils.setDataRoot("regular-data"); } }