[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Orekit Users] CelestialBody Position extrapolation

Hi quentin,

Quentin Nénon <q.nenon@gmail.com> a écrit :


Thank you very much Luc for your answer, it works for Venus ! Thanks also
for the advices. My problem with the propagation of the spacecraft was in
fact due to a frame problem, I was not using the good one...

This may help for the other post, I will see if this enable me to make it
work better. I will then continue the discussion on the approriated topic.
I must also admit that I compared the results of propagation in the Earth
and Sun inertially oriented frames to other results I have using another
software that is already validated. The orbit computed around the Earth
seems to be the good one (or, at least, better than the one around the Sun).
It seems to me the solution for the other problem is different. I improved the results there not by changing frames, but by changing the configuration of the integrator (well, maybe it was also important in the Venus propagation case because I fixed it there also).
The tolerances were much too tight, and reusing exactly the same tolerance for a propagation in Earth vicinity and considering full solar system scale does not work. The tolerances have a relative part and since the order of magnitude of the state vector components are not the same, this relative part must be adjusted differently in the two cases. I have also changed the min/max step sizes.
The attached file can be put in the distribution Junit test folder (in the propagation/numerical package) so it directly uses the distribution data for ephemerides and so on (assuming you use the development version, since I added JPL ephemerides for the 2000 period only two days ago, while fixing issue 171). When you run this, it will create in your home directory a file with the position differences (column 1 is time in seconds the other columns are dx, dy, dz). I have also attached a plot of the result.
You will see the error is now slightly reduced (to about 1400km), but more importantly it does not diverge, it appears to be centered around 0 with orbital period. With the previous settings, the error exploded after a few days with a clear quadratic tendency.
My guess (its only a guess), is that the integrator was completely lost in numerical noise and selected a wrong track then diverged and escaped Earth vicinity. 





Even if the solution you propose solve the problem for Venus and may solve
the issue of the other post, it does not explain why it doesn't work to
propagate in any inertially oriented frame. I expect to have errors with an
order of magnitude of numerical precision.
Don't be too confident on the accuracy of numerical integration of differential equations. It works well, but is sensitive to tuning in difficult cases. In the space flight dynamics domain, we don't often mett such "difficult cases". Typically, we don't have stiff problems like people working in chemical kinetics for example. However, when looking from a Sun point of view to a spacecraft orbiting Earth (at 13000km at start), then we put ourselves in a potential well and ask for accuracy. We ask too much to the integrator.
Another point is I think we should check acceleration composition more thoroughsly in Orekit. I have started some work while looking at your problem, in order to add an acceleration vector to PVCoordinates (this is trivial and done in a fex minutes) and to use it properly everywhere (this is not trivial and a long task). In particular implementing composition in the Transform class needs some care, even if the final equations will be only a few lines long, and should not be done in a hurry. If we manage to do that, we will probably be able to solve this kind of problems more elegantly, by allowing choosing an integration frame almost at will, and even non-inertial frames. This may also help for issue 22 <https://www.orekit.org/forge/issues/22>. 



But let me thank you very much again, you gave me a very fast and effective
answer, and I'm very happy that you find and take some times to answer on
the mailing list !


You're welcome
Luc



Quentin




2014-04-09 11:20 GMT+02:00 MAISONOBE Luc <luc.maisonobe@c-s.fr>:


Hi Quentin,


Quentin Nénon <q.nenon@gmail.com> a écrit :


 Hello Orekit users,


I find this post interesting because Simon tried to extrapolate the
position and velocity of celestial bodies in the solar system. We don't
need to do it as JPL already did it for us for instance, but I tried to
propagate the orbit of Venus (for instance) during 30 days and I have a
difference in the position, compare to the translated JPL ephemerides, of
more than 35,000 km !

I was interested by this when I was comparing results of propagation of
spacecrafts far from Earth to their "actual" positions, given by the
Horizon JPL website. I then decided to do the same comparison on a simple
body, like Venus. You will find enclosed a main class in which I am just
using keplerian attractions from the Sun and the 8 other planets of the
Solar System (assuming Pluto is a planet...). I know that I should add the
solar radiation pressure to the propagator to be more realistic, but this
doesn't work (or gives very strange results, such as Venus at a distance
of
100AU from the Sun, using a sphere model with the Venus radius). Even if I
do not take into account the solar radiation pressure, I do not expect to
have more than 35,000km of difference after 30 days !

