diff --git a/orekit-java-wrappers/src/main/java/org/orekit/python/PythonAbstractDetector.java b/java_additions/src/main/java/org/orekit/python/PythonAbstractDetector.java similarity index 100% rename from orekit-java-wrappers/src/main/java/org/orekit/python/PythonAbstractDetector.java rename to java_additions/src/main/java/org/orekit/python/PythonAbstractDetector.java diff --git a/orekit-java-wrappers/src/main/java/org/orekit/python/PythonEventDetector.java b/java_additions/src/main/java/org/orekit/python/PythonEventDetector.java similarity index 100% rename from orekit-java-wrappers/src/main/java/org/orekit/python/PythonEventDetector.java rename to java_additions/src/main/java/org/orekit/python/PythonEventDetector.java diff --git a/orekit-java-wrappers/src/main/java/org/orekit/python/PythonEventHandler.java b/java_additions/src/main/java/org/orekit/python/PythonEventHandler.java similarity index 100% rename from orekit-java-wrappers/src/main/java/org/orekit/python/PythonEventHandler.java rename to java_additions/src/main/java/org/orekit/python/PythonEventHandler.java diff --git a/orekit-java-wrappers/src/main/java/org/orekit/python/PythonOrekitFixedStepHandler.java b/java_additions/src/main/java/org/orekit/python/PythonOrekitFixedStepHandler.java similarity index 100% rename from orekit-java-wrappers/src/main/java/org/orekit/python/PythonOrekitFixedStepHandler.java rename to java_additions/src/main/java/org/orekit/python/PythonOrekitFixedStepHandler.java diff --git a/orekit-java-wrappers/src/main/java/org/orekit/python/PythonUnivariateFunction.java b/java_additions/src/main/java/org/orekit/python/PythonUnivariateFunction.java similarity index 100% rename from orekit-java-wrappers/src/main/java/org/orekit/python/PythonUnivariateFunction.java rename to java_additions/src/main/java/org/orekit/python/PythonUnivariateFunction.java diff --git a/orekit-java-wrappers/src/main/java/org/orekit/python/package-info.java b/java_additions/src/main/java/org/orekit/python/package-info.java similarity index 100% rename from orekit-java-wrappers/src/main/java/org/orekit/python/package-info.java rename to java_additions/src/main/java/org/orekit/python/package-info.java diff --git a/orekit-conda-recipe/meta.yaml b/orekit-conda-recipe/meta.yaml index 4c1a48f0be0afffe9015b154b336fd1f0e7005da..85fc3590126ded1d628c1d58642019f561873b60 100644 --- a/orekit-conda-recipe/meta.yaml +++ b/orekit-conda-recipe/meta.yaml @@ -15,7 +15,7 @@ source: build: skip: true # [win32 or linux32] - number: 4 + number: 3 rpaths: - lib/ diff --git a/python_files/README.md b/python_files/README.md new file mode 100644 index 0000000000000000000000000000000000000000..25f05ddb51bc9b34af2cd76482e6c2489008bf0f --- /dev/null +++ b/python_files/README.md @@ -0,0 +1,2 @@ +# orekit_python_artifacts +Primary jar's for python orekit diff --git a/python_files/pyhelpers.py b/python_files/pyhelpers.py new file mode 100644 index 0000000000000000000000000000000000000000..3df7e8ee2582f81417d74b37f89adff6fa1377e3 --- /dev/null +++ b/python_files/pyhelpers.py @@ -0,0 +1,123 @@ +# encoding: utf-8 + +# Copyright 2014 SSC +# Licensed 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. + +""" This document contains classes that are useful for using the orekit +library in Python. """ + +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function + +# Set up the orekit namespace +import orekit + +from java.io import File + +from org.orekit.data import DataProvidersManager, ZipJarCrawler +from org.orekit.time import TimeScalesFactory, AbsoluteDate +from org.orekit.utils import ElevationMask +from orekit import JArray + +import math +from datetime import datetime + + +def setup_orekit_curdir(filename='orekit-data.zip'): + """Setup the java engine with orekit. + + This function loads the orekit-data.zip from the current directory + and sets up the Orekit DataProviders to access it. + + The JVM needs to be initiated prior to calling this function: + + orekit.initVM() + + Args: + filename (str): Name of zip with orekit data. Default filename is 'orekit-data.zip' + + """ + + DM = DataProvidersManager.getInstance() + datafile = File(filename) + if not datafile.exists(): + print('File :', datafile.absolutePath, ' not found') + + crawler = ZipJarCrawler(datafile) + DM.clearProviders() + DM.addProvider(crawler) + + +def absolutedate_to_datetime(orekit_absolutedate): + """ Converts between orekit.AbsoluteDate objects + and python datetime objects (utc)""" + + utc = TimeScalesFactory.getUTC() + or_comp = orekit_absolutedate.getComponents(utc) + or_date = or_comp.getDate() + or_time = or_comp.getTime() + seconds = or_time.getSecond() + return datetime(or_date.getYear(), + or_date.getMonth(), + or_date.getDay(), + or_time.getHour(), + or_time.getMinute(), + int(math.floor(seconds)), + int(1000000.0 * (seconds - math.floor(seconds)))) + +def datetime_to_absolutedate(dt_date): + """ Converts between orekit.AbsoluteDate objects + and python datetime objects (utc) + + Args: + dt_date (datetime): python datetime object to convert + + Returns: + AbsoluteDate: time in orekit format""" + + utc = TimeScalesFactory.getUTC() + return AbsoluteDate(dt_date.year, + dt_date.month, + dt_date.day, + dt_date.hour, + dt_date.minute, + dt_date.second + dt_date.microsecond / 1000000., + utc) + + +def to_elevationmask(az, el): + ''' Converts an array of azimuths and elevations to a + orekit ElevationMask object. All unts in degrees. + + mask = to_elevationmask([0, 90, 180, 270], [5,10,8,5]) + + ''' + + mask = JArray('object')(len(az)) + + for i in range(len(az)): + mask[i] = JArray('double')([math.radians(az[i]), + math.radians(el[i])]) + + return ElevationMask(mask) + +def JArray_double2D(x, y): + '''Returns an JCC wrapped 2D double array''' + + arr = JArray('object')(x) + + for i in range(x): + arr[i] = JArray('double')(y) + + return arr \ No newline at end of file diff --git a/python_files/test/AltitudeDetectorTest.py b/python_files/test/AltitudeDetectorTest.py new file mode 100644 index 0000000000000000000000000000000000000000..28177b5a60172fb9b92ca6101d1ff3ccc8570b44 --- /dev/null +++ b/python_files/test/AltitudeDetectorTest.py @@ -0,0 +1,77 @@ +# -*- 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. + */ + +Python version translated from Java by Petrus Hyvönen, SSC 2014 + + """ + +#Python orekit specifics +import orekit +orekit.initVM() +from orekit.pyhelpers import setup_orekit_curdir +setup_orekit_curdir() # orekit-data.zip shall be in current dir +#from math import abs + +from org.hipparchus.util import FastMath +from org.orekit.bodies import CelestialBodyFactory +from org.orekit.bodies import OneAxisEllipsoid +from org.orekit.frames import FramesFactory +from org.orekit.orbits import KeplerianOrbit +from org.orekit.orbits import PositionAngle +from org.orekit.propagation import SpacecraftState +from org.orekit.propagation.analytical import KeplerianPropagator +from org.orekit.propagation.events.handlers import StopOnEvent +from org.orekit.time import AbsoluteDate +from org.orekit.time import TimeScalesFactory +from org.orekit.propagation.events import AltitudeDetector + +EME2000 = FramesFactory.getEME2000() +initialDate = AbsoluteDate(2009,1,1,TimeScalesFactory.getUTC()) +a = 8000000.0 +e = 0.1 +earthRadius = 6378137.0 +earthF = 1.0 / 298.257223563 +apogee = a*(1+e) +alt = apogee - earthRadius - 500 + +#// initial state is at apogee +initialOrbit = KeplerianOrbit(a,e,0.0,0.0,0.0,FastMath.PI,PositionAngle.MEAN,EME2000, + initialDate,CelestialBodyFactory.getEarth().getGM()) +initialState = SpacecraftState(initialOrbit) +kepPropagator = KeplerianPropagator(initialOrbit) +altDetector = AltitudeDetector(alt, + OneAxisEllipsoid(earthRadius, earthF, EME2000)).withHandler(StopOnEvent().of_(AltitudeDetector)) + +# altitudeDetector should stop propagation upon reaching required altitude +kepPropagator.addEventDetector(altDetector) + +#// propagation to the future +finalState = kepPropagator.propagate(initialDate.shiftedBy(1000.0)) +assert abs(finalState.getPVCoordinates().getPosition().getNorm()-earthRadius -alt)<1e-5 +assert abs(44.079 - finalState.getDate().durationFrom(initialDate))< 1.0e-3 + +#// propagation to the past +kepPropagator.resetInitialState(initialState) +finalState = kepPropagator.propagate(initialDate.shiftedBy(-1000.0)) +assert abs(finalState.getPVCoordinates().getPosition().getNorm()-earthRadius - alt)< 1e-5 +assert abs(-44.079 - finalState.getDate().durationFrom(initialDate))< 1.0e-3 + +print("AltitudeDetectorTest successfully run") diff --git a/python_files/test/BackAndForthDetectorTest.py b/python_files/test/BackAndForthDetectorTest.py new file mode 100644 index 0000000000000000000000000000000000000000..108e2a4644a3d645947efcbd8cd8d3863e01d543 --- /dev/null +++ b/python_files/test/BackAndForthDetectorTest.py @@ -0,0 +1,115 @@ +# -*- 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. + */ + +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() + +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 +from org.orekit.propagation.analytical import KeplerianPropagator +from org.orekit.propagation.events.handlers import EventHandler +from org.orekit.python import 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 +from org.orekit.propagation.events import ElevationDetector +import unittest +import sys +import math + + +class Visibility(PythonEventHandler): # implements EventHandler<ElevationDetector> { + + def __init__(self): + self._visiNb = 0 + super(Visibility, self).__init__() + + def getVisiNb(self): + return self._visiNb + + def eventOccurred(self, s, ed, increasing): + self._visiNb += 1 + return EventHandler.Action.CONTINUE +# + def resetState(self, detector, oldState): + return oldState + +class BackAndForthDetectorTest(unittest.TestCase): + + def testBackAndForth(self): + utc = TimeScalesFactory.getUTC() + date0 = AbsoluteDate(2006, 12, 27, 12, 0, 0.0, utc) + date1 = AbsoluteDate(2006, 12, 27, 22, 50, 0.0, utc) + date2 = AbsoluteDate(2006, 12, 27, 22, 58, 0.0, utc) + + # Orbit + a = 7274000.0 + e = 0.00127 + i = math.radians(90.) + w = math.radians(0.) + raan = math.radians(12.5) + lM = math.radians(60.) + iniOrb = KeplerianOrbit(a, e, i, w, raan, lM, + PositionAngle.MEAN, + FramesFactory.getEME2000(), date0, + Constants.WGS84_EARTH_MU) + + # Propagator + propagator = KeplerianPropagator(iniOrb) + + # Station + stationPosition = GeodeticPoint(math.radians(0.), math.radians(100.), 110.) + earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, + Constants.WGS84_EARTH_FLATTENING, + FramesFactory.getITRF(IERSConventions.IERS_2010, True)) + + stationFrame = TopocentricFrame(earth, stationPosition, "") + + # Detector + visi = Visibility() #.of_(ElevationDetector); + det = ElevationDetector(stationFrame).withConstantElevation(math.radians(10.0)).withHandler(visi) + propagator.addEventDetector(det) + + # Forward propagation (AOS + LOS) + propagator.propagate(date1) + propagator.propagate(date2) + # Backward propagation (AOS + LOS) + propagator.propagate(date1) + propagator.propagate(date0) + + self.assertEquals(4, visi.getVisiNb()) + +if __name__ == '__main__': + suite = unittest.TestLoader().loadTestsFromTestCase(BackAndForthDetectorTest) + ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() + sys.exit(ret) + diff --git a/python_files/test/EventDetectorTest.py b/python_files/test/EventDetectorTest.py new file mode 100644 index 0000000000000000000000000000000000000000..fc9e199eaebcd868db12d993caf8195262d051ba --- /dev/null +++ b/python_files/test/EventDetectorTest.py @@ -0,0 +1,129 @@ +# -*- coding: utf-8 -*- + +""" + +/* Copyright 2014 SSC + * 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. + */ + + + """ + +import orekit +orekit.initVM() +#orekit.initVM(vmargs='-Xcheck:jni,-verbose:jni,-verbose:class,-XX:+UnlockDiagnosticVMOptions') + +from org.orekit.frames import FramesFactory, TopocentricFrame +from org.orekit.bodies import OneAxisEllipsoid, GeodeticPoint +from org.orekit.time import AbsoluteDate, TimeScalesFactory +from org.orekit.orbits import KeplerianOrbit +from org.orekit.utils import Constants +from org.orekit.propagation.analytical import KeplerianPropagator +from org.orekit.utils import PVCoordinates, IERSConventions +from org.orekit.propagation.events.handlers import EventHandler +from org.hipparchus.geometry.euclidean.threed import Vector3D +from org.orekit.python import PythonEventHandler, PythonAbstractDetector, PythonEventDetector + +from math import radians +import math +import unittest +import sys + +from orekit.pyhelpers import setup_orekit_curdir +setup_orekit_curdir() + +class MyElevationDetector(PythonEventDetector): + passes = 0 + + def __init__(self, elevation, topo): + self.elevation = elevation + self.topo = topo + super(MyElevationDetector, self).__init__() + + def init(self, s, T): + pass + + def getThreshold(self): + return float(PythonAbstractDetector.DEFAULT_THRESHOLD) + + def getMaxCheckInterval(self): + return float(PythonAbstractDetector.DEFAULT_MAXCHECK) + + def getMaxIterationCount(self): + return PythonAbstractDetector.DEFAULT_MAX_ITER + + def g(self, s): + tmp = self.topo.getElevation(s.getPVCoordinates().getPosition(), s.getFrame(), s.getDate())-self.elevation + return tmp + + def eventOccurred(self, s, increasing): + if increasing: + self.passes = self.passes + 1 + + return EventHandler.Action.CONTINUE + + def resetState(self, oldState): + return oldState + + def getElevation(self): + return self.elevation + + def getTopocentricFrame(self): + return self.topo + + +class EventDetectorTest(unittest.TestCase): + + def testOwnElevationDetector(self): + + initialDate = AbsoluteDate(2014, 1, 1, 23, 30, 00.000, TimeScalesFactory.getUTC()) + inertialFrame = FramesFactory.getEME2000() # inertial frame for orbit definition + position = Vector3D(-6142438.668, 3492467.560, -25767.25680) + velocity = Vector3D(505.8479685, 942.7809215, 7435.922231) + pvCoordinates = PVCoordinates(position, velocity) + initialOrbit = KeplerianOrbit(pvCoordinates, + inertialFrame, + initialDate, + Constants.WGS84_EARTH_MU) + + kepler = KeplerianPropagator(initialOrbit) + + ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, True) + earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, + Constants.WGS84_EARTH_FLATTENING, + ITRF) + + # Station + longitude = radians(45.0) + latitude = radians(25.0) + altitude = 0 + station1 = GeodeticPoint(latitude, longitude, float (altitude)) + sta1Frame = TopocentricFrame(earth, station1, "station 1") + + elevation = math.radians(5.0) + + detector = MyElevationDetector(elevation, sta1Frame) + kepler.addEventDetector(detector) + + finalState = kepler.propagate(initialDate.shiftedBy(60*60*24.0*15)) + + print(detector.passes) + self.assertEquals(52, detector.passes) + +if __name__ == '__main__': + suite = unittest.TestLoader().loadTestsFromTestCase(EventDetectorTest) + ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() + sys.exit(ret) diff --git a/python_files/test/EventHandlerTest.py b/python_files/test/EventHandlerTest.py new file mode 100644 index 0000000000000000000000000000000000000000..8681a732cec762548782aef6cd732ecd3d3efaa8 --- /dev/null +++ b/python_files/test/EventHandlerTest.py @@ -0,0 +1,118 @@ +# -*- 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. + */ + +Python version translated from Java by Petrus Hyvönen, SSC 2014 + + """ + +import orekit +orekit.initVM() + +from org.orekit.frames import FramesFactory, TopocentricFrame +from org.orekit.bodies import OneAxisEllipsoid, GeodeticPoint +from org.orekit.time import AbsoluteDate, TimeScalesFactory +from org.orekit.orbits import KeplerianOrbit +from org.orekit.utils import Constants +from org.orekit.propagation.analytical import KeplerianPropagator +from org.orekit.utils import PVCoordinates, IERSConventions +from org.orekit.propagation.events import ElevationDetector +from org.orekit.propagation.events.handlers import EventHandler +from org.hipparchus.geometry.euclidean.threed import Vector3D +from org.orekit.python import PythonEventHandler +from org.orekit.propagation.events import EventsLogger + + +from math import radians +import math +import unittest +import sys + +#%% Setup Orekit +from orekit.pyhelpers import setup_orekit_curdir +setup_orekit_curdir() + +#%% +class EventHandlerTest(unittest.TestCase): + + def testOwnContinueOnEvent(self): + initialDate = AbsoluteDate(2014, 1, 1, 23, 30, 00.000, TimeScalesFactory.getUTC()) + inertialFrame = FramesFactory.getEME2000() # inertial frame for orbit definition + position = Vector3D(-6142438.668, 3492467.560, -25767.25680) + velocity = Vector3D(505.8479685, 942.7809215, 7435.922231) + pvCoordinates = PVCoordinates(position, velocity) + initialOrbit = KeplerianOrbit(pvCoordinates, + inertialFrame, + initialDate, + Constants.WGS84_EARTH_MU) + + # Propagator : consider a simple keplerian motion (could be more elaborate) + kepler = KeplerianPropagator(initialOrbit) + + #Earth and frame + ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, True) + earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, + Constants.WGS84_EARTH_FLATTENING, + ITRF) + + # Station + longitude = radians(45.0) + latitude = radians(25.0) + altitude = 0 + station1 = GeodeticPoint(latitude, longitude, float (altitude)) + sta1Frame = TopocentricFrame(earth, station1, "station 1") + + elevation = math.radians(5.0) + #%% + class myContinueOnEvent(PythonEventHandler): + + def eventOccurred(self, s, T, increasing): + return EventHandler.Action.CONTINUE + + def resetState(self, detector, oldState): + return oldState + + #%% detectors + detector = ElevationDetector(sta1Frame).withConstantElevation(elevation) + detector = detector.withHandler(myContinueOnEvent().of_(ElevationDetector)) + + logger = EventsLogger() + kepler.addEventDetector(logger.monitorDetector(detector)) + + #%%Propagate from the initial date to the first raising or for the fixed duration + finalState = kepler.propagate(initialDate.shiftedBy(60*60*24.0*15)) + + + taken_passes = 0 + + mylog = logger.getLoggedEvents() + for ev in mylog: + #print 'Date: ',ev.getState().getDate(), ' Start pass: ',ev.isIncreasing() + if ev.isIncreasing(): + taken_passes = taken_passes + 1 + + #print 'Taken passes:',taken_passes + self.assertEquals(52, taken_passes) + + +if __name__ == '__main__': + suite = unittest.TestLoader().loadTestsFromTestCase(EventHandlerTest) + ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() + sys.exit(ret) diff --git a/python_files/test/GroundFieldOfViewDetectorTest.py b/python_files/test/GroundFieldOfViewDetectorTest.py new file mode 100644 index 0000000000000000000000000000000000000000..c8967bccb169119725b68e1a6232ba9275637249 --- /dev/null +++ b/python_files/test/GroundFieldOfViewDetectorTest.py @@ -0,0 +1,111 @@ +# -*- coding: utf-8 -*- + +""" + +/* Copyright 2002-2016 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 2016 + + """ + +import orekit + +orekit.initVM() + +from org.orekit.frames import FramesFactory, TopocentricFrame +from org.orekit.bodies import OneAxisEllipsoid, GeodeticPoint +from org.orekit.time import AbsoluteDate +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.hipparchus.geometry.euclidean.threed import Vector3D + +from math import radians +import math +import unittest +import sys + +from orekit.pyhelpers import setup_orekit_curdir + +setup_orekit_curdir() + + +class GroundFieldOfViewDetectorTest(unittest.TestCase): + def testGroundFieldOfViewDetector(self): + date = AbsoluteDate.J2000_EPOCH # arbitrary date + endDate = date.shiftedBy(Constants.JULIAN_DAY) + eci = FramesFactory.getGCRF() + ecef = FramesFactory.getITRF(IERSConventions.IERS_2010, True) + earth = OneAxisEllipsoid( + Constants.WGS84_EARTH_EQUATORIAL_RADIUS, + Constants.WGS84_EARTH_FLATTENING, + ecef) + + gp = GeodeticPoint(radians(39), radians(77), 0.0) + topo = TopocentricFrame(earth, gp, "topo") + + # iss like orbit + orbit = KeplerianOrbit( + 6378137.0 + 400e3, 0.0, radians(51.65), 0.0, 0.0, 0.0, + PositionAngle.TRUE, eci, date, Constants.EGM96_EARTH_MU) + + prop = KeplerianPropagator(orbit) + + # compute expected result + elevationDetector = ElevationDetector(topo).withConstantElevation(math.pi / 6.0).withMaxCheck(5.0) + logger = EventsLogger() + prop.addEventDetector(logger.monitorDetector(elevationDetector)) + prop.propagate(endDate) + expected = logger.getLoggedEvents() + + # action + # construct similar FoV based detector + # 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) + + # 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()) + logger = EventsLogger() + + prop = KeplerianPropagator(orbit) + prop.addEventDetector(logger.monitorDetector(fovDetector)) + prop.propagate(endDate) + actual = logger.getLoggedEvents() + + # verify + self.assertEquals(2, expected.size()) + self.assertEquals(2, actual.size()) + + for i in range(0, 1): + expectedDate = expected.get(i).getState().getDate() + actualDate = actual.get(i).getState().getDate() + # same event times to within 1s. + self.assertAlmostEqual(expectedDate.durationFrom(actualDate), 0.0, delta=1.0) + + +if __name__ == '__main__': + suite = unittest.TestLoader().loadTestsFromTestCase(GroundFieldOfViewDetectorTest) + ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() + sys.exit(ret) diff --git a/python_files/test/ImpulseManeuverTest.py b/python_files/test/ImpulseManeuverTest.py new file mode 100644 index 0000000000000000000000000000000000000000..4ea62286afd91d09f420a23c1b176dfe4423b417 --- /dev/null +++ b/python_files/test/ImpulseManeuverTest.py @@ -0,0 +1,93 @@ +# -*- 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. + */ + +Python version translated from Java by Petrus Hyvönen, SSC 2014 + + """ + +#Python orekit specifics +import orekit +orekit.initVM() +from orekit.pyhelpers import setup_orekit_curdir +setup_orekit_curdir() + +from org.hipparchus.geometry.euclidean.threed import Vector3D +from org.hipparchus.util import FastMath +from org.orekit.attitudes import LofOffset +from org.orekit.frames import FramesFactory +from org.orekit.frames import LOFType +from org.orekit.orbits import KeplerianOrbit +from org.orekit.orbits import PositionAngle +from org.orekit.propagation.analytical import KeplerianPropagator +from org.orekit.propagation.events import NodeDetector +from org.orekit.time import AbsoluteDate +from org.orekit.time import DateComponents +from org.orekit.time import TimeComponents +from org.orekit.time import TimeScalesFactory +from org.orekit.forces.maneuvers import ImpulseManeuver +import unittest +import sys +import math + +class ImpulseManeuverTest(unittest.TestCase): + + def testInclinationManeuver(self): + initialOrbit = KeplerianOrbit(24532000.0, + 0.72, + 0.3, + FastMath.PI, + 0.4, + 2.0, + PositionAngle.MEAN, + FramesFactory.getEME2000(), + AbsoluteDate(DateComponents(2008, 6, 23), + TimeComponents(14, 18, 37.0), + TimeScalesFactory.getUTC()), + 3.986004415e14) + + a = initialOrbit.getA() + e = initialOrbit.getE() + i = initialOrbit.getI() + mu = initialOrbit.getMu() + vApo = math.sqrt(mu * (1 - e) / (a * (1 + e))) + dv = 0.99 * math.tan(i) * vApo + + propagator = KeplerianPropagator(initialOrbit, + LofOffset(initialOrbit.getFrame(), LOFType.VVLH)) + + det = ImpulseManeuver(NodeDetector(initialOrbit, + FramesFactory.getEME2000() ), + Vector3D(dv, Vector3D.PLUS_J), 400.0) + det = det.of_(NodeDetector) + + propagator.addEventDetector(det) + + propagated = propagator.propagate(initialOrbit.getDate().shiftedBy(8000.0)) + + self.assertAlmostEqual(0.0028257, propagated.getI(), delta=1.0e-6) + +if __name__ == '__main__': + #unittest.main() + + suite = unittest.TestLoader().loadTestsFromTestCase(ImpulseManeuverTest) + ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() + sys.exit(ret) + diff --git a/python_files/test/IodGibbsTest.py b/python_files/test/IodGibbsTest.py new file mode 100644 index 0000000000000000000000000000000000000000..ba418ce5bc1ef9fd156900e0592a2451dcadedea --- /dev/null +++ b/python_files/test/IodGibbsTest.py @@ -0,0 +1,144 @@ +# -*- coding: utf-8 -*- + +""" + +/** + * + * Source: http://ccar.colorado.edu/asen5050/projects/projects_2012/kemble/gibbs_derivation.htm + * + * @author Joris Olympio + * @since 7.1 + * + */ + +Java version translated to Python by Petrus Hyvönen, SSC 2017 + + """ + +import orekit + +orekit.initVM() +from orekit.pyhelpers import setup_orekit_curdir + +import java.util.List; + +import org.hipparchus.geometry.euclidean.threed.Vector3D; +import org.junit.Assert; +import org.junit.Test; +import org.orekit.errors.OrekitException; +import org.orekit.estimation.Context; +import org.orekit.estimation.EstimationTestUtils; +import org.orekit.estimation.measurements.ObservedMeasurement; +import org.orekit.estimation.measurements.PV; +import org.orekit.estimation.measurements.PVMeasurementCreator; +import org.orekit.frames.Frame; +import org.orekit.frames.FramesFactory; +import org.orekit.orbits.KeplerianOrbit; +import org.orekit.orbits.OrbitType; +import org.orekit.orbits.PositionAngle; +import org.orekit.propagation.Propagator; +import org.orekit.propagation.conversion.NumericalPropagatorBuilder; +import org.orekit.time.AbsoluteDate; +import org.orekit.time.TimeScalesFactory; + + + +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 IodGibbsTest(unittest.TestCase): + def testGibbs1(self): + context = EstimationTestUtils.eccentricContext("regular-data:potential:tides"); + final + double + mu = context.initialOrbit.getMu(); + final + Frame + frame = context.initialOrbit.getFrame(); + + final + NumericalPropagatorBuilder + propagatorBuilder = \ + context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, + 1.0e-6, 60.0, 0.001); + + // create + perfect + range + measurements + final + Propagator + propagator = EstimationTestUtils.createPropagator(context.initialOrbit, + propagatorBuilder); + + final + List < ObservedMeasurement <? >> measurements = + EstimationTestUtils.createMeasurements(propagator, + new + PVMeasurementCreator(), + 0.0, 1.0, 60.0); + + final + Vector3D + position1 = new + Vector3D(measurements.get(0).getObservedValue()[0], + measurements.get(0).getObservedValue()[1], + measurements.get(0).getObservedValue()[2]); + final + PV + pv1 = new + PV(measurements.get(0).getDate(), position1, Vector3D.ZERO, 0., 0., 1.); + + final + Vector3D + position2 = new + Vector3D(measurements.get(1).getObservedValue()[0], + measurements.get(1).getObservedValue()[1], + measurements.get(1).getObservedValue()[2]); + final + PV + pv2 = new + PV(measurements.get(1).getDate(), position2, Vector3D.ZERO, 0., 0., 1.); + + final + Vector3D + position3 = new + Vector3D(measurements.get(2).getObservedValue()[0], + measurements.get(2).getObservedValue()[1], + measurements.get(2).getObservedValue()[2]); + final + PV + pv3 = new + PV(measurements.get(2).getDate(), position3, Vector3D.ZERO, 0., 0., 1.); + + // instantiate + the + IOD + method + final + IodGibbs + gibbs = new + IodGibbs(mu); + final + KeplerianOrbit + orbit = gibbs.estimate(frame, pv1, pv2, pv3); + + Assert.assertEquals(context.initialOrbit.getA(), orbit.getA(), 1.0e-9 * context.initialOrbit.getA()); + Assert.assertEquals(context.initialOrbit.getE(), orbit.getE(), 1.0e-9 * context.initialOrbit.getE()); + Assert.assertEquals(context.initialOrbit.getI(), orbit.getI(), 1.0e-9 * context.initialOrbit.getI()); + } + + +if __name__ == '__main__': + suite = unittest.TestLoader().loadTestsFromTestCase(IodGibbsTest) + ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() + sys.exit(ret) diff --git a/python_files/test/KeplerianConverterTest.py b/python_files/test/KeplerianConverterTest.py new file mode 100644 index 0000000000000000000000000000000000000000..885dff70c7d425017caa8061b9398ad5f9e5e09d --- /dev/null +++ b/python_files/test/KeplerianConverterTest.py @@ -0,0 +1,125 @@ +# -*- 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. + */ + +Python version translated from Java by Petrus Hyvönen, SSC 2014 + + """ + +import orekit +orekit.initVM() +from orekit.pyhelpers import setup_orekit_curdir + +from org.orekit.frames import FramesFactory +from org.orekit.propagation.analytical import KeplerianPropagator +from org.orekit.time import AbsoluteDate +from org.hipparchus.geometry.euclidean.threed import Vector3D +from org.orekit.orbits import EquinoctialOrbit +from org.orekit.orbits import Orbit +from org.orekit.orbits import OrbitType +from org.orekit.orbits import PositionAngle +from org.orekit.utils import PVCoordinates +from org.orekit.propagation.conversion import FiniteDifferencePropagatorConverter +from org.orekit.propagation.conversion import KeplerianPropagatorBuilder +import unittest +import sys +from java.util import Arrays + + +class KeplerianConverterTest(unittest.TestCase): + + position = Vector3D(7.0e6, 1.0e6, 4.0e6) + velocity = Vector3D(-500.0, 8000.0, 1000.0) + mu = 3.9860047e14 + + def checkFit(self, orbit, + duration, + stepSize, + threshold, + positionOnly, + expectedRMS, + *args): + + p = KeplerianPropagator(orbit) + + sample = [] + dt = 0.0 + while dt < duration: + sample.append(p.propagate(orbit.getDate().shiftedBy(dt))) + dt += stepSize + + + builder = KeplerianPropagatorBuilder(OrbitType.KEPLERIAN.convertType(orbit), + PositionAngle.MEAN, + 1.0) + + fitter = FiniteDifferencePropagatorConverter(builder, threshold, 1000) + + fitter.convert(Arrays.asList(sample), positionOnly, []) + + self.assertAlmostEqual(fitter.getRMS(), 0.01 * expectedRMS, delta=expectedRMS) + + prop = fitter.getAdaptedPropagator() #(KeplerianPropagator) + fitted = prop.getInitialState().getOrbit() + + eps = 1.0e-12 + self.assertAlmostEqual(orbit.getPVCoordinates().getPosition().getX(), + fitted.getPVCoordinates().getPosition().getX(), + delta = eps * orbit.getPVCoordinates().getPosition().getX()) + self.assertAlmostEqual(orbit.getPVCoordinates().getPosition().getY(), + fitted.getPVCoordinates().getPosition().getY(), + delta = eps * orbit.getPVCoordinates().getPosition().getY()) + self.assertAlmostEqual(orbit.getPVCoordinates().getPosition().getZ(), + fitted.getPVCoordinates().getPosition().getZ(), + delta = eps * orbit.getPVCoordinates().getPosition().getZ()) + + self.assertAlmostEqual(orbit.getPVCoordinates().getVelocity().getX(), + fitted.getPVCoordinates().getVelocity().getX(), + delta = -eps * orbit.getPVCoordinates().getVelocity().getX()) + self.assertAlmostEqual(orbit.getPVCoordinates().getVelocity().getY(), + fitted.getPVCoordinates().getVelocity().getY(), + delta = eps * orbit.getPVCoordinates().getVelocity().getY()) + self.assertAlmostEqual(orbit.getPVCoordinates().getVelocity().getZ(), + fitted.getPVCoordinates().getVelocity().getZ(), + delta = eps * orbit.getPVCoordinates().getVelocity().getZ()) + + def testConversionPositionVelocity(self): + self.checkFit(self.orbit, 86400, 300, 1.0e-3, False, 1.89e-8) + + def testConversionPositionOnly(self): + self.checkFit(self.orbit, 86400, 300, 1.0e-3, True, 2.90e-8) + +# #@Test(expected = OrekitException.class) +# def testConversionWithFreeParameter(self): +# self.checkFit(self.orbit, 86400, 300, 1.0e-3, True, 2.65e-8, "toto"); + + def setUp(self): + setup_orekit_curdir() + + self.initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.) + self.orbit = EquinoctialOrbit(PVCoordinates(self.position, self.velocity), + FramesFactory.getEME2000(), self.initDate, self.mu) + + +if __name__ == '__main__': + + suite = unittest.TestLoader().loadTestsFromTestCase(KeplerianConverterTest) + ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() + sys.exit(ret) diff --git a/python_files/test/NodeDetectorTest.py b/python_files/test/NodeDetectorTest.py new file mode 100644 index 0000000000000000000000000000000000000000..7c7b139dc10e261071ee0bb431902b88948135db --- /dev/null +++ b/python_files/test/NodeDetectorTest.py @@ -0,0 +1,103 @@ +# -*- 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. + */ + +Python version translated from Java by Petrus Hyvönen, SSC 2014 + + """ + +# Python orekit specifics +import orekit +orekit.initVM() + +from orekit.pyhelpers import setup_orekit_curdir +setup_orekit_curdir() # orekit-data.zip shall be in current dir + +from org.orekit.propagation.events import EventsLogger +from org.orekit.propagation.events import NodeDetector +from org.hipparchus.ode.nonstiff import DormandPrince853Integrator +from org.hipparchus.util import FastMath +from org.orekit.frames import FramesFactory +from org.orekit.orbits import KeplerianOrbit +from org.orekit.orbits import PositionAngle +from org.orekit.propagation import SpacecraftState +from org.orekit.propagation.events.handlers import ContinueOnEvent +from org.orekit.propagation.numerical import NumericalPropagator +from org.orekit.time import AbsoluteDate +from org.orekit.time import TimeScalesFactory +from org.orekit.utils import Constants +from orekit import JArray_double + +# Floats are needed to be specific in the orekit interface +a = 800000.0 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS +e = 0.0001 +i = FastMath.toRadians(98.0) +w = -90.0 +raan = 0.0 +v = 0.0 + +inertialFrame = FramesFactory.getEME2000() +initialDate = AbsoluteDate(2014, 1, 1, 0, 0, 0.0, TimeScalesFactory.getUTC()) +finalDate = initialDate.shiftedBy(70*24*60*60.0) +initialOrbit = KeplerianOrbit(a, e, i, w, raan, v, PositionAngle.TRUE, inertialFrame, initialDate, Constants.WGS84_EARTH_MU) +initialState = SpacecraftState(initialOrbit, 1000.0) + +tol = NumericalPropagator.tolerances(10.0, initialOrbit, initialOrbit.getType()) + +# Double array of doubles needs to be retyped to work +integrator = DormandPrince853Integrator(0.001, 1000.0, + JArray_double.cast_(tol[0]), + JArray_double.cast_(tol[1])) + +propagator = NumericalPropagator(integrator) +propagator.setInitialState(initialState) + +# Define 2 instances of NodeDetector: +rawDetector = NodeDetector(1e-6, + initialState.getOrbit(), + initialState.getFrame()).withHandler(ContinueOnEvent().of_(NodeDetector)) + +logger1 = EventsLogger() +node1 = logger1.monitorDetector(rawDetector) +logger2 = EventsLogger() +node2 = logger2.monitorDetector(rawDetector) + +propagator.addEventDetector(node1) +propagator.addEventDetector(node2) + +# First propagation +propagator.setEphemerisMode() +propagator.propagate(finalDate) + +assert 1998==logger1.getLoggedEvents().size() +assert 1998== logger2.getLoggedEvents().size(); +logger1.clearLoggedEvents() +logger2.clearLoggedEvents() + +postpro = propagator.getGeneratedEphemeris() + +# Post-processing +postpro.addEventDetector(node1) +postpro.addEventDetector(node2) +postpro.propagate(finalDate) +assert 1998==logger1.getLoggedEvents().size() +assert 1998==logger2.getLoggedEvents().size() + +print("NodeDetectorTest Successfully run") diff --git a/python_files/test/SmallManeuverAnalyticalModelTest.py b/python_files/test/SmallManeuverAnalyticalModelTest.py new file mode 100644 index 0000000000000000000000000000000000000000..2c9713c9ff5521432a46fe5918006c868f3a6818 --- /dev/null +++ b/python_files/test/SmallManeuverAnalyticalModelTest.py @@ -0,0 +1,226 @@ +# -*- coding: utf-8 -*- + +""" + +/* Copyright 2002-2015 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 2015 + + """ +import orekit + +orekit.initVM() +from orekit.pyhelpers import setup_orekit_curdir + +setup_orekit_curdir() + +from orekit import JArray_double, JArray + +from org.orekit.forces.maneuvers import SmallManeuverAnalyticalModel, ConstantThrustManeuver + +from org.hipparchus.geometry.euclidean.threed import Vector3D +from org.hipparchus.ode.nonstiff import DormandPrince853Integrator +from org.hipparchus.util import FastMath +from org.orekit.utils import Constants +from org.orekit.attitudes import LofOffset + +from org.orekit.frames import FramesFactory +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.numerical import NumericalPropagator +from org.orekit.time import AbsoluteDate +from org.orekit.time import DateComponents +from org.orekit.time import TimeComponents +from org.orekit.time import TimeScalesFactory +from org.orekit.utils import Constants +from org.orekit.utils import PVCoordinates + +from math import radians + +# from org.orekit.forces.maneuvers import getEphemeris + +import unittest +import sys + + +class SmallManeuverAnalyticalModelTest(unittest.TestCase): + def testLowEarthOrbit1(self): + leo = CircularOrbit(7200000.0, -1.0e-5, 2.0e-4, + radians(98.0), + radians(123.456), + 0.0, PositionAngle.MEAN, + FramesFactory.getEME2000(), + AbsoluteDate(DateComponents(2004, 1, 1), + TimeComponents(23, 30, 00.000), + TimeScalesFactory.getUTC()), + Constants.EIGEN5C_EARTH_MU) + mass = 5600.0 + t0 = leo.getDate().shiftedBy(1000.0) + dV = Vector3D(-0.01, 0.02, 0.03) + f = 20.0 + isp = 315.0 + withoutManeuver = self.getEphemeris(leo, mass, t0, Vector3D.ZERO, f, isp) + withManeuver = self.getEphemeris(leo, mass, t0, dV, f, isp) + model = SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0), dV, isp) + + self.assertEquals(t0.toString(), model.getDate().toString()) + + t = withoutManeuver.getMinDate() + while t.compareTo(withoutManeuver.getMaxDate()) < 0: + pvWithout = withoutManeuver.getPVCoordinates(t, leo.getFrame()) + pvWith = withManeuver.getPVCoordinates(t, leo.getFrame()) + pvModel = model.apply(withoutManeuver.propagate(t)).getPVCoordinates(leo.getFrame()) + nominalDeltaP = PVCoordinates(pvWith, pvWithout).getPosition().getNorm() + modelError = PVCoordinates(pvWith, pvModel).getPosition().getNorm() + if t.compareTo(t0) < 0: + # before maneuver, all positions should be equal + self.assertEquals(0, nominalDeltaP, 1.0e-10) + self.assertEquals(0, modelError, 1.0e-10) + else: + # after maneuver, model error should be less than 0.8m, + # despite nominal deltaP exceeds 1 kilometer after less than 3 orbits + if t.durationFrom(t0) > 0.1 * leo.getKeplerianPeriod(): + self.assertTrue(modelError < 0.009 * nominalDeltaP) + + self.assertTrue(modelError < 0.8) + + t = t.shiftedBy(60.0) + + def testLowEarthOrbit2(self): + + leo = CircularOrbit(7200000.0, -1.0e-5, 2.0e-4, + radians(98.0), + radians(123.456), + 0.0, PositionAngle.MEAN, + FramesFactory.getEME2000(), + AbsoluteDate(DateComponents(2004, 1, 1), + TimeComponents(23, 30, 00.000), + TimeScalesFactory.getUTC()), + Constants.EIGEN5C_EARTH_MU) + mass = 5600.0 + t0 = leo.getDate().shiftedBy(1000.0) + dV = Vector3D(-0.01, 0.02, 0.03) + f = 20.0 + isp = 315.0 + withoutManeuver = self.getEphemeris(leo, mass, t0, Vector3D.ZERO, f, isp) + withManeuver = self.getEphemeris(leo, mass, t0, dV, f, isp) + model = SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0), dV, isp) + self.assertEquals(t0.toString(), model.getDate().toString()) + + t = withoutManeuver.getMinDate() + while t.compareTo(withoutManeuver.getMaxDate()) < 0: + pvWithout = withoutManeuver.getPVCoordinates(t, leo.getFrame()) + pvWith = withManeuver.getPVCoordinates(t, leo.getFrame()) + pvModel = model.apply(withoutManeuver.propagate(t).getOrbit()).getPVCoordinates(leo.getFrame()) + nominalDeltaP = PVCoordinates(pvWith, pvWithout).getPosition().getNorm() + modelError = PVCoordinates(pvWith, pvModel).getPosition().getNorm() + if t.compareTo(t0) < 0: + # before maneuver, all positions should be equal + self.assertEquals(0, nominalDeltaP, 1.0e-10) + self.assertEquals(0, modelError, 1.0e-10) + else: + # after maneuver, model error should be less than 0.8m, + # despite nominal deltaP exceeds 1 kilometer after less than 3 orbits + if t.durationFrom(t0) > 0.1 * leo.getKeplerianPeriod(): + self.assertTrue(modelError < 0.009 * nominalDeltaP) + + self.assertTrue(modelError < 0.8) + + t = t.shiftedBy(60.0) + + def testEccentricOrbit(self): + + heo = KeplerianOrbit(90000000.0, 0.92, FastMath.toRadians(98.0), + radians(12.3456), + radians(123.456), + radians(1.23456), PositionAngle.MEAN, + FramesFactory.getEME2000(), + AbsoluteDate(DateComponents(2004, 1, 1), + TimeComponents(23, 30, 00.000), + TimeScalesFactory.getUTC()), + Constants.EIGEN5C_EARTH_MU) + mass = 5600.0 + t0 = heo.getDate().shiftedBy(1000.0) + dV = Vector3D(-0.01, 0.02, 0.03) + f = 20.0 + isp = 315.0 + withoutManeuver = self.getEphemeris(heo, mass, t0, Vector3D.ZERO, f, isp) + withManeuver = self.getEphemeris(heo, mass, t0, dV, f, isp) + model = SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0), dV, isp) + + self.assertEquals(t0.toString(), model.getDate().toString()) + + t = withoutManeuver.getMinDate() + while t.compareTo(withoutManeuver.getMaxDate()) < 0: + pvWithout = withoutManeuver.getPVCoordinates(t, heo.getFrame()) + pvWith = withManeuver.getPVCoordinates(t, heo.getFrame()) + pvModel = model.apply(withoutManeuver.propagate(t)).getPVCoordinates(heo.getFrame()) + nominalDeltaP = PVCoordinates(pvWith, pvWithout).getPosition().getNorm() + modelError = PVCoordinates(pvWith, pvModel).getPosition().getNorm() + if t.compareTo(t0) < 0: + # before maneuver, all positions should be equal + self.assertEquals(0, nominalDeltaP, 1.0e-10) + self.assertEquals(0, modelError, 1.0e-10) + else: + # after maneuver, model error should be less than 1700m, + # despite nominal deltaP exceeds 300 kilometers at perigee, after 3 orbits + if t.durationFrom(t0) > 0.01 * heo.getKeplerianPeriod(): + self.assertTrue(modelError < 0.005 * nominalDeltaP) + + self.assertTrue(modelError < 1700) + t = t.shiftedBy(600.0) + + # Jacobian test removed due to the large amount of manual work currently with 2D java arrays + + def getEphemeris(self, orbit, mass, t0, dV, f, isp): + + law = LofOffset(orbit.getFrame(), LOFType.LVLH) + initialState = SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass) + + # set up numerical propagator + dP = 1.0 + tolerances = NumericalPropagator.tolerances(dP, orbit, orbit.getType()) + + integrator = DormandPrince853Integrator(0.001, 1000.0, JArray_double.cast_(tolerances[0]), + JArray_double.cast_(tolerances[1])) + + integrator.setInitialStepSize(orbit.getKeplerianPeriod() / 100.0) + propagator = NumericalPropagator(integrator) + propagator.setOrbitType(orbit.getType()) + propagator.setInitialState(initialState) + propagator.setAttitudeProvider(law) + + if dV.getNorm() > 1.0e-6: + # set up maneuver + vExhaust = Constants.G0_STANDARD_GRAVITY * isp + dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust) + maneuver = ConstantThrustManeuver(t0, dt, f, isp, dV.normalize()) + propagator.addForceModel(maneuver) + + propagator.setEphemerisMode() + propagator.propagate(t0.shiftedBy(5 * orbit.getKeplerianPeriod())) + return propagator.getGeneratedEphemeris() + + +if __name__ == '__main__': + suite = unittest.TestLoader().loadTestsFromTestCase(SmallManeuverAnalyticalModelTest) + ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() + sys.exit(ret) diff --git a/python_files/test/TLEConverterTest.py b/python_files/test/TLEConverterTest.py new file mode 100644 index 0000000000000000000000000000000000000000..83ab77b760d745302987717c97649a5806bef75e --- /dev/null +++ b/python_files/test/TLEConverterTest.py @@ -0,0 +1,141 @@ +# -*- 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.assertEquals(1, drivers.size()) + self.assertEquals("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() + + self.assertAlmostEquals(expectedRMS, fitter.getRMS(), delta=0.001 * expectedRMS) + + self.assertEquals(tle.getSatelliteNumber(), fitted.getSatelliteNumber()) + self.assertEquals(tle.getClassification(), fitted.getClassification()) + self.assertEquals(tle.getLaunchYear(), fitted.getLaunchYear()) + self.assertEquals(tle.getLaunchNumber(), fitted.getLaunchNumber()) + self.assertEquals(tle.getLaunchPiece(), fitted.getLaunchPiece()) + self.assertEquals(tle.getElementNumber(), fitted.getElementNumber()) + self.assertEquals(tle.getRevolutionNumberAtEpoch(), fitted.getRevolutionNumberAtEpoch()) + + eps = 1.0e-5 + self.assertAlmostEquals(tle.getMeanMotion(), fitted.getMeanMotion(), delta=eps * tle.getMeanMotion()) + self.assertAlmostEquals(tle.getE(), fitted.getE(), delta=eps * tle.getE()) + self.assertAlmostEquals(tle.getI(), fitted.getI(), delta=eps * tle.getI()) + self.assertAlmostEquals(tle.getPerigeeArgument(), fitted.getPerigeeArgument(), + delta=eps * tle.getPerigeeArgument()) + self.assertAlmostEquals(tle.getRaan(), fitted.getRaan(), delta=eps * tle.getRaan()) + self.assertAlmostEquals(tle.getMeanAnomaly(), fitted.getMeanAnomaly(), delta=eps * tle.getMeanAnomaly()) + + if withBStar: + self.assertAlmostEquals(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('regular-data.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/TestAbstractDetector.py b/python_files/test/TestAbstractDetector.py new file mode 100644 index 0000000000000000000000000000000000000000..cf99943c4348c3c23af3e5d4174ade2c247d3e20 --- /dev/null +++ b/python_files/test/TestAbstractDetector.py @@ -0,0 +1,115 @@ +# -*- coding: utf-8 -*- + +import orekit +orekit.initVM() + +from org.orekit.frames import FramesFactory, TopocentricFrame +from org.orekit.bodies import OneAxisEllipsoid, GeodeticPoint +from org.orekit.time import AbsoluteDate, TimeScalesFactory +from org.orekit.orbits import KeplerianOrbit +from org.orekit.utils import Constants +from org.orekit.propagation.analytical import KeplerianPropagator +from org.orekit.utils import PVCoordinates, IERSConventions +from org.orekit.propagation.events.handlers import EventHandler +from org.hipparchus.geometry.euclidean.threed import Vector3D +from org.orekit.python import PythonEventHandler, PythonAbstractDetector +from org.orekit.propagation.events.handlers import ContinueOnEvent, StopOnEvent + +from math import radians +import math +import unittest +import sys + +from orekit.pyhelpers import setup_orekit_curdir +setup_orekit_curdir() + +class PassCounter(PythonEventHandler): + """Eventhandler that counts positive events""" + passes = 0 + + def eventOccurred(self, s, T, increasing): + if increasing: + self.passes = self.passes + 1 + + return EventHandler.Action.CONTINUE + + def resetState(self, detector, oldState): + return oldState + + +class MyElevationDetector(PythonAbstractDetector): + + def __init__(self, elevation, topo, handler=None): + self.elevation = elevation + self.topo = topo + + dmax = float(PythonAbstractDetector.DEFAULT_MAXCHECK) + dthresh = float(PythonAbstractDetector.DEFAULT_THRESHOLD) + dmaxiter = PythonAbstractDetector.DEFAULT_MAX_ITER + if handler is None: + handler = StopOnEvent().of_(MyElevationDetector) + + super(MyElevationDetector, self).__init__(dmax, dthresh, dmaxiter, handler) #super(maxCheck, threshold, maxIter, handler); + + def init(self, *args, **kw): + pass + + def getElevation(self): + return self.elevation + + def getTopocentricFrame(self): + return self.topo + + def g(self, s): + tmp = self.topo.getElevation(s.getPVCoordinates().getPosition(), s.getFrame(), s.getDate())-self.elevation + return tmp + + def create(self, newMaxCheck, newThreshHold, newMaxIter, newHandler): + return MyElevationDetector(self.elevation, self.topo, handler=newHandler) + +class AbstractDetectorTest(unittest.TestCase): + + def testOwnElevationDetector(self): + + initialDate = AbsoluteDate(2014, 1, 1, 23, 30, 00.000, TimeScalesFactory.getUTC()) + inertialFrame = FramesFactory.getEME2000() # inertial frame for orbit definition + position = Vector3D(-6142438.668, 3492467.560, -25767.25680) + velocity = Vector3D(505.8479685, 942.7809215, 7435.922231) + pvCoordinates = PVCoordinates(position, velocity) + initialOrbit = KeplerianOrbit(pvCoordinates, + inertialFrame, + initialDate, + Constants.WGS84_EARTH_MU) + + kepler = KeplerianPropagator(initialOrbit) + + ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, True) + earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, + Constants.WGS84_EARTH_FLATTENING, + ITRF) + + # Station + longitude = radians(45.0) + latitude = radians(25.0) + altitude = 0 + station1 = GeodeticPoint(latitude, longitude, float (altitude)) + sta1Frame = TopocentricFrame(earth, station1, "station 1") + + elevation = math.radians(5.0) + + detector = MyElevationDetector(elevation, sta1Frame) + + mycounter = PassCounter().of_(MyElevationDetector) + detector = detector.withHandler(mycounter) + + kepler.addEventDetector(detector) + + finalState = kepler.propagate(initialDate.shiftedBy(60*60*24.0*15)) + + print(mycounter.passes) + self.assertEquals(52, mycounter.passes) + +if __name__ == '__main__': + suite = unittest.TestLoader().loadTestsFromTestCase(AbstractDetectorTest) + ret = not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() + sys.exit(ret) diff --git a/python_files/test/orekit-data.zip b/python_files/test/orekit-data.zip new file mode 100644 index 0000000000000000000000000000000000000000..d77c38d323d9ab92923e795dfdbe58f80407d902 Binary files /dev/null and b/python_files/test/orekit-data.zip differ diff --git a/python_files/test/pytest.ini b/python_files/test/pytest.ini new file mode 100644 index 0000000000000000000000000000000000000000..e0bdf3d18bc3b8d18f63bcce21b2a05e770b389e --- /dev/null +++ b/python_files/test/pytest.ini @@ -0,0 +1,6 @@ +# content of pytest.ini +# can also be defined in tox.ini or setup.cfg file, although the section +# name in setup.cfg files should be "tool:pytest" + +[pytest] +python_files=*.py diff --git a/python_files/test/regular-data.zip b/python_files/test/regular-data.zip new file mode 100644 index 0000000000000000000000000000000000000000..a3a45540349050f055d990871105f7b4367f1f47 Binary files /dev/null and b/python_files/test/regular-data.zip differ