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

Re: [Orekit Users] Problem with Eckstein-Hechler propagation - 50km off from STK/J2 after 24h



Paul,

Thank you for the support. The reason I stuck with that step is that when I ran a quick test to time my code, I found marginal improvements going from 60 to 120, 240 seconds, as my Fixed Step Handler is what really consumes most of the cpu time and this is determined by this handler's fixed time step (in the order of a few seconds). I did not give too much thought to my choice apart from this little experience and carried on.

However, I tested my code with much larger Runge-Kutta integration steps of a few hours as you suggested and indeed it is more efficient and makes more sense (as I am not interested in the short-periodic motion). Still, for my application it is the computational overhead of the FixedStepHandler that determines cpu time, so unfortunately I won't benefit from any practical reduction in cpu time.

I appreciate your interest and will reach out if I have any questions on DSST.

Kind regards,

Alvaro





Le Jeudi 9 Juin 2016 01:44 CEST, paulcefo <paulcefo@buffalo.edu> a écrit:

> Martin,
>
> Any questions about DSST, please ask.
>
> Also, may I ask the context for the 60/120 second step?  Does this
> include the short-periodic motion?
>
> In several LEO applications, 43,200 seconds is a useful Runge-Kutta 
> stepsize for the DSST mean element equations of motion.
>
> Paul
>
> --
> Dr. Paul J. Cefola
> Consultant in Aerospace Systems, Spaceflight Mechanics, & Astrodynamics
> Adjunct Faculty, Dept. of Mechanical and Aerospace Engineering,
> University at Buffalo (SUNY)
>
> 4 Moonstone Way
> Vineyard Haven, MA 02568
> USA
>
> 508-696-1884 (phone on Martha's Vineyard)
> 978-201-1393 (cell)
>
> paulcefo@buffalo.edu
> paul.cefola@gmail.com
>
> On 06/08/2016 9:11 am, MARTIN SAMPAYO Alvaro wrote:
> > Hi all,
> >
> > In case it may help future users (especially those less proficient on
> > the topic like me), I will detail how I have solved my problem.
> >
> > The difference in my solutions when comparing Eckstein-Hechler (EH)
> > propagation and the J2 model in STK was definitely coming from the
> > orbit definition, because EH considers osculating parameters and not
> > mean parameters. (Drifts of 50km per day)
> >
> > I finally went with the DSST propagator, as for me as a new user it
> > was easier to use DSST specifying that my parameters were mean and not
> > osculating rather than changing the propagator initialization of EH
> > (it is not obvious in the documentation!). I was inspired by the
> > source-code of the DSST tutorial made by CS and available in the
> > repository. If I had had to build my DSST model with just the Javadoc,
> > it would have been much more difficult for me.
> >
> > The DSST propagator is very fast. I actually have faster results than
> > with EH for my application (around 2x). With EH I propagate in master
> > mode and every n [1 to 10] seconds a FixedStepHandler collects data.
> > The same with DSST, where I am using a Runge Kutta integrator with
> > fixed step size of 60/120 seconds and a simple Zonal gravity model.
> > The results so far have been good.
> >
> >
> >
> > Thank you all for your kind help and support, especially Luc, with
> > whom I had the luck to personally discuss recently.
> >
> > Alvaro
> >
> >
> >
> >
> >
> > Le Jeudi 2 Juin 2016 17:01 CEST, MAISONOBE Luc <luc.maisonobe@c-s.fr> a
> > écrit:
> >
> >>
> >> MARTIN SAMPAYO Alvaro <Alvaro.MARTIN-SAMPAYO@etu.isae.fr> a écrit :
> >>
> >> > Hi, thank you for the answers.
> >> >
> >> > *@Christophe*:
> >> > I would like to test J2 only in orekit, but as far as I know
> >> > Eckstein-Hechler directly considers 6 harmonics and I think there is
> >> > no way to make it work as a J2, but I might be wrong.
> >>
> >> You can use a constructor with explicit values for the zonal terms
> >> and simply set the terms after J2 to 0.
> >>
> >> Beware that the constructor uses the un-normalized coefficients
> >> c20, c30, ... These coefficients are the *opposite* of J2.
> >> The amplitude is the same because they are un-normalized here,
> >> just the sign has to be set properly.
> >>
> >> Luc
> >>
> >> > I will look into the differences in orbit definitions, thanks for the hint.
> >> >
> >> >
> >> > *@Paul*:
> >> > My goal is to do a preliminary analysis of coverage, revisit times,
> >> > etc on a time-frame of up to 30 days. I am interested in 'nominal'
> >> > orbits, design orbits such as repeating ground tracks, etc. Hence I
> >> > thought a J2 model was ideal as it is the typical reference for the
> >> > design of these orbits as far as I know. Usually other smaller
> >> > perturbations are adjusted for by orbit control, so I thought very
> >> > accurate models would be counterproductive (for this reason and for
> >> > their longer computational time, as I am also interested in having a
> >> > quick propagator).
> >> >
> >> > While I knew Eckstein-Hechler would be slower than J2, I assumed the
> >> > impact of the extra computations would be small and EH and J2 would
> >> > not diverge that much. So my current situation is that I don't know
> >> > whether my EH propagator is really diverging that much from a J2 or
> >> > if my definition of orbits in STK/orekit is not consistent. I
> >> > considered DSST unnecessarily complex but I am no expert.
> >> >
> >> >
> >> > Thanks both for your help,
> >> >
> >> > Alvaro
> >> >
> >> >
> >> >
> >> > Le Jeudi 2 Juin 2016 15:56 CEST, paulcefo <paulcefo@buffalo.edu> a écrit:
> >> >
> >> >> Alvaro,
> >> >>
> >> >> May I ask you to give an approximate discussion of the accuracy goal and
> >> >> the time spans that you are requiring of the Orekit orbit propagator?
> >> >>
> >> >> The Orekit DSST orbit propagator may be on interest to you.
> >> >>
> >> >> Paul
> >> >>
> >> >> --
> >> >> Dr. Paul J. Cefola
> >> >> Consultant in Aerospace Systems, Spaceflight Mechanics, & Astrodynamics
> >> >> Adjunct Faculty, Dept. of Mechanical and Aerospace Engineering,
> >> >> University at Buffalo (SUNY)
> >> >>
> >> >> 4 Moonstone Way
> >> >> Vineyard Haven, MA 02568
> >> >> USA
> >> >>
> >> >> 508-696-1884 (phone on Martha's Vineyard)
> >> >> 978-201-1393 (cell)
> >> >>
> >> >> paulcefo@buffalo.edu
> >> >> paul.cefola@gmail.com
> >> >>
> >> >> On 06/02/2016 8:27 am, MARTIN SAMPAYO Alvaro wrote:
> >> >> > Hi all,
> >> >> >
> >> >> > This is my first post as I have recently started using orekit for my
> >> >> > master thesis, I hope the community can provide some advice, and I
> >> >> > thank you in advance for that.
> >> >> >
> >> >> > I am planning on propagating using Eckstein-hechler (EH) propagator,
> >> >> > my intention is to have a fast propagation with custom stephandlers
> >> >> > which I have already implemented. I could make-do with a simple J2,
> >> >> > but in orekit I found that EH was the simplest that included Earth
> >> >> > oblateness effect (which I definitely need to consider).
> >> >> >
> >> >> > I run Orekit from Matlab, and have tested my orbit propagation using
> >> >> > Keplerianpropagator against STK successfully. However, when I cross
> >> >> > check the results of orekit's EH propagator with STK's J2 or J4, I see
> >> >> > a huge discordance that I cannot attribute to nominal differences in
> >> >> > the propagators margins of error. I am talking of about 51 km of error
> >> >> > (ground track projection) after just 24 hours.
> >> >> >
> >> >> > Since the 2-body problem matches the results of STK, I assume my model
> >> >> > and my synchronization of orekit and STK are correct. I assume I am
> >> >> > having issues with orekit's EH propagator but I cannot see why or
> >> >> > where. There's some more info below. *Any suggestions will be
> >> >> > appreciated!*
> >> >> >
> >> >> > Thank you!
> >> >> >
> >> >> > Alvaro
> >> >> >
> >> >> > ------------------------------------
> >> >> >
> >> >> > I am attaching two images with the ground tracks after 5 days:
> >> >> >
> >> >> > - EH.png : the propagation of orekit/EH vs STK/J2.
> >> >> > - Kepler.png : the propagation of orekit/Kepler and STK/2-body
> >> >> > propagation
> >> >> >
> >> >> > Here are the parameters of my orbit.
> >> >> >
> >> >> > Propagation for 5 days starts on 2004-jan-01 at 01:30:00 UTC.
> >> >> >
> >> >> > semimajor axis 7578137 m.
> >> >> > eccentricity 0.
> >> >> > inclination 98.8 deg.
> >> >> > argument of perigee 90 deg.
> >> >> > RAAN 0.
> >> >> > true anomaly 0.
> >> >> >
> >> >> >
> >> >> > My initializations (note matlab's syntax to use java objects):
> >> >> >
> >> >> > initialDate = AbsoluteDate(2004, 01, 01, 01, 30, 00.000,...
> >> >> >     TimeScalesFactory.getUTC());
> >> >> >
> >> >> > inertialFrame = FramesFactory.getEME2000();
> >> >> >
> >> >> > orbit = KeplerianOrbit(a, e, i, omega, raan, lv,...
> >> >> >                 PositionAngle.TRUE,...
> >> >> >                 inertialFrame, initialDate, ...
> >> >> >                 Constants.EIGEN5C_EARTH_MU);
> >> >> >
> >> >> > EcksteinHechlerPropagator(orbit,...
> >> >> >                 Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,...
> >> >> >                 Constants.EIGEN5C_EARTH_MU, ...
> >> >> >                 Constants.EIGEN5C_EARTH_C20,...
> >> >> >                 Constants.EIGEN5C_EARTH_C30, ...
> >> >> >                 Constants.EIGEN5C_EARTH_C40,...
> >> >> >                 Constants.EIGEN5C_EARTH_C50,...
> >> >> >                 Constants.EIGEN5C_EARTH_C60);
> >> >> >
> >> >> >
> >> >> > and my transformations to ECEF (these happen in the fixedstephandler,
> >> >> > coded in java):
> >> >> >
> >> >> > this.ECEFframe = FramesFactory.getITRF(IERSConventions.IERS_2010, true)
> >> >> >
> >> >> > ...
> >> >> >
> >> >> > Transform fromJ2000toITRFtransform =
> >> >> > 					currentState.getFrame().getTransformTo(this.ECEFframe,
> >> >> > currentState.getDate());
> >> >> >
> >> >> > 			// Obtain satellite position vector
> >> >> > 			Vector3D satelliteVector = fromJ2000toITRFtransform.transformVector(
> >> >> > 					currentState.getPVCoordinates().getPosition());
> >> >> >
> >> >> > 			satLatitudes.add(satelliteVector.getDelta());
> >> >> > 			satLongitudes.add(satelliteVector.getAlpha() > FastMath.PI ?
> >> >> > satelliteVector.getAlpha()-2*FastMath.PI :
> >> >> > satelliteVector.getAlpha()); // To change representation from [0, 2pi]
> >> >> > to [-pi, pi].
> >>
> >>