I tried to do the same thing propagating around the Sun, but it almost
gives the same results (not exactly, as the propagation around the Sun and
the Earth are not equivalent for the moment, see the
"PropagationFromDifferentPointOfViewPb" post).

Does anyone have an idea of what I did wrong or what is not working ? This
may also help the discussion about the problem of propagation from
different points of view.



Here is a solution I propose to you for this specific problem. I have not
yet checked if this also correct the former problem you already submitted
last before.

What you are attempting here is not really propagating one body around
another one with a few less important perturbing bodies. It is rather
solving a full n-body problem in free space, where all bodies are consider
at the same level.

For this to work, you should rather work in a *fixed* frame, not really
bound to any specific body, and then consider the *absolute* attraction of
all bodies.
For the fixed frame, I suggest to use ICRF (using
FramesFactory.getICRF()), which is at the center of mass of the solar
system. Of course, it is close to Sun since Sun is the more massive body
around, but it is not exactly at the same position (due mainly to Jupiter
and Saturn).

In order to compute the absolute gravity pull from each other body with
respect to your propagated body, you need to build a custom BodyAttraction
force model (you will find one in the new Junit test case I committed a few
minutes ago, <https://www.orekit.org/forge/projects/orekit/repository/
revisions/e59497c723485eebb9e4b3e15100a50fb2a4c107/diff/src/test/java/
org/orekit/bodies/SolarBodyTest.java>, have a look at the
org.orekit.bodies.SolarBodyTest class and more precisely at the
testPropagationVsEphemeris method and BodyAttraction internal class.
BodyAttraction was simply created from the original ThirdBodyAttraction,
but removing the subtraction of the acceleration with respect to the frame,
in order to just compute the absolute acceleration. BEWARE, this works only
because we have selected the computation frame at the center of mass! We
use this force model for all bodies, including Sun (but obviuosly not the
body you propagate of course, which is Venus in your use case).

Now since all the gravity fields are considered in the new force model
instances, and since our central point is an empty point, we don't want to
have a "central" attraction. Unfortunately, Orekit does not allow you to
propagate without a central attraction coefficient. So we use a quick and
dirty hack and set the fictitious Venus orbit around IRF and the mu in the
propagator to be a negligible value, so this impossible to suppress force
will in fact be computed as very small.

All this has been set up in the testPropagationVsEphemeris test case. The
end result is that after one month propagation, the propagator and the
ephemeris are within 1400m to each other (much smaller than 35000 km!).

Could you look if this would also work in your other case?

On a side note, you will see that the new test case which is based on your
own test case has been heavily changed. Here are some random advices I can
give you and which I have applied:

 - don't set tool stringent accuracy in the integrator, propagating
   orbits with 1e-5 m tolerance on position is too much, even for
   spacecraft around Earth, so its even worse for planets. In your
   case, I have put 1000m
 - don't set the tolerances as litteral arrays which cannot be read,
   use the static NumericalPropagator.tolerances method
 - don't clutter the code with very low level computation of array
   lengths, just rely on List instances and let them be allocated
   on the fly, it's simpler to read (of course, if you think you may
   run low in memory because you have huge collections to store,
   computing the size beforehand can still be useful, so my comment
   is not an absolute rule)
 - instead of building arrays in the step handler and process them
   later (in your case to compute the difference) to build yet another
   array, just build the final one on the fly (I have even completely
   removed it in my test case, because I do not plot it, I just do some
   asserts, but if you plot it just fill up the positionDifferenceNorm
   from within the step handler and don't use an intermediate
   pvResultTable
 - replace litteral constants like 86400 with Constants.JULIAN_DAY

Hope this helps,
Luc




Thank you very much !

Quentin


2013-03-29 19:48 GMT+01:00 Spörri Simon (s) <simon.spoerri@students.fhnw.
ch>
:

  Hi Luc,


 thanks a lot for your help on even the most amateurish questions!
it works charmingly now.

 cheers,
simon

  From: MAISONOBE Luc <luc.maisonobe@c-s.fr>
Reply-To: <orekit-users@orekit.org>
Date: Fri, 29 Mar 2013 14:19:07 +0100
To: <orekit-users@orekit.org>
Subject: Re: [Orekit Users] CelestialBody Position extrapolation

  simon.spoerri@students.fhnw.ch a écrit :

 Dear all,


 Hi Simon,


 i'd like to create a discrete sampling of celestial body positions for a
complete orbit of all bodies. I tried both getting positions from the
celestialBody object directly and using a Propagator to get positions of
the
orbits. However i run into ephemerides constraints and cannot obtain
coordinates past around 2030 (Error: out of range date for EARTH_MOON
ephemerides: 2030-05-20T14:36:30.712).


 I'm not sure about your needs. Orekit uses JPL ephemerides (or INPOP
ephemerides) to compute position/velocity of celestial bodies. It does
not do any propagation by itself for these bodies.

 The limitation you get is simply the limitation of the JPL ephemerides
