[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Orekit Users] DSSTPropagator and NumericalPropagator differences
Elisabet <elisabet.cid-borobia@thalesaleniaspace.com> a écrit :
Hello,
Hi Elisabet,
I have been using your tutorial which compares DSSTPropagation with
NumericalPropagation (DSSTPropagation.java). The between the
different orbital
parameters seems ok but the right ascending node presents a growing drift
along the time. Is it normal?
Sorry for the delay, I have been quite busy these days.
I think the behaviour you observe is normal given the configuration.
DSST handles tesseral terms using a regular rotation model for Earth,
which is specified using a single angular rotation rate. In your input
file, you selected the standard WGS84 rate:
central.body.rotation.rate = 7.292115e-5
The numerical propagator, on the other hand relies on the body frame
to get Earth orientation. Accurate body frames will take a lot of
irregularities into account (precession, nutation, Earth Orientation
Parameters, tides) and the rotation axis slightly changes, it is not
fixed.
This implies that both the inertial frame and body frame settings in
the dsst-propagation.in input file from the DSST tutorial have an
effect when comparing DSST with respect to numerical propagator. Slight
offsets of Earth orientation (both axis and phase angle) will introduce
drifts.
I changed the frames to
body.frame = GTOD/2010 simple EOP
and this drastically reduced the drift (which I observed in inclination
rather than ascending node by the way).
Beware! If you change this setting right now with the current Orekit
tutorial, you will get an error stating that the body frame is not
an Earth frame. This is due to a bug in the tutorial, I will fix it in
the upcoming minutes. If you want to give it a try, just apply the
following patch:
diff --git a/src/tutorials/java/fr/cs/examples/KeyValueFileParser.java
b/src/tutorials/java/fr/cs/examples/KeyValueFileParser.java
index f82b8bf..d3b11c9 100644
--- a/src/tutorials/java/fr/cs/examples/KeyValueFileParser.java
+++ b/src/tutorials/java/fr/cs/examples/KeyValueFileParser.java
@@ -438,7 +438,8 @@ public class KeyValueFileParser<Key extends Enum<Key>> {
// check the name against predefined frames
for (Predefined predefined : Predefined.values()) {
if (frameName.equals(predefined.getName())) {
- if (predefined.toString().startsWith("ITRF")) {
+ if (predefined.toString().startsWith("ITRF") ||
+ predefined.toString().startsWith("GTOD")) {
return FramesFactory.getFrame(predefined);
} else {
throw new OrekitException(NOT_EARTH_FRAME, frameName);
Also I would not recommend using any frame with "accurate EOP" in this
tutorial.
Accurate EOP means that the smooth EOP read from IERS files are post-processed
to add high frequency tide effects at sub-daily rates. These effects change
Earth orientation by about 20mm only, so they are meaningful only for geodesy,
not really for traditional space flight dynamics, and they are really, really
heavy to compute. they slow down computation by a huge factor. I have seen a
20 times slowdown in an application that needed a lot of frames transforms.
As in your case you propagate at a future date (2019), you event do not have
the EOP available at that time, so adding the 20mm tide effect when you don't
have the up to 475m effect of DUT1 seems a waste of computing time.
Another thing that may be worth changing is the propagation duration. This is
due to an inefficiency in the tutorial. When duration exceeds about 8-9 days,
the computation speed slow down a lot as memory gets full (the tutorial keeps
all the short periodic terms in memory throughout computation). At the 8-9
days level with your current force model setting, memory consumption reaches
1.6 Gbytes and the classical settings of the Java virtual machines starts to
do a lot of garbage collection more and more often. At some points, it is
triggered almost at each step. Reaching 34 days propagation would take
forever.
Note that this effect is not due to DSST, it is a specific feature of the
tutorial and the way it is written. There is clearly a lot of room for
improvement there, but we did not start optimizing DSST yet. We will do this
once we are sure our implementation is fully validated.
Hope this helps,
Luc
Here attached I send you my input file. As you would see is a LEO SSO orbit.
## Input file for DSSTPropagation
## The input file syntax is a set of key=value lines.
## Blank lines and lines starting with '#' (after whitespace trimming) are
## silently ignored.
## The equal sign may be surrounded by space characters.
## Keys must correspond to the ParameterKey enumerate constants, given that
## matching is not case sensitive and that '_' characters may appear as '.'
## characters in the file.
## This file must contain one orbit defined as keplerian, equinoctial,
circular
## or cartesian.
## Some parameters are optional, default values are shown below between [].
## All dates are treated in UTC timescale.
## The inertial frame for orbit definition and propagation is EME2000.
## Physical data are read from the
src/tutorial/resources/tutorial-orekit-data
## directory.
### Frames definition
# body, choose one from the list (default is CIO/2010-based ITRF simple EOP)
# The "accurate EOP" ones are not really adapted for this
# program, they are more related to geodesy when ground accuracy
# below 20mm is needed
# body.frame = CIO/1996-based ITRF simple EOP
# body.frame = CIO/1996-based ITRF accurate EOP
# body.frame = CIO/2003-based ITRF simple EOP
# body.frame = CIO/2003-based ITRF accurate EOP
# body.frame = CIO/2010-based ITRF simple EOP
body.frame = CIO/2010-based ITRF accurate EOP
# body.frame = Equinox/1996-based ITRF simple EOP
# body.frame = Equinox/1996-based ITRF accurate EOP
# body.frame = Equinox/2003-based ITRF simple EOP
# body.frame = Equinox/2003-based ITRF accurate EOP
# body.frame = Equinox/2010-based ITRF simple EOP
# body.frame = Equinox/2010-based ITRF accurate EOP
# body.frame = TIRF simple EOP
# body.frame = TIRF accurate EOP
# body.frame = TIRF simple EOP
# body.frame = TIRF accurate EOP
# body.frame = TIRF simple EOP
# body.frame = TIRF accurate EOP
# body.frame = GTOD/1996 simple EOP
# body.frame = GTOD/1996 accurate EOP
# body.frame = GTOD/2003 simple EOP
# body.frame = GTOD/2003 accurate EOP
# body.frame = GTOD/2010 simple EOP
# body.frame = GTOD/2010 accurate EOP
# Inertial frame, choose one from the list (default is EME2000).
# The "accurate EOP" ones are not really adapted for this
# program, they are more related to geodesy when ground accuracy
# below 20mm is needed
# inertial.frame = GCRF
# inertial.frame = TOD/1996 simple EOP
# inertial.frame = TOD/1996 accurate EOP
# inertial.frame = TOD/2003 simple EOP
# inertial.frame = TOD/2003 accurate EOP
# inertial.frame = TOD/2010 simple EOP
# inertial.frame = TOD/2010 accurate EOP
# inertial.frame = MOD/1996
# inertial.frame = MOD/2003
# inertial.frame = MOD/2010
inertial.frame = EME2000
# inertial.frame = CIRF/1996 simple EOP
# inertial.frame = CIRF/1996 accurate EOP
# inertial.frame = CIRF/2003 simple EOP
# inertial.frame = CIRF/2003 accurate EOP
# inertial.frame = CIRF/2010 simple EOP
# inertial.frame = CIRF/2010 accurate EOP
# inertial.frame = VEIS1950
### Orbit definition
## date of the orbital parameters (UTC)
orbit.date = 2019-03-01T22:30:00.000
### Keplerian elements
## Semi-major axis (km)
orbit.keplerian.a = 7050.177655
## Eccentricity
orbit.keplerian.e = 0.001229434
## Inclination (degrees)
orbit.keplerian.i = 98.00760166
## Right Ascension of Ascending Node (degrees)
orbit.keplerian.raan = -20.81071488
## Perigee Argument (degrees)
orbit.keplerian.pa = 66.47954373
## Anomaly (degrees)
orbit.keplerian.anomaly = -66.71320157
### Equinoctial elements
## Semi-major axis (km)
# orbit.equinoctial.a = 0.0
## ex/k component of eccentricity vector
# orbit.equinoctial.ex = 0.0
## ey/h component of eccentricity vector
# orbit.equinoctial.ey = 0.0
## hx/q component of inclination vector
# orbit.equinoctial.hx = 0.0
## hy/p component of inclination vector
# orbit.equinoctial.hy = 0.0
## Longitude Argument (degrees)
# orbit.equinoctial.lambda = 0.0
### Circular elements
## Semi-major axis (km)
# orbit.circular.a = 0.0
## ex component of eccentricity vector
# orbit.circular.ex = 0.0
## ey component of eccentricity vector
# orbit.circular.ey = 0.0
## Inclination (degrees)
# orbit.circular.i = 0.0
## Right Ascension of Ascending Node (degrees)
# orbit.circular.raan = 0.0
## Latitude Argument (degrees)
# orbit.circular.alpha = 0.0
### Angle type for anomaly, alpha or lambda (ECCENTRIC/MEAN/TRUE) [MEAN]
orbit.angle.type = TRUE
### Cartesian elements
## Position along X in inertial frame (km)
# orbit.cartesian.px = 0.0
## Position along Y in inertial frame (km)
# orbit.cartesian.py = 0.0
## Position along Z in inertial frame (km)
# orbit.cartesian.pz = 0.0
## Velocity along X in inertial frame (km/s)
# orbit.cartesian.vx = 0.0
## Velocity along Y in inertial frame (km/s)
# orbit.cartesian.vy = 0.0
## Velocity along Z in inertial frame (km/s)
# orbit.cartesian.vz = 0.0
### Do we consider initial orbital elements as osculating? (true/false)
[false]
initial.orbit.is.osculating = true
### Do we output osculating orbital elements? (true/false) [inherited from
initial.orbit.is.osculating by default]
output.orbit.is.osculating = true
### Force models
## Spacecraft mass (kg) [1000.]
# mass = 1000.
## Central body rotation rate (beware, it is in rad/s because it is the
traditional unit for this) [7.292115e-5]
# the default value (7.292115e-5 rad/s) corresponds to WGS84 standard
central.body.rotation.rate = 7.292115e-5
## Central body gravity potential degree
central.body.degree = 6
## Central body gravity potential order
central.body.order = 6
## short period limits
max.degree.zonal.short.periods = 6
max.degree.tesseral.short.periods = 6
max.order.tesseral.short.periods = 6
max.degree.tesseral.m.dailies.short.periods = 6
max.order.tesseral.m.dailies.short.periods = 6
## 3rd body Moon (true/false) [false]
third.body.moon = false
## 3rd body Sun (true/false) [false]
third.body.sun = false
## Atmospheric drag (true/false) [false]
drag = false
## Drag coefficient
# drag.cd = 2.0
## Drag area (m^2)
# drag.sf = 10.0
## Solar Radiation Pressure (true/false) [false]
solar.radiation.pressure = false
## SRP coefficient
# solar.radiation.pressure.cr = 1.2
## SRP area (m^2)
# solar.radiation.pressure.sf = 25.0
### Simulation parameters
## Start date (UTC) [orbit.date]
# start.date = 2011-12-12T11:57:20.000
## One duration must be set, whether in seconds or in days
## If both durations are set, duration.in.days will overwrite duration
## Duration (seconds)
# duration = 31536000.0
## Duration (days)
duration.in.days = 34.0
## Time step between printed elements (seconds)
output.step = 3600
### Output parameters
## Output Keplerian parameters in output file (true/false) [true]
output.Keplerian = true
## Output equinoctial parameters in output file (true/false) [true]
output.equinoctial = false
## Output Cartesian parameters in output file (true/false) [true]
output.Cartesian = false
## Short period coefficients to output (comma separated list of coefficients
names, or "none", or "all") [none]
## in order to get the list of available coefficients names, one should
perform a first
## run with a value set to all and then look at the output file header. Note
that the
## (a), (h), (k), (p), (q) and (L) parts do not belong to the coefficient
names, they
## describe the 6 columns that correspond to one coefficient.
## Beware that the list is huge ...
output.short.period.coefficients = none
## Fixed integration step size for propagation (seconds) [-1.]
## If fixed.integration.step > 0, the ClassicalRungeKutta integrator
## will be used with the given value for the step size, otherwise
## the adaptative step size integrator DormandPrince853 will be used
# fixed.integration.step = 43200.
## If integration step is not fixed, three parameters are used to
## configure the DormandPrince853 integrator, min and max steps in
## seconds, and a position tolerance in meters (beware, it is not
## really the global error, it correspond to a local estimated error,
## at each step, often difficult to relate to final global error)
min.variable.integration.step = 10.0
max.variable.integration.step = 60.0
position.tolerance.variable.integration.step = 60.0
## interpolation grid specification
## only one of the two following options can be set: either fixed number
## of points or maximum gap. If both options are specified, an error
## will be generated.
## if neither option is set, a fixed number of 3 points will be used
# fixed number of interpolation points at each mean elements integration step
# fixed.number.of.interpolation.points = 3
# maximum gap between interpolation points
max.time.gap.between.interpolation.points = 86400
## Do we want to compare the DSST with the numerical propagator ? [false]
numerical.comparison = true
Best regards,
Elisabet