diff --git a/python_files/test/AEMTest.py b/python_files/test/AEMTest.py deleted file mode 100644 index 9e73c26b241798088db1df37d9abdb04b9ae888d..0000000000000000000000000000000000000000 --- a/python_files/test/AEMTest.py +++ /dev/null @@ -1,216 +0,0 @@ -# -*- coding: utf-8 -*- - -""" - -/* Copyright 2002-2019 CS Syst��mes d'Information - * Licensed to CS Syst��mes 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. - */ - -Python version translated from Java by Olivier Podevin, CS Group 2020 - - """ - -import orekit - -orekit.initVM() -from orekit.pyhelpers import setup_orekit_curdir, datetime_to_absolutedate - -setup_orekit_curdir(filename="resources.zip") - -import unittest -import sys -import os -import math -from datetime import datetime - -from java.util import ArrayList, List, HashMap, Map -from java.lang import Double, StringBuffer -from java.io import BufferedReader, StringReader - -from org.hipparchus.geometry.euclidean.threed import Rotation, Vector3D - -from org.orekit.attitudes import AttitudeProvider, InertialProvider -from org.orekit.bodies import CelestialBodyFactory -from org.orekit.files.ccsds import AEMParser, AEMWriter, AEMFile, Keyword, StreamingAemWriter -from org.orekit.frames import FramesFactory -from org.orekit.orbits import CartesianOrbit -from org.orekit.propagation.analytical import KeplerianPropagator -from org.orekit.time import AbsoluteDate, TimeScalesFactory, TimeScale -from org.orekit.utils import IERSConventions, PVCoordinates - -class AEMTest(unittest.TestCase): - - def testAemParsing(self): - """Test the parsing of an AEM file functionality""" - inFilename = "./resources/ccsds/AEMExample.txt" - aemParser = AEMParser() - aemParser = aemParser.withMu(CelestialBodyFactory.getEarth().getGM()). \ - withConventions(IERSConventions.IERS_2010). \ - withSimpleEOP(True) - aemFile = aemParser.parse(inFilename) - - ####### Checks - - # Start date - startRef = datetime_to_absolutedate(datetime(1996, 11, 28, 22, 8, 2, 555500)) - startFile = aemFile.getSatellites().get("1996-062A").getStart() - # print("start ref=" + startRef.toString() + ", type=" + str(type(startRef))) - # print("start file=" + startFile.toString() + ", type=" + str(type(startFile))) - # print("delta=" + str(startRef.durationFrom(startFile))) - self.assertTrue(math.isclose(startRef.durationFrom(startFile), 0, rel_tol=0.0, abs_tol=Double.MIN_VALUE)) - - # Stop date - stopRef = datetime_to_absolutedate(datetime(1996, 12, 28, 21, 23, 0, 555500)) - stopFile = aemFile.getSatellites().get("1996-062A").getStop() - self.assertTrue(math.isclose(stopRef.durationFrom(stopFile), 0, rel_tol=0.0, abs_tol=Double.MIN_VALUE)) - - def testAemWriting(self): - """Test the writing of an AEM file functionality""" - inFilename = "./resources/ccsds/AEMExample.txt" - - aemParser = AEMParser() - aemParser = aemParser.withMu(CelestialBodyFactory.getEarth().getGM()) \ - .withConventions(IERSConventions.IERS_2010) - - aemFile = aemParser.parse(inFilename) - - ####### Read values - originator = aemFile.getOriginator() - objectName = aemFile.getAttitudeBlocks().get(0).getMetaData().getObjectName() - objectID = aemFile.getAttitudeBlocks().get(0).getMetaData().getObjectID() - - ####### Create AEM file - outputDir = "./output/ccsds" - tempFilename = "TestWriteAEM1.aem" - tempAEMFilePath = os.path.join(outputDir, tempFilename) - - # Create output dir if needed - if (not os.path.isdir(outputDir)): - os.makedirs(outputDir) - - aemWriter = AEMWriter(originator, objectID, objectName) - aemWriter.write(tempAEMFilePath, aemFile) - - ####### Checks - generatedAemFile = aemParser.parse(tempAEMFilePath) - originatorOut = aemFile.getOriginator() - objectNameOut = aemFile.getAttitudeBlocks().get(0).getMetaData().getObjectName() - objectIDOut = aemFile.getAttitudeBlocks().get(0).getMetaData().getObjectID() - - self.assertEqual(originator, originatorOut) - self.assertEqual(objectName, objectNameOut) - self.assertEqual(objectID, objectIDOut) - - def testAemStreamWriting(self): - """Test the stream writing of an AEM file functionality""" - - QUATERNION_PRECISION = 1e-5 - DATE_PRECISION = 1e-3 - - # Time scale - utc = TimeScalesFactory.getUTC() - - inFilename = "./resources/ccsds/AEMExample.txt" - - # Create a list of files - files = ["./resources/ccsds/AEMExample7.txt"] - - for ex in files: - # Reference AEM file - aemParser = AEMParser() - aemParser = aemParser.withMu(CelestialBodyFactory.getEarth().getGM()) \ - .withConventions(IERSConventions.IERS_2010) - - aemFile = aemParser.parse(ex) - - # Satellite attitude ephemeris as read from the reference file - satellite = aemFile.getSatellites().values().iterator().next() - - # Meta data are extracted from the reference file - originator = aemFile.getOriginator() - objectName = satellite.getSegments().get(0).getMetaData().getObjectName() - objectID = satellite.getSegments().get(0).getMetaData().getObjectID() - headerCmt = aemFile.getHeaderComment().get(0) - attitudeDir = satellite.getSegments().get(0).getAttitudeDirection() - refFrameA = satellite.getSegments().get(0).getRefFrameAString() - refFrameB = satellite.getSegments().get(0).getRefFrameBString() - attitudeType = satellite.getSegments().get(0).getAttitudeType() - isFirst = "LAST" - - metadata = HashMap() - metadata.put(Keyword.ORIGINATOR, originator) - metadata.put(Keyword.OBJECT_NAME, "will be overwritten") - metadata.put(Keyword.OBJECT_ID, objectID) - metadata.put(Keyword.COMMENT, headerCmt) - - segmentData = HashMap() - segmentData.put(Keyword.OBJECT_NAME, objectName) - segmentData.put(Keyword.ATTITUDE_DIR, attitudeDir) - segmentData.put(Keyword.QUATERNION_TYPE, isFirst) - segmentData.put(Keyword.ATTITUDE_TYPE, attitudeType) - segmentData.put(Keyword.REF_FRAME_A, refFrameA) - segmentData.put(Keyword.REF_FRAME_B, refFrameB.replace(' ', '_')) - - # Initialize a Keplerian propagator with an Inertial attitude provider - # It is expected that all attitude data lines will have the same value - buffer = StringBuffer() - streamingAemWriter = StreamingAemWriter(buffer, utc, metadata) - streamingAemWriter.writeHeader() - segment = streamingAemWriter.newSegment(segmentData) - - #### Create a Keplerian propagator. - startingDate = satellite.getSegments().get(0).getStart() - attitudeProvider = InertialProvider(satellite.getSegments().get(0).getAngularCoordinates().get(0).getRotation(), \ - FramesFactory.getEME2000()) - position = Vector3D(-29536113.0, 30329259.0, -100125.0) - velocity = Vector3D(-2194.0, -2141.0, -8.0) - pvCoordinates = PVCoordinates( position, velocity) - mu = 3.9860047e14 - cartesianOrbit = CartesianOrbit(pvCoordinates, FramesFactory.getEME2000(), startingDate, mu) - propagator = KeplerianPropagator(cartesianOrbit, attitudeProvider) - - # We propagate 60 seconds after the start date with a step equals to 10.0 seconds - # It is expected to have an attitude data block containing 7 data lines - step = 10.0 - propagator.setMasterMode(step, segment) - propagator.propagate(satellite.getSegments().get(0).getStart().shiftedBy(60.0)) - - # Generated AEM file - reader = BufferedReader(StringReader(buffer.toString())) - generatedAemFile = aemParser.parse(reader, "buffer") - - # There is only one attitude ephemeris block - self.assertEqual(1, generatedAemFile.getAttitudeBlocks().size()) - # There are 7 data lines in the attitude ephemeris block - ac = generatedAemFile.getAttitudeBlocks().get(0).getAngularCoordinates() - self.assertEqual(7, ac.size()) - - # Verify - i = 0 - while (i < 7): - #self.assertEqual(step * i, ac.get(i).getDate().durationFrom(satellite.getSegments().get(0).getStart()), DATE_PRECISION) - self.assertTrue(math.isclose(step * i, ac.get(i).getDate().durationFrom(satellite.getSegments().get(0).getStart()), rel_tol=0.0, abs_tol=DATE_PRECISION)) - rot = ac.get(i).getRotation() - self.assertTrue(math.isclose(0.68427, rot.getQ0(), rel_tol=0.0, abs_tol=QUATERNION_PRECISION)) - self.assertTrue(math.isclose(0.56748, rot.getQ1(), rel_tol=0.0, abs_tol=QUATERNION_PRECISION)) - self.assertTrue(math.isclose(0.03146, rot.getQ2(), rel_tol=0.0, abs_tol=QUATERNION_PRECISION)) - self.assertTrue(math.isclose(0.45689, rot.getQ3(), rel_tol=0.0, abs_tol=QUATERNION_PRECISION)) - i += 1 - -if __name__ == '__main__': - suite = unittest.TestLoader().loadTestsFromTestCase(AEMTest) - ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() - sys.exit(ret) diff --git a/python_files/test/AbstractDetectorTest.py b/python_files/test/AbstractDetectorTest.py index 5e0621a94fc62895991a9069c0660751c5abc5b7..64c60f15414bf1692c1b19d1c0c83aefb2097b9d 100644 --- a/python_files/test/AbstractDetectorTest.py +++ b/python_files/test/AbstractDetectorTest.py @@ -23,16 +23,16 @@ import unittest import sys from orekit.pyhelpers import setup_orekit_curdir -setup_orekit_curdir() +setup_orekit_curdir("resources") class PassCounter(PythonEventHandler): """Eventhandler that counts positive events""" passes = 0 - def init(self, initialstate, target): + def init(self, initialstate, target, detector): pass - def eventOccurred(self, s, T, increasing): + def eventOccurred(self, s, detector, increasing): if increasing: self.passes = self.passes + 1 diff --git a/python_files/test/AdditionalEquationsTest.py b/python_files/test/AdditionalEquationsTest.py index 12f4437f8300a711b5a036969b0d89816798bbb0..6056afef80713e35ed7c01682cf79cfefeb808cd 100644 --- a/python_files/test/AdditionalEquationsTest.py +++ b/python_files/test/AdditionalEquationsTest.py @@ -31,7 +31,7 @@ import orekit orekit.initVM() from orekit import JArray_double -from org.orekit.data import DataProvidersManager, ZipJarCrawler +from org.orekit.data import DataProvidersManager, ZipJarCrawler, DataContext, DirectoryCrawler from org.orekit.propagation.integration import PythonAdditionalEquations from org.orekit.propagation.integration import AdditionalEquations from org.orekit.forces.gravity.potential import GravityFieldFactory @@ -83,12 +83,12 @@ class InitCheckerEquations(PythonAdditionalEquations): # implements AdditionalE class AdditionalEquationsTest(unittest.TestCase): def setUp(self): - DM = DataProvidersManager.getInstance() - datafile = File('resources.zip') + DM = DataContext.getDefault().getDataProvidersManager() + datafile = File('resources') if not datafile.exists(): print('File :', datafile.absolutePath, ' not found') - crawler = ZipJarCrawler(datafile) + crawler = DirectoryCrawler(datafile) DM.clearProviders() DM.addProvider(crawler) System.setProperty(DataProvidersManager.OREKIT_DATA_PATH, 'potential/shm-format') diff --git a/python_files/test/AltitudeDetectorTest.py b/python_files/test/AltitudeDetectorTest.py index 8308159bfdc0faa4544321bd90f3fb25bc1fdd45..623bcf55818602be254f668888d5133e52f72c34 100644 --- a/python_files/test/AltitudeDetectorTest.py +++ b/python_files/test/AltitudeDetectorTest.py @@ -27,7 +27,7 @@ Python version translated from Java by Petrus Hyvönen, SSC 2014 import orekit orekit.initVM() from orekit.pyhelpers import setup_orekit_curdir -setup_orekit_curdir() # orekit-data.zip shall be in current dir +setup_orekit_curdir("resources") import unittest import sys diff --git a/python_files/test/BackAndForthDetectorTest.py b/python_files/test/BackAndForthDetectorTest.py index afe000bee2ea31a724eb3b0fb1eb78b2261d7663..7cfb3dfb095fe4105e954ea991db4492704b35af 100644 --- a/python_files/test/BackAndForthDetectorTest.py +++ b/python_files/test/BackAndForthDetectorTest.py @@ -26,7 +26,7 @@ Python version translated from Java by Petrus Hyvönen, SSC 2014 import orekit orekit.initVM() from orekit.pyhelpers import setup_orekit_curdir -setup_orekit_curdir() +setup_orekit_curdir("resources") from org.orekit.bodies import GeodeticPoint from org.orekit.bodies import OneAxisEllipsoid @@ -54,7 +54,7 @@ class Visibility(PythonEventHandler): # implements EventHandler<ElevationDetecto self._visiNb = 0 super(Visibility, self).__init__() - def init(self, initialstate, target): + def init(self, initialstate, target, detector): pass def getVisiNb(self): diff --git a/python_files/test/BrouwerLyddanePropagatorTest.py b/python_files/test/BrouwerLyddanePropagatorTest.py new file mode 100644 index 0000000000000000000000000000000000000000..7d455fb8a4d16211c3ee14067aa6706578dacc5a --- /dev/null +++ b/python_files/test/BrouwerLyddanePropagatorTest.py @@ -0,0 +1,147 @@ +# -*- coding: utf-8 -*- + +""" + +Python version translated from Java orekit 11.1 by Petrus Hyvönen, SSC 2022 + + """ + +import orekit +orekit.initVM() +from orekit.pyhelpers import setup_orekit_curdir + +from org.hipparchus.geometry.euclidean.threed import Vector3D +from org.hipparchus.ode.nonstiff import DormandPrince853Integrator +from org.orekit.bodies import GeodeticPoint +from org.orekit.bodies import OneAxisEllipsoid +from org.orekit.frames import TopocentricFrame +from org.orekit.orbits import KeplerianOrbit +from org.orekit.frames import FramesFactory +from org.orekit.orbits import PositionAngle, EquinoctialOrbit, OrbitType +from org.orekit.propagation.analytical import KeplerianPropagator +from org.orekit.propagation import SpacecraftState +from org.orekit.propagation.numerical import NumericalPropagator +from org.orekit.propagation.events.handlers import EventHandler, PythonEventHandler +from org.orekit.time import AbsoluteDate +from org.orekit.time import TimeScalesFactory +from org.orekit.utils import Constants +from org.orekit.utils import IERSConventions, PVCoordinates +from org.orekit.propagation.events import ElevationDetector +from org.hipparchus.ode.events import Action +from org.orekit.forces.gravity.potential import GravityFieldFactory +from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel +from org.orekit.propagation.analytical import BrouwerLyddanePropagator +from org.hipparchus.util import FastMath, MathUtils + +from orekit import JArray_double + +import unittest +import sys +import math + + +class BrouwerLyddanePropagatorTest(unittest.TestCase): + + def setUp(self): + setup_orekit_curdir("resources") + self.provider = GravityFieldFactory.getNormalizedProvider(5, 0) + + def test_sameDateCartesian(self): + # Definition of initial conditions with position and velocity + # ------------------------------------------------------------ + # e = 0.04152500499523033 and i = 1.705015527659039 + + initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.) + position = Vector3D(3220103., 69623., 6149822.) + velocity = Vector3D(6414.7, -2006., -3180.) + + initialOrbit = EquinoctialOrbit(PVCoordinates(position, velocity), + FramesFactory.getEME2000(), initDate, self.provider.getMu()) + + # Extrapolation at the initial date + # --------------------------------- + extrapolator = BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(self.provider), + BrouwerLyddanePropagator.M2) + + finalOrbit = extrapolator.propagate(initDate) + + # positions velocity and semi major axis match perfectly + self.assertAlmostEquals(0.0, Vector3D.distance(initialOrbit.getPVCoordinates().getPosition(), + finalOrbit.getPVCoordinates().getPosition()), delta=1.0e-8) + + self.assertAlmostEquals(0.0, Vector3D.distance(initialOrbit.getPVCoordinates().getVelocity(), + finalOrbit.getPVCoordinates().getVelocity()), delta= 1.0e-11) + self.assertAlmostEquals(0.0, finalOrbit.getA() - initialOrbit.getA(), delta=0.0) + + def test_compareToNumericalPropagation(self): + + inertialFrame = FramesFactory.getEME2000() + initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.) + timeshift = 60000.0 + + # Initial orbit + a = 24396159.0 # semi major axis in meters + e = 0.01 # eccentricity + i = FastMath.toRadians(7.0) # inclination + omega = FastMath.toRadians(180.0) # perigee argument + raan = FastMath.toRadians(261.0) # right ascention of ascending node + lM = 0.0 # mean anomaly + initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.TRUE, + inertialFrame, initDate, self.provider.getMu()) + # Initial state definition + initialState = SpacecraftState(initialOrbit) + + #_______________________________________________________________________________________________ + # SET UP A REFERENCE NUMERICAL PROPAGATION + #_______________________________________________________________________________________________ + + # Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000 + minStep = 0.001 + maxstep = 1000.0 + positionTolerance = 10.0 + propagationType = OrbitType.KEPLERIAN + tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType) + integrator = DormandPrince853Integrator(minStep, maxstep, + JArray_double.cast_(tolerances[0]), + JArray_double.cast_(tolerances[1])) + + # Numerical Propagator + NumPropagator = NumericalPropagator(integrator) + NumPropagator.setOrbitType(propagationType) + + holmesFeatherstone = HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, True), self.provider) + NumPropagator.addForceModel(holmesFeatherstone) + + # Set up initial state in the propagator + NumPropagator.setInitialState(initialState) + + # Extrapolate from the initial to the date + NumFinalState = NumPropagator.propagate(initDate.shiftedBy(timeshift)) + NumOrbit = KeplerianOrbit.cast_(OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit())) + + #_______________________________________________________________________________________________ + # SET UP A BROUWER LYDDANE PROPAGATION + #_______________________________________________________________________________________________ + + BLextrapolator = BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(self.provider), + BrouwerLyddanePropagator.M2) + + BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift)) + BLOrbit = KeplerianOrbit.cast_(OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit())) + + self.assertAlmostEquals(NumOrbit.getA(), BLOrbit.getA(), delta=0.072) + self.assertAlmostEquals(NumOrbit.getE(), BLOrbit.getE(), delta=0.00000028) + self.assertAlmostEquals(NumOrbit.getI(), BLOrbit.getI(), delta=0.000004) + self.assertAlmostEquals(MathUtils.normalizeAngle(NumOrbit.getPerigeeArgument(), FastMath.PI), + MathUtils.normalizeAngle(BLOrbit.getPerigeeArgument(), FastMath.PI), delta=0.119) + self.assertAlmostEquals(MathUtils.normalizeAngle(NumOrbit.getRightAscensionOfAscendingNode(), FastMath.PI), + MathUtils.normalizeAngle(BLOrbit.getRightAscensionOfAscendingNode(), FastMath.PI), delta=0.000072) + self.assertAlmostEquals(MathUtils.normalizeAngle(NumOrbit.getTrueAnomaly(), FastMath.PI), + MathUtils.normalizeAngle(BLOrbit.getTrueAnomaly(), FastMath.PI), delta=0.12) + + +if __name__ == '__main__': + suite = unittest.TestLoader().loadTestsFromTestCase(BrouwerLyddanePropagatorTest) + ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() + sys.exit(ret) + diff --git a/python_files/test/Context.py b/python_files/test/Context.py new file mode 100644 index 0000000000000000000000000000000000000000..2950981bab671f3f89e53a4c52e81024d20f7a54 --- /dev/null +++ b/python_files/test/Context.py @@ -0,0 +1,162 @@ +# -*- coding: utf-8 -*- + +""" +/* Copyright 2002-2018 CS Systèmes d'Information + * Licensed to CS Systèmes 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. + */ + + Python version translated from Java by Petrus Hyvönen, SSC 2018 + +""" + +import math +import sys +import unittest + +# Python orekit specifics +import orekit +orekit.initVM() + + + +from org.hipparchus.geometry.euclidean.threed import Vector3D +from org.hipparchus.ode.nonstiff import DormandPrince853Integrator +from orekit import JArray_double +from org.orekit.data import DataProvidersManager, ZipJarCrawler +from org.orekit.propagation.integration import PythonAdditionalEquations +from org.orekit.forces.gravity.potential import GravityFieldFactory +from org.orekit.forces.gravity.potential import SHMFormatReader +from java.io import File + + +from org.hipparchus.geometry.euclidean.threed import Vector3D +from org.hipparchus.ode.nonstiff import DormandPrince853Integrator + +from org.orekit.frames import FramesFactory +from org.orekit.orbits import EquinoctialOrbit +from org.orekit.orbits import OrbitType +from org.orekit.propagation import SpacecraftState +from org.orekit.propagation.numerical import NumericalPropagator +from org.orekit.propagation.semianalytical.dsst import DSSTPropagator +from org.orekit.time import AbsoluteDate +from org.orekit.utils import PVCoordinates + +from org.orekit.data import DataProvidersManager, ZipJarCrawler +from java.io import File + +from org.hipparchus.geometry.euclidean.threed import Rotation +from org.hipparchus.geometry.euclidean.threed import Vector3D +from org.hipparchus.util import Decimal64Field +from org.hipparchus.util import FastMath + +from org.orekit.bodies import CelestialBodyFactory + +from org.orekit.frames import FramesFactory +from org.orekit.orbits import KeplerianOrbit +from org.orekit.orbits import PositionAngle +from org.orekit.propagation import FieldSpacecraftState +from org.orekit.propagation import SpacecraftState +from org.orekit.propagation.analytical import KeplerianPropagator +from org.orekit.time import AbsoluteDate +from org.orekit.time import DateComponents +from org.orekit.time import FieldAbsoluteDate +from org.orekit.time import TimeComponents +from org.orekit.time import TimeScalesFactory +from org.orekit.utils import AngularCoordinates +from org.orekit.utils import PVCoordinates +from org.orekit.utils import PVCoordinatesProvider +from org.orekit.attitudes import CelestialBodyPointed, SpinStabilized, InertialProvider + +# import java.util.List; +# import java.util.Map; +# +# import org.hipparchus.geometry.euclidean.threed.Vector3D; +# import org.hipparchus.util.FastMath; +# import org.orekit.bodies.CelestialBody; +from org.orekit.bodies import GeodeticPoint +# import org.orekit.bodies.OneAxisEllipsoid; +from org.orekit.estimation.measurements import GroundStation +# import org.orekit.forces.drag.DragSensitive; +# import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider; +# import org.orekit.forces.radiation.RadiationSensitive; +from org.orekit.frames import TopocentricFrame +# import org.orekit.models.earth.displacement.StationDisplacement; +from org.orekit.orbits import CartesianOrbit +# import org.orekit.orbits.Orbit; +# import org.orekit.orbits.OrbitType; +# import org.orekit.orbits.PositionAngle; +from org.orekit.propagation.conversion import DormandPrince853IntegratorBuilder +from org.orekit.propagation.conversion import NumericalPropagatorBuilder +# import org.orekit.time.TimeScale; +# import org.orekit.time.UT1Scale; +from org.orekit.utils import IERSConventions +from org.orekit.utils import PVCoordinates + +class Context(): + + + # public class Context { + # public IERSConventions conventions; + # public OneAxisEllipsoid earth; + # public CelestialBody sun; + # public CelestialBody moon; + # public RadiationSensitive radiationSensitive; + # public DragSensitive dragSensitive; + # public NormalizedSphericalHarmonicsProvider gravity; + # public TimeScale utc; + # public UT1Scale ut1; + # public Orbit initialOrbit; + # public StationDisplacement[] displacements; + # public List<GroundStation> stations; + # // Stations for turn-around range + # // Map entry = master station + # // Map value = slave station associated + # public Map<GroundStation, GroundStation> TARstations; + + def createBuilder(self, orbitType, positionAngle, perfectStart, minStep, maxStep, dP, *forces): + + if perfectStart: + # orbit estimation will start from a perfect orbit + startOrbit = self.initialOrbit + else: + # orbit estimation will start from a wrong point + initialPosition = self.initialOrbit.getPVCoordinates().getPosition() + initialVelocity = self.initialOrbit.getPVCoordinates().getVelocity() + wrongPosition = self.initialPosition.add(Vector3D(1000.0, 0.0, 0.0)) + wrongVelocity = initialVelocity.add(Vector3D(0.0, 0.0, 0.01)) + startOrbit = CartesianOrbit(PVCoordinates(wrongPosition, wrongVelocity), + self.initialOrbit.getFrame(), + self.initialOrbit.getDate(), + self.initialOrbit.getMu()) + propagatorBuilder = NumericalPropagatorBuilder(orbitType.convertType(startOrbit), + DormandPrince853IntegratorBuilder(minStep, maxStep, dP), + positionAngle, dP) + + for force in forces: + propagatorBuilder.addForceModel(force.getForceModel(self)) # self?? + + return propagatorBuilder + + + def createStation(self, latitudeInDegrees, longitudeInDegrees, altitude, name): + gp = GeodeticPoint(FastMath.toRadians(latitudeInDegrees), + FastMath.toRadians(longitudeInDegrees), + altitude) + return GroundStation(TopocentricFrame(self.earth, gp, name), + self.ut1.getEOPHistory(), self.displacements) + + + diff --git a/python_files/test/EstimationTestUtils.py b/python_files/test/EstimationTestUtils.py new file mode 100644 index 0000000000000000000000000000000000000000..c3aee8aecfdb5e753479d2af6b78c45dc1f16af1 --- /dev/null +++ b/python_files/test/EstimationTestUtils.py @@ -0,0 +1,578 @@ +""" +/* Copyright 2002-2018 CS Systèmes d'Information + * Licensed to CS Systèmes 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. + */ + +Python version translated from Java by Petrus Hyvönen, SSC 2018 + +""" + +import unittest +import sys + +# Python orekit specifics +import orekit + +orekit.initVM() + +#from orekit.pyhelpers import setup_orekit_curdir +#setup_orekit_curdir() + +from orekit import JArray_double +from org.orekit.data import DataProvidersManager, ZipJarCrawler, DataContext, DirectoryCrawler +from java.util import Arrays, HashMap +from java.io import File + +# from org.hipparchus.RealFieldElement; +# from org.hipparchus.geometry.euclidean.threed.FieldRotation; +# from org.hipparchus.geometry.euclidean.threed.FieldVector3D; +# from org.hipparchus.geometry.euclidean.threed.Rotation; +# from org.hipparchus.geometry.euclidean.threed.RotationConvention; +# from org.hipparchus.geometry.euclidean.threed.Vector3D; +# from org.hipparchus.linear.MatrixUtils; +# from org.hipparchus.linear.RealMatrix; +# from org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer.Optimum; +from org.hipparchus.util import FastMath +# from org.junit.Assert; +# from org.orekit.Utils; +from org.orekit.bodies import CelestialBodyFactory +from org.orekit.bodies import OneAxisEllipsoid +# from org.orekit.data.DataProvidersManager; +# from org.orekit.estimation.leastsquares.BatchLSEstimator; +# from org.orekit.estimation.measurements.EstimatedMeasurement; +from org.orekit.estimation.measurements import GroundStation +# from org.orekit.estimation.measurements.MeasurementCreator; +# from org.orekit.estimation.measurements.ObservedMeasurement; +# from org.orekit.estimation.sequential.KalmanEstimator; +from org.orekit.forces.drag import IsotropicDrag +from org.orekit.forces.gravity.potential import AstronomicalAmplitudeReader +from org.orekit.forces.gravity.potential import FESCHatEpsilonReader +from org.orekit.forces.gravity.potential import GRGSFormatReader +from org.orekit.forces.gravity.potential import GravityFieldFactory; +from org.orekit.forces.gravity.potential import OceanLoadDeformationCoefficients +from org.orekit.forces.radiation import IsotropicRadiationClassicalConvention +# from org.orekit.frames.EOPHistory; +# from org.orekit.frames.FieldTransform; +# from org.orekit.frames.Frame; +from org.orekit.frames import FramesFactory +# from org.orekit.frames.Transform; +# from org.orekit.frames.TransformProvider; +from org.orekit.models.earth.displacement import StationDisplacement +from org.orekit.models.earth.displacement import TidalDisplacement; +from org.orekit.orbits import KeplerianOrbit +# from org.orekit.orbits.Orbit; +from org.orekit.orbits import PositionAngle +# from org.orekit.propagation.Propagator; +# from org.orekit.propagation.conversion.PropagatorBuilder; +# from org.orekit.propagation.numerical.NumericalPropagator; +from org.orekit.time import AbsoluteDate +# from org.orekit.time.FieldAbsoluteDate; +from org.orekit.time import TimeScalesFactory +from org.orekit.utils import Constants +from org.orekit.utils import IERSConventions +# from org.orekit.utils.PVCoordinates; +# from org.orekit.utils.ParameterDriver; +from org.orekit.propagation import Propagator + +from Context import Context + + +class EstimationTestUtils(): + def eccentricContext(self, dataRoot): + DM = DataContext.getDefault().getDataProvidersManager() + datafile = File(dataRoot) + if not datafile.exists(): + print('File :', datafile.absolutePath, ' not found') + + crawler = DirectoryCrawler(datafile) + #DM.clearProviders() + DM.addProvider(crawler) + #DataProvidersManager.OREKIT_DATA_PATH = dataRoot + context = Context() + + context.conventions = IERSConventions.IERS_2010 + context.earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, + Constants.WGS84_EARTH_FLATTENING, + FramesFactory.getITRF(context.conventions, True)) + context.sun = CelestialBodyFactory.getSun() + context.moon = CelestialBodyFactory.getMoon() + context.radiationSensitive = IsotropicRadiationClassicalConvention(2.0, 0.2, 0.8) + + context.dragSensitive = IsotropicDrag(2.0, 1.2) + eopHistory = FramesFactory.getEOPHistory(context.conventions, True) + context.utc = TimeScalesFactory.getUTC() + context.ut1 = TimeScalesFactory.getUT1(eopHistory) + context.displacements = [TidalDisplacement(Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS, + # StationDisplacement[] { #NOT UNDERSTOOD, WHAT DOES THIS DO?? + Constants.JPL_SSD_SUN_EARTH_PLUS_MOON_MASS_RATIO, + Constants.JPL_SSD_EARTH_MOON_MASS_RATIO, + context.sun, context.moon, + context.conventions, False)] + + GravityFieldFactory.addPotentialCoefficientsReader(GRGSFormatReader("grim4s4_gr", True)) + aaReader = AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0) + + DataContext.getDefault().getDataProvidersManager().feed(aaReader.getSupportedNames(), aaReader) + map = aaReader.getAstronomicalAmplitudesMap() + GravityFieldFactory.addOceanTidesReader(FESCHatEpsilonReader("fes2004-7x7.dat", + 0.01, FastMath.toRadians(1.0), + OceanLoadDeformationCoefficients.IERS_2010, + map)) + context.gravity = GravityFieldFactory.getNormalizedProvider(20, 20) + context.initialOrbit = KeplerianOrbit(15000000.0, 0.125, 1.25, + 0.250, 1.375, 0.0625, PositionAngle.TRUE, + FramesFactory.getEME2000(), + AbsoluteDate(2000, 2, 24, 11, 35, 47.0, context.utc), + context.gravity.getMu()) + + context.stations = Arrays.asList([context.createStation(-53.05388, -75.01551, 1750.0, "Isla Desolación"), + context.createStation(62.29639, -7.01250, 880.0, "Slættaratindur")]) + + # Turn-around range stations + # Map entry = master station + # Map value = slave station associated + context.TARstations = HashMap().of_(GroundStation, GroundStation) + + context.TARstations.put(context.createStation(-53.05388, -75.01551, 1750.0, "Isla Desolación"), + context.createStation(-54.815833, -68.317778, 6.0, "Ushuaïa")) + + context.TARstations.put(context.createStation( 62.29639, -7.01250, 880.0, "Slættaratindur"), + context.createStation( 61.405833, -6.705278, 470.0, "Sumba")) + + return context + + def createPropagator(self, initialOrbit, propagatorBuilder): + + # override orbital parameters + orbitArray = JArray_double(6) + propagatorBuilder.getOrbitType().mapOrbitToArray(initialOrbit, + propagatorBuilder.getPositionAngle(), + orbitArray, None); + for i in range(0, len(orbitArray)-1): + propagatorBuilder.getOrbitalParametersDrivers().getDrivers().get(i).setValue(orbitArray[i]) + + return propagatorBuilder.buildPropagator(propagatorBuilder.getSelectedNormalizedParameters()) + + def createMeasurements(propagator, creator, startPeriod, endPeriod, step): + Propagator.cast_(propagator).setStepHandler(step, creator) + period = propagator.getInitialState().getKeplerianPeriod() + start = propagator.getInitialState().getDate().shiftedBy(startPeriod * period) + end = propagator.getInitialState().getDate().shiftedBy(endPeriod * period) + propagator.propagate(start, end) + + measurements = creator.getMeasurements() + + for measurement in measurements: + for driver in measurement.getParametersDrivers(): + if driver.getReferenceDate is None: + driver.setReferenceDate(propagator.getInitialState().getDate()) + + return measurements + +""" +/** Utility class for orbit determination tests. */ +public class EstimationTestUtils { + + public static Context eccentricContext(final String dataRoot) { + + Utils.setDataRoot(dataRoot); + Context context = new Context(); + context.conventions = IERSConventions.IERS_2010; + context.earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, + Constants.WGS84_EARTH_FLATTENING, + FramesFactory.getITRF(context.conventions, true)); + context.sun = CelestialBodyFactory.getSun(); + context.moon = CelestialBodyFactory.getMoon(); + context.radiationSensitive = new IsotropicRadiationClassicalConvention(2.0, 0.2, 0.8); + context.dragSensitive = new IsotropicDrag(2.0, 1.2); + final EOPHistory eopHistory = FramesFactory.getEOPHistory(context.conventions, true); + context.utc = TimeScalesFactory.getUTC(); + context.ut1 = TimeScalesFactory.getUT1(eopHistory); + context.displacements = new StationDisplacement[] { + new TidalDisplacement(Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS, + Constants.JPL_SSD_SUN_EARTH_PLUS_MOON_MASS_RATIO, + Constants.JPL_SSD_EARTH_MOON_MASS_RATIO, + context.sun, context.moon, + context.conventions, false) + }; + GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true)); + AstronomicalAmplitudeReader aaReader = + new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0); + DataProvidersManager.getInstance().feed(aaReader.getSupportedNames(), aaReader); + Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap(); + GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat", + 0.01, FastMath.toRadians(1.0), + OceanLoadDeformationCoefficients.IERS_2010, + map)); + context.gravity = GravityFieldFactory.getNormalizedProvider(20, 20); + context.initialOrbit = new KeplerianOrbit(15000000.0, 0.125, 1.25, + 0.250, 1.375, 0.0625, PositionAngle.TRUE, + FramesFactory.getEME2000(), + new AbsoluteDate(2000, 2, 24, 11, 35, 47.0, context.utc), + context.gravity.getMu()); + + context.stations = Arrays.asList(//context.createStation(-18.59146, -173.98363, 76.0, "Leimatu`a"), + context.createStation(-53.05388, -75.01551, 1750.0, "Isla Desolación"), + context.createStation( 62.29639, -7.01250, 880.0, "Slættaratindur") + //context.createStation( -4.01583, 103.12833, 3173.0, "Gunung Dempo") + ); + + // Turn-around range stations + // Map entry = master station + // Map value = slave station associated + context.TARstations = new HashMap<GroundStation, GroundStation>(); + + context.TARstations.put(context.createStation(-53.05388, -75.01551, 1750.0, "Isla Desolación"), + context.createStation(-54.815833, -68.317778, 6.0, "Ushuaïa")); + + context.TARstations.put(context.createStation( 62.29639, -7.01250, 880.0, "Slættaratindur"), + context.createStation( 61.405833, -6.705278, 470.0, "Sumba")); + + return context; + + } + + public static Context geoStationnaryContext(final String dataRoot) { + + Utils.setDataRoot(dataRoot); + Context context = new Context(); + context.conventions = IERSConventions.IERS_2010; + context.utc = TimeScalesFactory.getUTC(); + context.ut1 = TimeScalesFactory.getUT1(context.conventions, true); + context.displacements = new StationDisplacement[0]; + String Myframename = "MyEarthFrame"; + final AbsoluteDate datedef = new AbsoluteDate(2000, 1, 1, 12, 0, 0.0, context.utc); + final double omega = Constants.WGS84_EARTH_ANGULAR_VELOCITY; + final Vector3D rotationRate = new Vector3D(0.0, 0.0, omega); + + TransformProvider MyEarthFrame = new TransformProvider() { + private static final long serialVersionUID = 1L; + public Transform getTransform(final AbsoluteDate date) { + final double rotationduration = date.durationFrom(datedef); + final Vector3D alpharot = new Vector3D(rotationduration, rotationRate); + final Rotation rotation = new Rotation(Vector3D.PLUS_K, -alpharot.getZ(), + RotationConvention.VECTOR_OPERATOR); + return new Transform(date, rotation, rotationRate); + } + public <T extends RealFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) { + final T rotationduration = date.durationFrom(datedef); + final FieldVector3D<T> alpharot = new FieldVector3D<>(rotationduration, rotationRate); + final FieldRotation<T> rotation = new FieldRotation<>(FieldVector3D.getPlusK(date.getField()), + alpharot.getZ().negate(), + RotationConvention.VECTOR_OPERATOR); + return new FieldTransform<>(date, rotation, new FieldVector3D<>(date.getField(), rotationRate)); + } + }; + Frame FrameTest = new Frame(FramesFactory.getEME2000(), MyEarthFrame, Myframename, true); + + // Earth is spherical, rotating in one sidereal day + context.earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 0.0, FrameTest); + context.sun = CelestialBodyFactory.getSun(); + context.moon = CelestialBodyFactory.getMoon(); + context.radiationSensitive = new IsotropicRadiationClassicalConvention(2.0, 0.2, 0.8); + context.dragSensitive = new IsotropicDrag(2.0, 1.2); + GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true)); + AstronomicalAmplitudeReader aaReader = + new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0); + DataProvidersManager.getInstance().feed(aaReader.getSupportedNames(), aaReader); + Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap(); + GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat", + 0.01, FastMath.toRadians(1.0), + OceanLoadDeformationCoefficients.IERS_2010, + map)); + context.gravity = GravityFieldFactory.getNormalizedProvider(20, 20); + + // semimajor axis for a geostationnary satellite + double da = FastMath.cbrt(context.gravity.getMu() / (omega * omega)); + + //context.stations = Arrays.asList(context.createStation( 0.0, 0.0, 0.0, "Lat0_Long0"), + // context.createStation( 62.29639, -7.01250, 880.0, "Slættaratindur") + // ); + context.stations = Arrays.asList(context.createStation(0.0, 0.0, 0.0, "Lat0_Long0") ); + + // Station position & velocity in EME2000 + final Vector3D geovelocity = new Vector3D (0., 0., 0.); + + // Compute the frames transformation from station frame to EME2000 + Transform topoToEME = + context.stations.get(0).getBaseFrame().getTransformTo(FramesFactory.getEME2000(), new AbsoluteDate(2000, 1, 1, 12, 0, 0.0, context.utc)); + + // Station position in EME2000 at reference date + Vector3D stationPositionEME = topoToEME.transformPosition(Vector3D.ZERO); + + // Satellite position and velocity in Station Frame + final Vector3D sat_pos = new Vector3D(0., 0., da-stationPositionEME.getNorm()); + final Vector3D acceleration = new Vector3D(-context.gravity.getMu(), sat_pos); + final PVCoordinates pv_sat_topo = new PVCoordinates(sat_pos, geovelocity, acceleration); + + // satellite position in EME2000 + final PVCoordinates pv_sat_iner = topoToEME.transformPVCoordinates(pv_sat_topo); + + // Geo-stationary Satellite Orbit, tightly above the station (l0-L0) + context.initialOrbit = new KeplerianOrbit(pv_sat_iner, + FramesFactory.getEME2000(), + new AbsoluteDate(2000, 1, 1, 12, 0, 0.0, context.utc), + context.gravity.getMu()); + + context.stations = Arrays.asList(context.createStation(10.0, 45.0, 0.0, "Lat10_Long45") ); + + // Turn-around range stations + // Map entry = master station + // Map value = slave station associated + context.TARstations = new HashMap<GroundStation, GroundStation>(); + + context.TARstations.put(context.createStation( 41.977, 13.600, 671.354, "Fucino"), + context.createStation( 43.604, 1.444, 263.0 , "Toulouse")); + + context.TARstations.put(context.createStation( 49.867, 8.65 , 144.0 , "Darmstadt"), + context.createStation( -25.885, 27.707, 1566.633, "Pretoria")); + + return context; + + } + + public static Propagator createPropagator(final Orbit initialOrbit, + final PropagatorBuilder propagatorBuilder) + { + + // override orbital parameters + double[] orbitArray = new double[6]; + propagatorBuilder.getOrbitType().mapOrbitToArray(initialOrbit, + propagatorBuilder.getPositionAngle(), + orbitArray, null); + for (int i = 0; i < orbitArray.length; ++i) { + propagatorBuilder.getOrbitalParametersDrivers().getDrivers().get(i).setValue(orbitArray[i]); + } + + return propagatorBuilder.buildPropagator(propagatorBuilder.getSelectedNormalizedParameters()); + + } + + public static List<ObservedMeasurement<?>> createMeasurements(final Propagator propagator, + final MeasurementCreator creator, + final double startPeriod, final double endPeriod, + final double step) + { + + propagator.setMasterMode(step, creator); + final double period = propagator.getInitialState().getKeplerianPeriod(); + final AbsoluteDate start = propagator.getInitialState().getDate().shiftedBy(startPeriod * period); + final AbsoluteDate end = propagator.getInitialState().getDate().shiftedBy(endPeriod * period); + propagator.propagate(start, end); + + final List<ObservedMeasurement<?>> measurements = creator.getMeasurements(); + + for (final ObservedMeasurement<?> measurement : measurements) { + for (final ParameterDriver driver : measurement.getParametersDrivers()) { + if (driver.getReferenceDate() == null) { + driver.setReferenceDate(propagator.getInitialState().getDate()); + } + } + } + + return measurements; + + } + + /** + * Checker for batch LS estimator validation + * @param context Context used for the test + * @param estimator Batch LS estimator + * @param iterations Number of iterations expected + * @param evaluations Number of evaluations expected + * @param expectedRMS Expected RMS value + * @param rmsEps Tolerance on expected RMS + * @param expectedMax Expected weighted residual maximum + * @param maxEps Tolerance on weighted residual maximum + * @param expectedDeltaPos Expected position difference between estimated orbit and initial orbit + * @param posEps Tolerance on expected position difference + * @param expectedDeltaVel Expected velocity difference between estimated orbit and initial orbit + * @param velEps Tolerance on expected velocity difference + */ + public static void checkFit(final Context context, final BatchLSEstimator estimator, + final int iterations, final int evaluations, + final double expectedRMS, final double rmsEps, + final double expectedMax, final double maxEps, + final double expectedDeltaPos, final double posEps, + final double expectedDeltaVel, final double velEps) + { + + final Orbit estimatedOrbit = estimator.estimate()[0].getInitialState().getOrbit(); + final Vector3D estimatedPosition = estimatedOrbit.getPVCoordinates().getPosition(); + final Vector3D estimatedVelocity = estimatedOrbit.getPVCoordinates().getVelocity(); + + Assert.assertEquals(iterations, estimator.getIterationsCount()); + Assert.assertEquals(evaluations, estimator.getEvaluationsCount()); + Optimum optimum = estimator.getOptimum(); + Assert.assertEquals(iterations, optimum.getIterations()); + Assert.assertEquals(evaluations, optimum.getEvaluations()); + + int k = 0; + double sum = 0; + double max = 0; + for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry : + estimator.getLastEstimations().entrySet()) { + final ObservedMeasurement<?> m = entry.getKey(); + final EstimatedMeasurement<?> e = entry.getValue(); + final double[] weight = m.getBaseWeight(); + final double[] sigma = m.getTheoreticalStandardDeviation(); + final double[] observed = m.getObservedValue(); + final double[] theoretical = e.getEstimatedValue(); + for (int i = 0; i < m.getDimension(); ++i) { + final double weightedResidual = weight[i] * (theoretical[i] - observed[i]) / sigma[i]; + ++k; + sum += weightedResidual * weightedResidual; + max = FastMath.max(max, FastMath.abs(weightedResidual)); + } + } + + final double rms = FastMath.sqrt(sum / k); + final double deltaPos = Vector3D.distance(context.initialOrbit.getPVCoordinates().getPosition(), estimatedPosition); + final double deltaVel = Vector3D.distance(context.initialOrbit.getPVCoordinates().getVelocity(), estimatedVelocity); + Assert.assertEquals(expectedRMS, + rms, + rmsEps); + Assert.assertEquals(expectedMax, + max, + maxEps); + Assert.assertEquals(expectedDeltaPos, + deltaPos, + posEps); + Assert.assertEquals(expectedDeltaVel, + deltaVel, + velEps); + + } + + /** + * Checker for Kalman estimator validation + * @param context context used for the test + * @param kalman Kalman filter + * @param measurements List of observed measurements to be processed by the Kalman + * @param refOrbit Reference orbits at last measurement date + * @param expectedDeltaPos Expected position difference between estimated orbit and reference orbit + * @param posEps Tolerance on expected position difference + * @param expectedDeltaVel Expected velocity difference between estimated orbit and reference orbit + * @param velEps Tolerance on expected velocity difference + * @param expectedSigmasPos Expected values for covariance matrix on position + * @param sigmaPosEps Tolerance on expected covariance matrix on position + * @param expectedSigmasVel Expected values for covariance matrix on velocity + * @param sigmaVelEps Tolerance on expected covariance matrix on velocity + */ + public static void checkKalmanFit(final Context context, final KalmanEstimator kalman, + final List<ObservedMeasurement<?>> measurements, + final Orbit refOrbit, final PositionAngle positionAngle, + final double expectedDeltaPos, final double posEps, + final double expectedDeltaVel, final double velEps, + final double[] expectedSigmasPos,final double sigmaPosEps, + final double[] expectedSigmasVel,final double sigmaVelEps) + { + checkKalmanFit(context, kalman, measurements, + new Orbit[] { refOrbit }, + new PositionAngle[] { positionAngle }, + new double[] { expectedDeltaPos }, new double[] { posEps }, + new double[] { expectedDeltaVel }, new double[] { velEps }, + new double[][] { expectedSigmasPos }, new double[] { sigmaPosEps }, + new double[][] { expectedSigmasVel }, new double[] { sigmaVelEps }); + } + + /** + * Checker for Kalman estimator validation + * @param context context used for the test + * @param kalman Kalman filter + * @param measurements List of observed measurements to be processed by the Kalman + * @param refOrbit Reference orbits at last measurement date + * @param expectedDeltaPos Expected position difference between estimated orbit and reference orbits + * @param posEps Tolerance on expected position difference + * @param expectedDeltaVel Expected velocity difference between estimated orbit and reference orbits + * @param velEps Tolerance on expected velocity difference + * @param expectedSigmasPos Expected values for covariance matrix on position + * @param sigmaPosEps Tolerance on expected covariance matrix on position + * @param expectedSigmasVel Expected values for covariance matrix on velocity + * @param sigmaVelEps Tolerance on expected covariance matrix on velocity + */ + public static void checkKalmanFit(final Context context, final KalmanEstimator kalman, + final List<ObservedMeasurement<?>> measurements, + final Orbit[] refOrbit, final PositionAngle[] positionAngle, + final double[] expectedDeltaPos, final double[] posEps, + final double[] expectedDeltaVel, final double []velEps, + final double[][] expectedSigmasPos,final double[] sigmaPosEps, + final double[][] expectedSigmasVel,final double[] sigmaVelEps) + { + + // Add the measurements to the Kalman filter + NumericalPropagator[] estimated = kalman.processMeasurements(measurements); + + // Check the number of measurements processed by the filter + Assert.assertEquals(measurements.size(), kalman.getCurrentMeasurementNumber()); + + for (int k = 0; k < refOrbit.length; ++k) { + // Get the last estimation + final Orbit estimatedOrbit = estimated[k].getInitialState().getOrbit(); + final Vector3D estimatedPosition = estimatedOrbit.getPVCoordinates().getPosition(); + final Vector3D estimatedVelocity = estimatedOrbit.getPVCoordinates().getVelocity(); + + // Get the last covariance matrix estimation + final RealMatrix estimatedP = kalman.getPhysicalEstimatedCovarianceMatrix(); + + // Convert the orbital part to Cartesian formalism + // Assuming all 6 orbital parameters are estimated by the filter + final double[][] dCdY = new double[6][6]; + estimatedOrbit.getJacobianWrtParameters(positionAngle[k], dCdY); + final RealMatrix Jacobian = MatrixUtils.createRealMatrix(dCdY); + final RealMatrix estimatedCartesianP = + Jacobian. + multiply(estimatedP.getSubMatrix(0, 5, 0, 5)). + multiply(Jacobian.transpose()); + + // Get the final sigmas (ie.sqrt of the diagonal of the Cartesian orbital covariance matrix) + final double[] sigmas = new double[6]; + for (int i = 0; i < 6; i++) { + sigmas[i] = FastMath.sqrt(estimatedCartesianP.getEntry(i, i)); + } +// // FIXME: debug print values +// final double dPos = Vector3D.distance(refOrbit[k].getPVCoordinates().getPosition(), estimatedPosition); +// final double dVel = Vector3D.distance(refOrbit[k].getPVCoordinates().getVelocity(), estimatedVelocity); +// System.out.println("Nb Meas = " + kalman.getCurrentMeasurementNumber()); +// System.out.println("dPos = " + dPos + " m"); +// System.out.println("dVel = " + dVel + " m/s"); +// System.out.println("sigmas = " + sigmas[0] + " " +// + sigmas[1] + " " +// + sigmas[2] + " " +// + sigmas[3] + " " +// + sigmas[4] + " " +// + sigmas[5]); +// //debug + + // Check the final orbit estimation & PV sigmas + final double deltaPosK = Vector3D.distance(refOrbit[k].getPVCoordinates().getPosition(), estimatedPosition); + final double deltaVelK = Vector3D.distance(refOrbit[k].getPVCoordinates().getVelocity(), estimatedVelocity); + Assert.assertEquals(expectedDeltaPos[k], deltaPosK, posEps[k]); + Assert.assertEquals(expectedDeltaVel[k], deltaVelK, velEps[k]); + + for (int i = 0; i < 3; i++) { + Assert.assertEquals(expectedSigmasPos[k][i], sigmas[i], sigmaPosEps[k]); + Assert.assertEquals(expectedSigmasVel[k][i], sigmas[i+3], sigmaVelEps[k]); + } + } + } + +}""" + +if __name__ == '__main__': + a = EstimationTestUtils() + Context = a.eccentricContext("../test/resources") + print(Context) diff --git a/python_files/test/EventDetectorTest.py b/python_files/test/EventDetectorTest.py index 5d2914f775fc6b1bbd3917a15984be996d14a6d0..334152811e5809d530e8133a91080f4d35ee862b 100644 --- a/python_files/test/EventDetectorTest.py +++ b/python_files/test/EventDetectorTest.py @@ -48,7 +48,7 @@ import sys from orekit.pyhelpers import setup_orekit_curdir -setup_orekit_curdir() +setup_orekit_curdir("resources") class MyElevationDetector(PythonEventDetector): diff --git a/python_files/test/EventHandlerTest.py b/python_files/test/EventHandlerTest.py index 4589183ccf05545a0e67a94af23fce2cfb692c5c..9580237d8d0c2a4e46424432a307d0d670fa0dd1 100644 --- a/python_files/test/EventHandlerTest.py +++ b/python_files/test/EventHandlerTest.py @@ -46,7 +46,7 @@ import sys #%% Setup Orekit from orekit.pyhelpers import setup_orekit_curdir -setup_orekit_curdir() +setup_orekit_curdir("resources") class EventHandlerTest(unittest.TestCase): @@ -82,7 +82,7 @@ class EventHandlerTest(unittest.TestCase): class myContinueOnEvent(PythonEventHandler): - def init(self, initialstate, target): + def init(self, initialstate, target, detector): pass def eventOccurred(self, s, T, increasing): diff --git a/python_files/test/FieldAdditionalEquationsTest.py b/python_files/test/FieldAdditionalEquationsTest.py index afa7f373f46425ee2a3510f8c09ca34b7edf3d5e..440507e84cbd88bc38afe147a3a06b173e4707b4 100644 --- a/python_files/test/FieldAdditionalEquationsTest.py +++ b/python_files/test/FieldAdditionalEquationsTest.py @@ -32,7 +32,7 @@ import orekit orekit.initVM() from orekit import JArray_double -from org.orekit.data import DataProvidersManager, ZipJarCrawler +from org.orekit.data import DataProvidersManager, ZipJarCrawler, DataContext, DirectoryCrawler from org.orekit.propagation.integration import PythonFieldAdditionalEquations from org.orekit.forces.gravity.potential import GravityFieldFactory from org.orekit.forces.gravity.potential import SHMFormatReader @@ -41,7 +41,7 @@ from java.lang import System from org.hipparchus.geometry.euclidean.threed import Vector3D -from org.hipparchus import RealFieldElement +from org.hipparchus import CalculusFieldElement from org.hipparchus.ode.nonstiff import DormandPrince853Integrator from org.hipparchus.ode.nonstiff import DormandPrince853FieldIntegrator @@ -76,7 +76,7 @@ class InitCheckerEquations(PythonFieldAdditionalEquations): # implements Additi # Part of AdditionalEquations interface def computeDerivatives(self, s, pDot): - zerovalue = RealFieldElement.cast_(s.getDate().getField().getZero()) # Not completely sure why this needs casting but becomes just an "Object" otherwise + zerovalue = CalculusFieldElement.cast_(s.getDate().getField().getZero()) # Not completely sure why this needs casting but becomes just an "Object" otherwise pDot[0] = zerovalue.add(1.5) return MathArrays.buildArray(s.getDate().getField(), 7) @@ -91,12 +91,12 @@ class InitCheckerEquations(PythonFieldAdditionalEquations): # implements Additi class FieldAdditionalEquationsTest(unittest.TestCase): def setUp(self): # Initialize the data sources - DM = DataProvidersManager.getInstance() - datafile = File('resources.zip') + DM = DataContext.getDefault().getDataProvidersManager() + datafile = File('resources') if not datafile.exists(): print('File :', datafile.absolutePath, ' not found') - crawler = ZipJarCrawler(datafile) + crawler = DirectoryCrawler(datafile) DM.clearProviders() DM.addProvider(crawler) System.setProperty(DataProvidersManager.OREKIT_DATA_PATH, 'potential/shm-format') @@ -131,7 +131,7 @@ class FieldAdditionalEquationsTest(unittest.TestCase): # action integrator = DormandPrince853FieldIntegrator(field, 0.001, 200.0, JArray_double.cast_(self.tolerance[0]), JArray_double.cast_(self.tolerance[1])) - integrator.setInitialStepSize(field.getZero().add(60.0)) + integrator.setInitialStepSize(60.0) propagatorNumerical = FieldNumericalPropagator(field, integrator) propagatorNumerical.setInitialState( FieldSpacecraftState(field, self.initialState).addAdditionalState(checker.getName(), diff --git a/python_files/test/GroundFieldOfViewDetectorTest.py b/python_files/test/GroundFieldOfViewDetectorTest.py index d909a1b09aa416610b10414735e45838f951c32e..b95e042dab5c51a5717da6478566b5150ca7a10b 100644 --- a/python_files/test/GroundFieldOfViewDetectorTest.py +++ b/python_files/test/GroundFieldOfViewDetectorTest.py @@ -34,7 +34,10 @@ from org.orekit.orbits import KeplerianOrbit, PositionAngle from org.orekit.utils import Constants from org.orekit.propagation.analytical import KeplerianPropagator from org.orekit.utils import IERSConventions -from org.orekit.propagation.events import ElevationDetector, EventsLogger, FieldOfView, GroundFieldOfViewDetector +from org.orekit.propagation.events import ElevationDetector, EventsLogger, GroundFieldOfViewDetector +from org.orekit.geometry.fov import FieldOfView +from org.orekit.geometry.fov import PolygonalFieldOfView +# from org.orekit.geometry.fov import PolygonalFieldOfView.DefiningConeType from org.hipparchus.geometry.euclidean.threed import Vector3D from math import radians @@ -44,7 +47,7 @@ import sys from orekit.pyhelpers import setup_orekit_curdir -setup_orekit_curdir() +setup_orekit_curdir("resources") class GroundFieldOfViewDetectorTest(unittest.TestCase): @@ -80,13 +83,15 @@ class GroundFieldOfViewDetectorTest(unittest.TestCase): # half width of 60 deg pointed along +Z in antenna frame # not a perfect small circle b/c FoV makes a polygon with great circles - fov = FieldOfView(Vector3D.PLUS_K, Vector3D.PLUS_I, math.pi / 3.0, 16, 0.0) + + fov = PolygonalFieldOfView(Vector3D.PLUS_K, PolygonalFieldOfView.DefiningConeType.INSIDE_CONE_TOUCHING_POLYGON_AT_EDGES_MIDDLE, + Vector3D.PLUS_I, math.pi / 3.0, 16, 0.0) # simple case for fixed pointing to be similar to elevation detector. # could define new frame with varying rotation for slewing antenna. fovDetector = GroundFieldOfViewDetector(topo, fov).withMaxCheck(5.0) self.assertEqual(topo, fovDetector.getFrame()) - self.assertEqual(fov, fovDetector.getFieldOfView()) + self.assertEqual(fov, fovDetector.getFOV()) logger = EventsLogger() prop = KeplerianPropagator(orbit) diff --git a/python_files/test/ImpulseManeuverTest.py b/python_files/test/ImpulseManeuverTest.py index 4ea62286afd91d09f420a23c1b176dfe4423b417..1eb26c2b25300f46593e750cf0312155fef65269 100644 --- a/python_files/test/ImpulseManeuverTest.py +++ b/python_files/test/ImpulseManeuverTest.py @@ -27,7 +27,7 @@ Python version translated from Java by Petrus Hyvönen, SSC 2014 import orekit orekit.initVM() from orekit.pyhelpers import setup_orekit_curdir -setup_orekit_curdir() +setup_orekit_curdir("resources") from org.hipparchus.geometry.euclidean.threed import Vector3D from org.hipparchus.util import FastMath diff --git a/python_files/test/InterSatDirectViewDetectorTest.py b/python_files/test/InterSatDirectViewDetectorTest.py index b4d9aa60ccf4fd792c929ce3afa1c57177357ad2..a42852574508736fba4fc073eb7cd3bbe2068086 100644 --- a/python_files/test/InterSatDirectViewDetectorTest.py +++ b/python_files/test/InterSatDirectViewDetectorTest.py @@ -30,7 +30,7 @@ import orekit orekit.initVM() from orekit.pyhelpers import setup_orekit_curdir -setup_orekit_curdir() # orekit-data.zip shall be in current dir +setup_orekit_curdir("resources") from org.orekit.propagation.events import EventsLogger from org.hipparchus.util import FastMath @@ -47,6 +47,7 @@ from org.orekit.bodies import OneAxisEllipsoid from org.orekit.frames import TopocentricFrame from org.orekit.orbits import CircularOrbit from org.orekit.propagation.analytical import KeplerianPropagator +from org.orekit.propagation import Propagator from org.orekit.propagation.events.handlers import EventHandler, PythonEventHandler from org.orekit.utils import IERSConventions @@ -55,7 +56,7 @@ from org.hipparchus.ode.events import Action class GrazingHandler(PythonEventHandler): - def init(self, initialstate, target): + def init(self, initialstate, target, detector): pass def eventOccurred(self, s, detector, increasing): @@ -68,7 +69,7 @@ class GrazingHandler(PythonEventHandler): grazingDate = s.getDate().shiftedBy(dt) pMaster = s.shiftedBy(dt).getPVCoordinates(frame).getPosition() - pSlave = detector.getSlave().getPVCoordinates(grazingDate, frame).getPosition() + pSlave = detector.getSecondary().getPVCoordinates(grazingDate, frame).getPosition() grazing = earth.getCartesianIntersectionPoint(Line(pMaster, pSlave, 1.0), pMaster, frame, grazingDate) @@ -119,14 +120,17 @@ class InterSatDirectViewDetectorTest(unittest.TestCase): def init(self, s0, t, step): pass - def handleStep(self, state, isLast): + def handleStep(self, state): pos1 = state.getPVCoordinates().getPosition() pos2 = o2.getPVCoordinates(state.getDate(), state.getFrame()).getPosition() assert Vector3D.distance(pos1, pos2) > 8100.0 assert Vector3D.distance(pos1, pos2) < 16400.0 - p.setMasterMode(10.0, StepHandler()) + def finish(self,s): + pass + + Propagator.cast_(p).setStepHandler(10.0, StepHandler()) p.propagate(o1.getDate().shiftedBy(o1.getKeplerianPeriod())) self.assertEqual(0, logger.getLoggedEvents().size()) diff --git a/python_files/test/IodGibbsTest.py b/python_files/test/IodGibbsTest.py new file mode 100644 index 0000000000000000000000000000000000000000..5a1341af12767d7400d6f0437819ceff34b4cc3e --- /dev/null +++ b/python_files/test/IodGibbsTest.py @@ -0,0 +1,132 @@ +# -*- coding: utf-8 -*- + +""" +/* Copyright 2002-2018 CS Systèmes d'Information + * Licensed to CS Systèmes 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. + */ + + Python version translated from Java by Petrus Hyvönen, SSC 2018 + +""" + +import math +import sys +import unittest + +# Python orekit specifics +import orekit +orekit.initVM() + +from org.orekit.data import DataProvidersManager, ZipJarCrawler, DataContext +from java.io import File + +# import java.util.List; +# +from org.hipparchus.geometry.euclidean.threed import Vector3D +# import org.junit.Assert; +# import org.junit.Test; +# import org.orekit.estimation.Context; +# from org.orekit.estimation import EstimationTestUtils +# import org.orekit.estimation.measurements.ObservedMeasurement; +# import org.orekit.estimation.measurements.PV; +# import org.orekit.estimation.measurements.PVMeasurementCreator; +# import org.orekit.frames.Frame; +from org.orekit.frames import FramesFactory +# import org.orekit.orbits.KeplerianOrbit; +from org.orekit.orbits import OrbitType +from org.orekit.orbits import PositionAngle +# import org.orekit.propagation.Propagator; +# import org.orekit.propagation.conversion.NumericalPropagatorBuilder; +from org.orekit.time import AbsoluteDate, TimeScalesFactory +# import org.orekit.time.TimeScalesFactory; +from org.orekit.estimation.measurements import ObservableSatellite, PV +from PVMeasurementCreator import PVMeasurementCreator +from org.orekit.estimation.iod import IodGibbs + +import pathlib, os +curdir = pathlib.Path(__file__).parent.resolve() + +from EstimationTestUtils import EstimationTestUtils + +class IodGibbsTest(unittest.TestCase): + + def testGibbs1(self): + context = EstimationTestUtils().eccentricContext(os.path.join(curdir, "resources")) + mu = context.initialOrbit.getMu() + frame = context.initialOrbit.getFrame() + + propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, True, 1.0e-6, 60.0, 0.001) + + # create perfect range measurements + propagator = EstimationTestUtils().createPropagator(context.initialOrbit, propagatorBuilder) + satellite = ObservableSatellite(0) + measurements = EstimationTestUtils.createMeasurements(propagator, PVMeasurementCreator(), 0.0, 1.0, 60.0) + + position1 = Vector3D(*[x for x in measurements[0].getObservedValue()[0:3]]) + pv1 = PV(measurements[0].getDate(), position1, Vector3D.ZERO, 0., 0., 1., satellite) + + position2 = Vector3D(*[x for x in measurements[1].getObservedValue()[0:3]]) + pv2 = PV(measurements[1].getDate(), position2, Vector3D.ZERO, 0., 0., 1., satellite) + + position3 = Vector3D(*[x for x in measurements[2].getObservedValue()[0:3]]) + pv3 = PV(measurements[2].getDate(), position3, Vector3D.ZERO, 0., 0., 1., satellite) + + # instantiate the IOD method + + gibbs = IodGibbs(mu) + orbit = gibbs.estimate(frame, pv1, pv2, pv3) + + self.assertAlmostEquals(context.initialOrbit.getA(), orbit.getA(), delta=1.0e-9 * context.initialOrbit.getA()) + self.assertAlmostEquals(context.initialOrbit.getE(), orbit.getE(), delta=1.0e-9 * context.initialOrbit.getA()) + self.assertAlmostEquals(context.initialOrbit.getI(), orbit.getI(), delta=1.0e-9 * context.initialOrbit.getA()) + + pass + + + def testGibbs2(self): + + # test extracted from "Fundamentals of astrodynamics & applications", D. Vallado, 3rd ed, chap Initial Orbit Determination, Exple 7-3, p457 + context = EstimationTestUtils().eccentricContext(os.path.join(curdir, "resources")) + mu = context.initialOrbit.getMu() + + # Initialization + gibbs = IodGibbs(mu) + + # Observation vector (EME2000) + posR1 = Vector3D(0.0, 0.0, 6378137.0) + posR2 = Vector3D(0.0, -4464696.0, -5102509.0) + posR3 = Vector3D(0.0, 5740323.0, 3189068.0) + + # epoch corresponding to the observation vector + dateRef = AbsoluteDate(2000, 1, 1, 0, 0, 0.0, TimeScalesFactory.getUTC()) + date2 = dateRef.shiftedBy(76.48) + date3 = dateRef.shiftedBy(153.04) + + # Reference result (cf. Vallado) + velR2 = Vector3D(0.0, 5531.148, -5191.806) + + # Gibbs IOD + orbit = gibbs.estimate(FramesFactory.getEME2000(), + posR1, dateRef, posR2, date2, posR3, date3) + + # test + self.assertAlmostEquals(0.0, orbit.getPVCoordinates().getVelocity().getNorm() - velR2.getNorm(), delta=1e-3) + + +if __name__ == '__main__': + suite = unittest.TestLoader().loadTestsFromTestCase(IodGibbsTest) + ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() + sys.exit(ret) \ No newline at end of file diff --git a/python_files/test/IodLaplaceTest.py b/python_files/test/IodLaplaceTest.py new file mode 100644 index 0000000000000000000000000000000000000000..ae9360426047422bd4c500d2c6ee77ad2ac387dc --- /dev/null +++ b/python_files/test/IodLaplaceTest.py @@ -0,0 +1,200 @@ +# -*- coding: utf-8 -*- + +""" +/* Copyright 2002-2018 CS Systèmes d'Information + * Licensed to CS Systèmes 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. + */ + + Python version translated from Java by Petrus Hyvönen, SSC 2022 + +""" + +import math +import sys +import unittest + +# Python orekit specifics +import orekit +orekit.initVM() + +from orekit.pyhelpers import setup_orekit_curdir + + +from org.orekit.data import DataProvidersManager, ZipJarCrawler, DataContext +from java.io import File + +# import java.util.List; +# +from org.hipparchus.geometry.euclidean.threed import Vector3D +# import org.junit.Assert; +# import org.junit.Test; +# import org.orekit.estimation.Context; +from org.hipparchus.util import FastMath + +# from org.orekit.estimation import EstimationTestUtils +# import org.orekit.estimation.measurements.ObservedMeasurement; +# import org.orekit.estimation.measurements.PV; +# import org.orekit.estimation.measurements.PVMeasurementCreator; +# import org.orekit.frames.Frame; +from org.orekit.frames import FramesFactory, TopocentricFrame +# import org.orekit.orbits.KeplerianOrbit; +from org.orekit.orbits import OrbitType, KeplerianOrbit +from org.orekit.orbits import PositionAngle +# import org.orekit.propagation.Propagator; +from org.orekit.propagation.analytical import KeplerianPropagator +from org.orekit.propagation.analytical.tle import TLE, TLEPropagator +from org.orekit.propagation import SpacecraftState +# import org.orekit.propagation.conversion.NumericalPropagatorBuilder; +from org.orekit.time import AbsoluteDate, TimeScalesFactory +# import org.orekit.time.TimeScalesFactory; +from org.orekit.estimation.measurements import ObservableSatellite, PV, GroundStation, AngularRaDec +from PVMeasurementCreator import PVMeasurementCreator +from org.orekit.estimation.iod import IodLaplace +from org.orekit.utils import IERSConventions, Constants, AbsolutePVCoordinates +from org.orekit.bodies import OneAxisEllipsoid, GeodeticPoint +from orekit import JArray_double + +#from .EstimationTestUtils import EstimationTestUtils + +# Private class to calculate the errors between truth and estimated orbits at +# the central observation time. + +class Result(): + def __init__(this, truth, estOrbit): + this.errorNorm = JArray_double(2) + + this.errorNorm[0] = Vector3D.distance(truth.getPosition(), + estOrbit.getPVCoordinates().getPosition()) + + this.errorNorm[1] = Vector3D.distance(truth.getVelocity(), + estOrbit.getPVCoordinates().getVelocity()) + def getErrorNorm(this): + return this.errorNorm + + + +class IodLaplaceTest(unittest.TestCase): + + def setUp(self): + setup_orekit_curdir("resources") + self.gcrf = FramesFactory.getGCRF() + self.itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, False) + + # The ground station is set to Austin, Texas, U.S.A + body = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, + Constants.WGS84_EARTH_FLATTENING, self.itrf) + + self.observer = GroundStation(TopocentricFrame(body, GeodeticPoint(0.528253, -1.705768, 0.0), "Austin")) + + self.observer.getPrimeMeridianOffsetDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH) + self.observer.getPolarOffsetXDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH) + self.observer.getPolarOffsetYDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH) + + + def getLOSAngles(self, prop, date, observer): + # Helper function to generate a Line of Sight angles measurement for the given + # observer and date using the TLE propagator. + + satPvc = prop.getPVCoordinates(date, self.gcrf) + raDec = AngularRaDec(observer, self.gcrf, date, JArray_double([0.0, 0.0]), + JArray_double([1.0, 1.0]), + JArray_double([1.0, 1.0]), + ObservableSatellite(0)) + + angular = raDec.estimate(0, 0, [SpacecraftState( + AbsolutePVCoordinates(self.gcrf, satPvc))]).getEstimatedValue() + ra = angular[0] + dec = angular[1] + + return Vector3D(FastMath.cos(dec)*FastMath.cos(ra), + FastMath.cos(dec)*FastMath.sin(ra), FastMath.sin(dec)) + + + def estimateOrbit(self, prop, obsDate1, t2, t3): + # Helper function to generate measurements and estimate orbit for the given propagator + + # Generate 3 Line Of Sight angles measurements + los1 = self.getLOSAngles(prop, obsDate1, self.observer) + + obsDate2 = obsDate1.shiftedBy(t2) + los2 = self.getLOSAngles(prop, obsDate2, self.observer) + + obsDate3 = obsDate1.shiftedBy(t3) + los3 = self.getLOSAngles(prop, obsDate3, self.observer) + + obsPva = self.observer.getBaseFrame().getPVCoordinates(obsDate2, self.gcrf) + + # Estimate the orbit using the classical Laplace method + truth = prop.getPVCoordinates(obsDate2, self.gcrf) + estOrbit = IodLaplace(Constants.EGM96_EARTH_MU).estimate(self.gcrf, obsPva, obsDate1, los1, obsDate2, los2, obsDate3, + los3) + return Result(truth, estOrbit) + + def testLaplaceKeplerian1(self): + # Estimate the orbit of ISS (ZARYA) based on Keplerian motion + + date = AbsoluteDate(2019, 9, 29, 22, 0, 2.0, TimeScalesFactory.getUTC()) + kep = KeplerianOrbit(6798938.970424857, 0.0021115522920270016, 0.9008866630545347, + 1.8278985811406743, -2.7656136723308524, + 0.8823034512437679, PositionAngle.MEAN, self.gcrf, + date, Constants.EGM96_EARTH_MU) + + prop = KeplerianPropagator(kep) + + # With only 3 measurements, we can expect ~400 meters error in position and ~1 m/s in velocity + error = self.estimateOrbit(prop, date, 30.0, 60.0).getErrorNorm() + self.assertAlmostEquals(0.0, error[0], delta=275.0) + self.assertAlmostEquals(0.0, error[1], delta=0.8) + + + def testLaplaceKeplerian2(self): + # Estimate the orbit of Galaxy 15 based on Keplerian motion + date = AbsoluteDate(2019, 9, 29, 22, 0, 2.0, TimeScalesFactory.getUTC()) + kep = KeplerianOrbit(42165414.60406032, 0.00021743441091199163, 0.0019139259842569903, + 1.8142608912728584, 1.648821262690012, + 0.11710513241172144, PositionAngle.MEAN, self.gcrf, + date, Constants.EGM96_EARTH_MU) + + prop = KeplerianPropagator(kep) + + error = self.estimateOrbit(prop, date, 60.0, 120.0).getErrorNorm() + self.assertAlmostEquals(0.0, error[0], delta=395.0) + self.assertAlmostEquals(0.0, error[1], delta=0.03) + + + def testLaplaceTLE1(self): + # Estimate the orbit of ISS (ZARYA) based on TLE propagation + + tle1 = "1 25544U 98067A 19271.53261574 .00000256 00000-0 12497-4 0 9993" + tle2 = "2 25544 51.6447 208.7465 0007429 92.6235 253.7389 15.50110361191281" + + tleParser = TLE(tle1, tle2) + tleProp = TLEPropagator.selectExtrapolator(tleParser) + obsDate1 = tleParser.getDate() + + error = self.estimateOrbit(tleProp, obsDate1, 30.0, 60.0).getErrorNorm() + + # With only 3 measurements, an error of 5km in position and 10 m/s in velocity is acceptable + # because the Laplace method uses only two-body dynamics + self.assertAlmostEquals(0.0, error[0], delta=5000.0) + self.assertAlmostEquals(0.0, error[1], delta=10.0) + + + +if __name__ == '__main__': + suite = unittest.TestLoader().loadTestsFromTestCase(IodLaplaceTest) + ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() + sys.exit(ret) \ No newline at end of file diff --git a/python_files/test/KeplerianConverterTest.py b/python_files/test/KeplerianConverterTest.py index 885dff70c7d425017caa8061b9398ad5f9e5e09d..c8d5d7d037dfff75c5ccfa4011f9389cb68180e3 100644 --- a/python_files/test/KeplerianConverterTest.py +++ b/python_files/test/KeplerianConverterTest.py @@ -111,7 +111,7 @@ class KeplerianConverterTest(unittest.TestCase): # self.checkFit(self.orbit, 86400, 300, 1.0e-3, True, 2.65e-8, "toto"); def setUp(self): - setup_orekit_curdir() + setup_orekit_curdir("resources") self.initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.) self.orbit = EquinoctialOrbit(PVCoordinates(self.position, self.velocity), diff --git a/python_files/test/KlobucharModelTest.py b/python_files/test/KlobucharModelTest.py index 8ad4c407a0e50ac999034a383bbcacf5e9f7ac14..e9147e46dab928d4bd97f3449be27edc3d922b32 100644 --- a/python_files/test/KlobucharModelTest.py +++ b/python_files/test/KlobucharModelTest.py @@ -31,7 +31,7 @@ import orekit orekit.initVM() from orekit import JArray_double -from org.orekit.data import DataProvidersManager, ZipJarCrawler +from org.orekit.data import DataProvidersManager, ZipJarCrawler, DataContext, DirectoryCrawler from org.orekit.propagation.integration import PythonFieldAdditionalEquations from org.orekit.forces.gravity.potential import GravityFieldFactory from org.orekit.forces.gravity.potential import SHMFormatReader @@ -40,7 +40,6 @@ from java.lang import System from org.hipparchus.geometry.euclidean.threed import Vector3D -from org.hipparchus import RealFieldElement from org.hipparchus.ode.nonstiff import DormandPrince853Integrator from org.hipparchus.ode.nonstiff import DormandPrince853FieldIntegrator @@ -77,12 +76,12 @@ class KlobucharModelTest(unittest.TestCase): [.1430e+06, 0.0, -.3280e+06, .1130e+06]) # Initialize the data sources - DM = DataProvidersManager.getInstance() - datafile = File('resources.zip') + DM = DataContext.getDefault().getDataProvidersManager() + datafile = File('resources') if not datafile.exists(): print('File :', datafile.absolutePath, ' not found') - crawler = ZipJarCrawler(datafile) + crawler = DirectoryCrawler(datafile) DM.clearProviders() DM.addProvider(crawler) diff --git a/python_files/test/MeasurementCreator.py b/python_files/test/MeasurementCreator.py new file mode 100644 index 0000000000000000000000000000000000000000..fbab9fe86d5398d86a177abc15756cc7a5e98ab9 --- /dev/null +++ b/python_files/test/MeasurementCreator.py @@ -0,0 +1,47 @@ +"""/* Copyright 2002-2022 CS GROUP + * 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. + */ + + Python version translated from Java by Petrus Hyvönen, SSC 2022 + """ + +import math +import sys +import unittest + +# Python orekit specifics +import orekit +orekit.initVM() + +#from org.orekit.estimation.measurements +from org.orekit.propagation.sampling import PythonOrekitFixedStepHandler + +class MeasurementCreator(PythonOrekitFixedStepHandler): + + def init(self, s0, t, step): + self.measurements = [] + + def addMeasurement(self, measurement): + self.measurements.append(measurement) + + def getMeasurements(self): + return self.measurements + + #def handleStep(self, currentState): + # pass + + def finish(self, finalState): + pass diff --git a/python_files/test/NodeDetectorTest.py b/python_files/test/NodeDetectorTest.py index e9fb5613c75f624fee9aa8495536f433bebcfc09..e1db903b55296c531ac3917ae50b9d7bd2585255 100644 --- a/python_files/test/NodeDetectorTest.py +++ b/python_files/test/NodeDetectorTest.py @@ -33,7 +33,7 @@ import unittest from orekit.pyhelpers import setup_orekit_curdir -setup_orekit_curdir() # orekit-data.zip shall be in current dir +setup_orekit_curdir("resources") # orekit-data.zip shall be in current dir from org.orekit.propagation.events import EventsLogger from org.orekit.propagation.events import NodeDetector @@ -92,7 +92,7 @@ class NodeDetectorTest(unittest.TestCase): propagator.addEventDetector(node2) # First propagation - propagator.setEphemerisMode() + generator = propagator.getEphemerisGenerator() propagator.propagate(finalDate) assert 1998 == logger1.getLoggedEvents().size() @@ -100,7 +100,7 @@ class NodeDetectorTest(unittest.TestCase): logger1.clearLoggedEvents() logger2.clearLoggedEvents() - postpro = propagator.getGeneratedEphemeris() + postpro = generator.getGeneratedEphemeris() # Post-processing postpro.addEventDetector(node1) diff --git a/python_files/test/OrekitFixedStepHandlerMultiplexerTest.py b/python_files/test/OrekitFixedStepHandlerMultiplexerTest.py deleted file mode 100644 index a75c03d7e35e3e426ac4d3d84e80bf4b49468dcf..0000000000000000000000000000000000000000 --- a/python_files/test/OrekitFixedStepHandlerMultiplexerTest.py +++ /dev/null @@ -1,130 +0,0 @@ -# -*- coding: utf-8 -*- - -""" - -/* Copyright 2002-2019 CS Syst��mes d'Information - * Licensed to CS Syst��mes 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. - */ - -Python version translated from Java by Olivier Podevin, CS Group 2020 - - """ - -import orekit - -orekit.initVM() -from orekit.pyhelpers import setup_orekit_curdir, datetime_to_absolutedate - -setup_orekit_curdir() - -import unittest -import sys -import os -import math -from datetime import datetime - -from java.util import ArrayList, List, HashMap, Map -from java.lang import Double, StringBuffer -from java.io import BufferedReader, StringReader - -from org.hipparchus.geometry.euclidean.threed import Rotation, Vector3D -from org.hipparchus.util import FastMath - -from org.orekit.attitudes import AttitudeProvider, InertialProvider -from org.orekit.bodies import CelestialBodyFactory -from org.orekit.files.ccsds import AEMParser, AEMWriter, AEMFile, Keyword, StreamingAemWriter -from org.orekit.frames import FramesFactory -from org.orekit.orbits import CartesianOrbit, KeplerianOrbit, PositionAngle -from org.orekit.propagation import SpacecraftState -from org.orekit.propagation.analytical import KeplerianPropagator -from org.orekit.propagation.sampling import OrekitFixedStepHandler, OrekitFixedStepHandlerMultiplexer, PythonOrekitFixedStepHandler -from org.orekit.time import AbsoluteDate, TimeScalesFactory, TimeScale -from org.orekit.utils import Constants, IERSConventions, PVCoordinates - - - -class InitCheckerHandler(PythonOrekitFixedStepHandler): - - def __init__(self, expected, initialized = False): - super(InitCheckerHandler, self).__init__() - self.expected = expected - self.initialized = initialized - - def init(self, s0, t, step): - self.initialized = True - - def handleStep(self, currentState, isLast): - self.expected = 2.0 - - def isInitialized(self): - return self.initialized - - def getExpected(self): - return self.expected - -class IncrementationHandler(PythonOrekitFixedStepHandler): - - def __init__(self): - super(IncrementationHandler, self).__init__() - self.value = 0 - - def init(self, s0, t, step): - self.value = 1 - - def handleStep(self, currentState, isLast): - self.value += 1 - - def getValue(self): - return self.value - -class OrekitFixedStepHandlerMultiplexerTest(unittest.TestCase): - - def testOrekitFixedStepHandlerMultiplexer(self): - """Test the Orekit fixed StepHandler Multiplexer functionality""" - - # init - initDate = AbsoluteDate(2020, 2, 28, 16, 15, 0.0, TimeScalesFactory.getUTC()) - - initHandler = InitCheckerHandler(1.0) - incrementationHandler = IncrementationHandler() - - multiplexer = OrekitFixedStepHandlerMultiplexer() - multiplexer.add(initHandler) - multiplexer.add(incrementationHandler) - - ic = KeplerianOrbit(6378137 + 500e3, 1e-3, 0.0, 0.0, 0.0, 0.0, - PositionAngle.TRUE, - FramesFactory.getGCRF(), - initDate, - Constants.WGS84_EARTH_MU) - - propagator = KeplerianPropagator(ic) - propagator.setMasterMode(60.0, multiplexer) - - self.assertFalse(initHandler.isInitialized()) - self.assertTrue(math.isclose(1.0, initHandler.getExpected(), rel_tol = 0.0, abs_tol = Double.MIN_VALUE)) - - propagator.propagate(initDate.shiftedBy(90.0)) - - # verify - self.assertTrue(initHandler.isInitialized()) - self.assertTrue(math.isclose(2.0, initHandler.getExpected(), rel_tol = 0.0, abs_tol = Double.MIN_VALUE)) - self.assertEqual(3, incrementationHandler.getValue()) - -if __name__ == '__main__': - suite = unittest.TestLoader().loadTestsFromTestCase(OrekitFixedStepHandlerMultiplexerTest) - ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() - sys.exit(ret) diff --git a/python_files/test/OrekitStepHandlerTest.py b/python_files/test/OrekitStepHandlerTest.py index 7e6fd6c3a132540420c53933249581f6fee45098..a4ff2e53d1119def20243ed4c2097031135b1689 100644 --- a/python_files/test/OrekitStepHandlerTest.py +++ b/python_files/test/OrekitStepHandlerTest.py @@ -37,7 +37,7 @@ jcc = orekit.initVM() # import static org.junit.Assert.assertNotNull # import static org.junit.Assert.assertTrue -from org.orekit.data import DataProvidersManager, ZipJarCrawler +from org.orekit.data import DataProvidersManager, ZipJarCrawler, DataContext, DirectoryCrawler from java.io import File from java.util import Arrays @@ -55,6 +55,7 @@ from org.orekit.orbits import PositionAngle from org.orekit.propagation import Propagator from org.orekit.propagation import SpacecraftState from org.orekit.propagation.analytical import KeplerianPropagator +from org.orekit.propagation import Propagator from org.orekit.propagation.events import DateDetector from org.orekit.propagation.events.handlers import ContinueOnEvent from org.orekit.propagation.numerical import NumericalPropagator @@ -98,10 +99,13 @@ class OrekitStepHandlerTest(unittest.TestCase): def init(self, s0, t, step): # All native defined calls needs to be implemented pass - def handleStep(self, currentState, isLast): + def handleStep(self, currentState): + pass + + def finish(self, s): pass - kepler.setMasterMode(fixedStepSize, MyFixedHandler()) + Propagator.cast_(kepler).setStepHandler(fixedStepSize, MyFixedHandler()) kepler.propagate(initialDate.shiftedBy(propagationTime)) stepSizeInSeconds = 120 @@ -143,11 +147,14 @@ class OrekitStepHandlerTest(unittest.TestCase): def init(self, s0, t): pass - def handleStep(self, interpolator, isLast): + def handleStep(self, interpolator): assert (expected.popleft() == interpolator.isPreviousStateInterpolated()) assert (expected.popleft() == interpolator.isCurrentStateInterpolated()) - propagator.setMasterMode(MyHandler()) + def finish(self, s): + pass + + Propagator.cast_(propagator).setStepHandler(MyHandler()) end = date.shiftedBy(120.0) prop_end = propagator.propagate(end) self.assertEqual(end, prop_end.getDate()) @@ -156,12 +163,12 @@ class OrekitStepHandlerTest(unittest.TestCase): def setUp(self): # Utils.setDataRoot("regular-data") - DM = DataProvidersManager.getInstance() - datafile = File('resources.zip') + DM = DataContext.getDefault().getDataProvidersManager() + datafile = File('resources') if not datafile.exists(): print('File :', datafile.absolutePath, ' not found') - crawler = ZipJarCrawler(datafile) + crawler = DirectoryCrawler(datafile) DM.clearProviders() DM.addProvider(crawler) diff --git a/python_files/test/PVMeasurementCreator.py b/python_files/test/PVMeasurementCreator.py new file mode 100644 index 0000000000000000000000000000000000000000..efe9e3bfa6ff4d12731cd7d14b3318b608da68f2 --- /dev/null +++ b/python_files/test/PVMeasurementCreator.py @@ -0,0 +1,42 @@ +"""/* Copyright 2002-2022 CS GROUP + * 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. + */ + + Python version translated from Java by Petrus Hyvönen, SSC 2022 + """ + +import math +import sys +import unittest + +# Python orekit specifics +import orekit +orekit.initVM() + +from MeasurementCreator import MeasurementCreator +from org.orekit.estimation.measurements import ObservableSatellite, PV + +class PVMeasurementCreator(MeasurementCreator): + def __init__(self): + super().__init__() + self.satellite = ObservableSatellite(0) + + def handleStep(self, current_state): + p = current_state.getPVCoordinates().getPosition() + v = current_state.getPVCoordinates().getVelocity() + measurement = PV(current_state.getDate(), p, v, 1.0, 0.001, 1.0, self.satellite) + self.addMeasurement(measurement) + diff --git a/python_files/test/SmallManeuverAnalyticalModelTest.py b/python_files/test/SmallManeuverAnalyticalModelTest.py index 0fd306d934c72b0a1c7ea05024fbe6f7991d927c..605aae8181f89e568e39c36218849e90c700b243 100644 --- a/python_files/test/SmallManeuverAnalyticalModelTest.py +++ b/python_files/test/SmallManeuverAnalyticalModelTest.py @@ -27,7 +27,7 @@ import orekit orekit.initVM() from orekit.pyhelpers import setup_orekit_curdir -setup_orekit_curdir() +setup_orekit_curdir("resources") from orekit import JArray_double, JArray @@ -44,7 +44,7 @@ from org.orekit.frames import LOFType from org.orekit.orbits import CircularOrbit from org.orekit.orbits import KeplerianOrbit from org.orekit.orbits import PositionAngle -from org.orekit.propagation import SpacecraftState +from org.orekit.propagation import SpacecraftState, Propagator from org.orekit.propagation.numerical import NumericalPropagator from org.orekit.time import AbsoluteDate from org.orekit.time import DateComponents @@ -213,9 +213,9 @@ class SmallManeuverAnalyticalModelTest(unittest.TestCase): maneuver = ConstantThrustManeuver(t0, dt, f, isp, dV.normalize()) propagator.addForceModel(maneuver) - propagator.setEphemerisMode() + generator = propagator.getEphemerisGenerator() propagator.propagate(t0.shiftedBy(5 * orbit.getKeplerianPeriod())) - return propagator.getGeneratedEphemeris() + return generator.getGeneratedEphemeris() if __name__ == '__main__': diff --git a/python_files/test/SpinStabilizedTest.py b/python_files/test/SpinStabilizedTest.py index 60bcb40f3e301214611f45a5c5c6b21e42dee9de..4f479417a54c58fb7a34cf0eb3ccb83974692044 100644 --- a/python_files/test/SpinStabilizedTest.py +++ b/python_files/test/SpinStabilizedTest.py @@ -30,7 +30,7 @@ import unittest import orekit orekit.initVM() -from org.orekit.data import DataProvidersManager, ZipJarCrawler +from org.orekit.data import DataProvidersManager, ZipJarCrawler, DataContext, DirectoryCrawler from java.io import File from org.hipparchus.geometry.euclidean.threed import Rotation @@ -60,12 +60,12 @@ from org.orekit.attitudes import CelestialBodyPointed, SpinStabilized, InertialP class SpinStabilizedTest(unittest.TestCase): def setUp(self): - DM = DataProvidersManager.getInstance() - datafile = File('resources.zip') + DM = DataContext.getDefault().getDataProvidersManager() + datafile = File('resources') if not datafile.exists(): print('File :', datafile.absolutePath, ' not found') - crawler = ZipJarCrawler(datafile) + crawler = DirectoryCrawler(datafile) DM.clearProviders() DM.addProvider(crawler) diff --git a/python_files/test/TLEConverterTest.py b/python_files/test/TLEConverterTest.py deleted file mode 100644 index 06e5baaeda079581d9a28cc98085c164dcd110d2..0000000000000000000000000000000000000000 --- a/python_files/test/TLEConverterTest.py +++ /dev/null @@ -1,143 +0,0 @@ -# -*- coding: utf-8 -*- - -""" - -/* Copyright 2002-2013 CS Syst��mes d'Information - * Licensed to CS Syst��mes 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. - */ - -Java version translated to Python by Petrus Hyvönen, SSC 2017 - - """ - -import orekit - -orekit.initVM() -from orekit.pyhelpers import setup_orekit_curdir - -from org.orekit.orbits import PositionAngle -from org.orekit.propagation.conversion import FiniteDifferencePropagatorConverter -from org.orekit.propagation.conversion import TLEPropagatorBuilder -from org.orekit.propagation.analytical.tle import TLE -from org.orekit.propagation.analytical.tle import TLEPropagator -from org.orekit.data import DataProvidersManager, ZipJarCrawler -from java.util import Arrays -from java.io import File - -import unittest -import sys - -class TLEConverterTest(unittest.TestCase): - - def checkFit(self, tle, - duration, - stepSize, - threshold, - positionOnly, - withBStar, - expectedRMS): - - p = TLEPropagator.selectExtrapolator(tle) - - sample = [] - dt = 0.0 - while dt < duration: - sample.append(p.propagate(tle.getDate().shiftedBy(dt))) - dt += stepSize - - builder = TLEPropagatorBuilder(tle, PositionAngle.TRUE, 1.0) - - drivers = builder.getPropagationParametersDrivers().getDrivers() - - # there should *not *be any drivers for central attraction coefficient (see issue # 313) - self.assertEqual(1, drivers.size()) - self.assertEqual("BSTAR", drivers.get(0).getName()) - - fitter = FiniteDifferencePropagatorConverter(builder, threshold, 1000) - sample = Arrays.asList(sample) - - if withBStar: - fitter.convert(sample, positionOnly, TLEPropagatorBuilder.B_STAR) - else: - fitter.convert(sample, positionOnly, []) - - prop = TLEPropagator.cast_(fitter.getAdaptedPropagator()) - fitted = prop.getTLE() - - tolerance = max(threshold, 0.001 * expectedRMS) - - self.assertAlmostEqual(expectedRMS, fitter.getRMS(), delta=tolerance) - - self.assertEqual(tle.getSatelliteNumber(), fitted.getSatelliteNumber()) - self.assertEqual(tle.getClassification(), fitted.getClassification()) - self.assertEqual(tle.getLaunchYear(), fitted.getLaunchYear()) - self.assertEqual(tle.getLaunchNumber(), fitted.getLaunchNumber()) - self.assertEqual(tle.getLaunchPiece(), fitted.getLaunchPiece()) - self.assertEqual(tle.getElementNumber(), fitted.getElementNumber()) - self.assertEqual(tle.getRevolutionNumberAtEpoch(), fitted.getRevolutionNumberAtEpoch()) - - eps = 1.0e-5 - self.assertAlmostEqual(tle.getMeanMotion(), fitted.getMeanMotion(), delta=eps * tle.getMeanMotion()) - self.assertAlmostEqual(tle.getE(), fitted.getE(), delta=eps * tle.getE()) - self.assertAlmostEqual(tle.getI(), fitted.getI(), delta=eps * tle.getI()) - self.assertAlmostEqual(tle.getPerigeeArgument(), fitted.getPerigeeArgument(), - delta=eps * tle.getPerigeeArgument()) - self.assertAlmostEqual(tle.getRaan(), fitted.getRaan(), delta=eps * tle.getRaan()) - self.assertAlmostEqual(tle.getMeanAnomaly(), fitted.getMeanAnomaly(), delta=eps * tle.getMeanAnomaly()) - - if withBStar: - self.assertAlmostEqual(tle.getBStar(), fitted.getBStar(), delta=eps * tle.getBStar()) - - def testConversionGeoPositionVelocity(self): - self.checkFit(self.geoTLE, 86400, 300, 1.0e-3, False, False, 9.350e-8) - - def testConversionGeoPositionOnly(self): - self.checkFit(self.geoTLE, 86400, 300, 1.0e-3, True, False, 1.328e-7) - - def testConversionLeoPositionVelocityWithoutBStar(self): - self.checkFit(self.leoTLE, 86400, 300, 1.0e-3, False, False, 10.77) - - def testConversionLeoPositionOnlyWithoutBStar(self): - self.checkFit(self.leoTLE, 86400, 300, 1.0e-3, True, False, 15.23) - - def testConversionLeoPositionVelocityWithBStar(self): - self.checkFit(self.leoTLE, 86400, 300, 1.0e-3, False, True, 2.646e-8) - - def testConversionLeoPositionOnlyWithBStar(self): - self.checkFit(self.leoTLE, 86400, 300, 1.0e-3, True, True, 4.102e-8) - - def setUp(self): - #setup_orekit_curdir() - - DM = DataProvidersManager.getInstance() - datafile = File('resources.zip') - if not datafile.exists(): - print('File :', datafile.absolutePath, ' not found') - - crawler = ZipJarCrawler(datafile) - DM.clearProviders() - DM.addProvider(crawler) - - self.geoTLE = TLE("1 27508U 02040A 12021.25695307 -.00000113 00000-0 10000-3 0 7326", - "2 27508 0.0571 356.7800 0005033 344.4621 218.7816 1.00271798 34501") - self.leoTLE = TLE("1 31135U 07013A 11003.00000000 .00000816 00000+0 47577-4 0 11", - "2 31135 2.4656 183.9084 0021119 236.4164 60.4567 15.10546832 15") - - -if __name__ == '__main__': - suite = unittest.TestLoader().loadTestsFromTestCase(TLEConverterTest) - ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() - sys.exit(ret) diff --git a/python_files/test/TransformTest.py b/python_files/test/TransformTest.py index c973605eb1d9cfb1feab8effa4fbaee3d927a259..a07ebb60a5f5ea9a4e0806e5138b628fdcd92e2a 100644 --- a/python_files/test/TransformTest.py +++ b/python_files/test/TransformTest.py @@ -28,7 +28,7 @@ import orekit orekit.initVM() from orekit.pyhelpers import setup_orekit_curdir, datetime_to_absolutedate -setup_orekit_curdir() +setup_orekit_curdir("resources") from org.orekit.time import AbsoluteDate import unittest diff --git a/python_files/test/__init__.py b/python_files/test/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/python_files/test/orekit-data.zip b/python_files/test/orekit-data.zip deleted file mode 100644 index 564a2ad7eb4d223c5a532737d1960e5a5a2bee97..0000000000000000000000000000000000000000 Binary files a/python_files/test/orekit-data.zip and /dev/null differ diff --git a/python_files/test/resources.zip b/python_files/test/resources.zip deleted file mode 100644 index f48326c55e155ed95e5c7254a76d24ed25d67466..0000000000000000000000000000000000000000 Binary files a/python_files/test/resources.zip and /dev/null differ