you use. I guess you simply relied on the orekit-data.zip file
provided for convenience in the download section of the forge and did
not change the default set of data contained in this zip archive. If
you want to go past 2030, you can look inside the archive and you will
see it contains a file named unxp1962.406. This is a reduce JPL
ephemeris file. Remove this file and replace it by the original files
from JPL that you will find here:
<ftp://ssd.jpl.nasa.gov/pub/eph/planets/SunOS/de406/>. Each file
contains ephemerides for 300 years. If you retrieve all 20 files, you
will be able to propagate from year -3000 to year +3000. You can look
at the wiki page
<https://www.orekit.org/forge/projects/orekit/wiki/Configuration> for
further information on how Orekit does load such data.

 Please note that JPL ephemerides are already sampled. The correspond
to large sets of Chebyshev polynomials. Adding another sampling layer
on top of that is probably not worth the trouble. Just using the
celestial body itself is probably sufficient for most purposes. The
current version of celestial bodies that use JPL ephemerides is also
thread safe (it was not thread safe in the 5.x series of the Orekit
library), and uses efficient cache mechanisms even when several
threads uses dates very wide apart from each other. This may be useful
in a web application as you cannot know beforehand if two different
clients will request positions around the same dates.

 best regards,
Luc

 I try to obtain the positions in the SolarSystemBarycenter Frame.
My question is if it is possible at all in Orekit to extrapolate
body positions
around 150 years into the future/past? I know that this is not the
main purpose
of Orekit and I do not require very accurate results, nevertheless i
wanted to
use Orekit for that for convenience...
And if it is possible, how should i do it?

 thanks a lot in advance,
Simon





 ----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.










----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.









----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.

package org.orekit.propagation.numerical;

import java.io.File;
import java.io.FileNotFoundException;
import java.io.IOException;
import java.io.PrintStream;
import java.util.Locale;

import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
import org.apache.commons.math3.ode.AbstractIntegrator;
import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
import org.apache.commons.math3.util.FastMath;
import org.junit.Test;
import org.orekit.Utils;
import org.orekit.bodies.CelestialBody;
import org.orekit.bodies.CelestialBodyFactory;
import org.orekit.errors.OrekitException;
import org.orekit.errors.PropagationException;
import org.orekit.forces.gravity.ThirdBodyAttraction;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.orbits.CartesianOrbit;
import org.orekit.orbits.OrbitType;
import org.orekit.propagation.BoundedPropagator;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.sampling.OrekitFixedStepHandler;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.Constants;
import org.orekit.utils.PVCoordinates;

public class propagationFromDifferentPOVPb {

    @Test
	public void testDifferentPOVPb() throws OrekitException, IOException {

        Utils.setDataRoot("regular-data");

        //We want to propagate from dateInit to dateEnd
		final AbsoluteDate dateInit = AbsoluteDate.J2000_EPOCH;
		AbsoluteDate dateEnd = dateInit.shiftedBy(100 * Constants.JULIAN_DAY);
		
		// We load the sun, earth and moon celestial bodies
		CelestialBody sun   = CelestialBodyFactory.getSun();
        CelestialBody earth = CelestialBodyFactory.getEarth();
        CelestialBody moon  = CelestialBodyFactory.getMoon();

		//Definition of the EME2000 Frame for the Earth point of view propagation
		// and definition of the sunFrame for the Sun point of view propagation
		final Frame eme2000 = FramesFactory.getEME2000();
		final Frame sunFrame = sun.getInertiallyOrientedFrame();
		
		//Coordinates of the spacecraft at the dateInit date in the EME2000 and ICRF frames
		PVCoordinates pvEME2000  = new PVCoordinates(new Vector3D(13092474.5414837,
		                                                           4713640.31787418,
		                                                           -194611.161764848),
		                                             new Vector3D( 338.108445229555,
		                                                           661.852963447466,
		                                                           129.590404183624));
		PVCoordinates pvSunFrame = eme2000.getTransformTo(sunFrame, dateInit).transformPVCoordinates(pvEME2000);
		
		//First propagation : we create a SpacecraftState using the EME2000 Frame and the
		// Earth GM and we add the perturbations due to the Sun and the Moon
		SpacecraftState stateEarthPOV = new SpacecraftState(new CartesianOrbit(pvEME2000, eme2000, dateInit, earth.getGM()));

		final double[][] tolE = NumericalPropagator.tolerances(1000.0, stateEarthPOV.getOrbit(), OrbitType.CARTESIAN);
        AbstractIntegrator dopE = new DormandPrince853Integrator(1.0e-3, 1e5, tolE[0], tolE[1]);
		NumericalPropagator propagEarthPOV = new NumericalPropagator(dopE);
		propagEarthPOV.setInitialState(stateEarthPOV);
		propagEarthPOV.setOrbitType(OrbitType.CARTESIAN);
		propagEarthPOV.setEphemerisMode();
		propagEarthPOV.setMu(earth.getGM());
		propagEarthPOV.addForceModel(new ThirdBodyAttraction(sun));
        propagEarthPOV.addForceModel(new ThirdBodyAttraction(moon));
		propagEarthPOV.propagate(dateEnd);
		final BoundedPropagator earthPOVEphemeris = propagEarthPOV.getGeneratedEphemeris();
		
        //Second propagation : we create a SpacecraftState using the Sun frame, Sun mu
        // and then we add the perturbation of all remaining bodies
        SpacecraftState stateSunPOV = new SpacecraftState(new CartesianOrbit(pvSunFrame, sunFrame, dateInit, sun.getGM()));
        final double[][] tolS = NumericalPropagator.tolerances(100.0, stateSunPOV.getOrbit(), OrbitType.CARTESIAN);
        AbstractIntegrator dopS = new DormandPrince853Integrator(1.0e-3, 1e5, tolS[0], tolS[1]);
        NumericalPropagator propagSunPOV = new NumericalPropagator(dopS);
        propagSunPOV.setInitialState(stateSunPOV);
        propagSunPOV.setOrbitType(OrbitType.CARTESIAN);
        propagSunPOV.setMasterMode(900.0, new OrekitFixedStepHandler() {

            PrintStream out;
            double maxError;

            public void init(SpacecraftState s0, AbsoluteDate t) throws PropagationException {
                try {
                    out      = new PrintStream(new File(System.getProperty("user.home"), "earth-sun.dat"));
                    maxError = Double.NEGATIVE_INFINITY;
                } catch (FileNotFoundException fnfe) {
                    throw new PropagationException(fnfe, LocalizedFormats.SIMPLE_MESSAGE, fnfe.getLocalizedMessage());
                }
            }
            
            public void handleStep(SpacecraftState currentState, boolean isLast)
                throws PropagationException {
                try {
                    final Vector3D pE = earthPOVEphemeris.getPVCoordinates(currentState.getDate(), eme2000).getPosition();
                    final Vector3D pS = currentState.getPVCoordinates(eme2000).getPosition();
                    final Vector3D delta = pE.subtract(pS);
                    out.format(Locale.US, "%6.0f %12.5e %12.5e %12.5ee%n",
                               currentState.getDate().durationFrom(dateInit),
                               delta.getX(), delta.getY(), delta.getZ());
                    maxError = FastMath.max(maxError, delta.getNorm());
                    if (isLast) {
                        out.close();
                        System.out.println("max error: " + maxError + " m");
                    }
                } catch (OrekitException oe) {
                    throw new PropagationException(oe);
                }
            }

        });

        propagSunPOV.setMu(sun.getGM());
        propagSunPOV.addForceModel(new ThirdBodyAttraction(earth));
        propagSunPOV.addForceModel(new ThirdBodyAttraction(moon));
        propagSunPOV.propagate(dateEnd);

	}

}

Attachment: earth-sun.png
Description: PNG image