Bug #121

osculating to mean elements

Added by Daniel Aguilar Taboada about 5 years ago. Updated over 4 years ago.

Status:ClosedStart date:2012-10-16
Priority:NormalDue date:
Assignee:-% Done:

0%

Category:-
Target version:-

Description

It seems that there is an issue within the OsculatingToMeanElementsConverter class. I was getting some errors and after investigating a bit, here is where I see the problem:

Function "value", within QuadratureComputation class, accepts time offset (sec) as input. However, angular offset (rad) is given as inputs to "quadrature.solve".

issue_121.patch Magnifier (4.81 KB) Thomas Neidhart, 2012-10-16 22:34

OsculatingToMeanElementsConverterTest.java Magnifier (13 KB) Daniel Aguilar Taboada, 2012-10-17 22:51

OsculatingToMeanElementsConverter.java Magnifier (13.9 KB) Daniel Aguilar Taboada, 2012-10-17 22:51

OsculatingToMeanElementsConverter.java Magnifier (5.43 KB) Daniel Aguilar Taboada, 2012-10-30 21:10

OsculatingToMeanElementsConverterTest.java Magnifier (11.1 KB) Daniel Aguilar Taboada, 2012-10-30 21:10

History

#1 Updated by Thomas Neidhart about 5 years ago

Daniel Aguilar Taboada wrote:

It seems that there is an issue within the OsculatingToMeanElementsConverter class. I was getting some errors and after investigating a bit, here is where I see the problem:

Function "value", within QuadratureComputation class, accepts time offset (sec) as input. However, angular offset (rad) is given as inputs to "quadrature.solve".

Hi Daniel,

could you provide a simple testcase to illustrate the problems you encountered?
I did some investigations myself and can confirm the observation: the integration interval is defined in radians, while the QuadratureComputation class indeed expects the x value in seconds (for the propagation).

I tried to adapt the class with the following changes:

  • convert the input radians into seconds taking the keplerian period into account
  • adapt the propagation calls from (-delta, delta) to (lat-delta, lat+delta), as the computation would otherwise be out of range of the generated ephemerides this was not visible before, as the integrated interval was much smaller than the actual generated ephemerides due to the wrong unit
  • the tolerances for the Simpson integrator were too strict after the change, resulting in max evaluation exceptions, after looking in other uses changed the estimation of tolerances to this: NumericalPropagator.tolerances(1.0, state.getOrbit(), OrbitType.EQUINOCTIAL);
  • the resulting l value was exactly twice the expected value, thus changing the creation of the EquinoctialOrbit to the following return new SpacecraftState(new EquinoctialOrbit(aMean, results[0], results[1], results[2], results[3], results[4] / 2, PositionAngle.ECCENTRIC, state.getFrame(), state.getDate(), state.getMu())); The unit tests run through with one exception: java.lang.AssertionError: expected: but was: at org.junit.Assert.fail(Assert.java:91) at org.junit.Assert.failNotEquals(Assert.java:645) at org.junit.Assert.assertEquals(Assert.java:441) at org.junit.Assert.assertEquals(Assert.java:510) at org.orekit.propagation.OsculatingToMeanElementsConverterTest.checkResult(OsculatingToMeanElementsConverterTest.java:161) at org.orekit.propagation.OsculatingToMeanElementsConverterTest.testLowEarthOrbit(OsculatingToMeanElementsConverterTest.java:106) In the paper of H.G. Walter about Converting osculating to mean elements ([[http://adsabs.harvard.edu/full/1967AJ.....72..994W]]) it is written that the error for near-earth satellites may be of magnitude 10e-3 so this may not be wrong after all, but I am not sure.

Thomas

#2 Updated by Daniel Aguilar Taboada about 5 years ago

Hi Thomas,

Thanks for your reply.

Find attached a modified OsculatingToMeanElementsConverterTest.java which includes a new test ("testLowEarthOrbit2") where I try to compute mean elements averaged over 27 days (Metop repeat cycle).

I have used the attached OsculatingToMeanElementsConverter.java. Suggests you compare it against that in git repository to see my changes. Note that I believe that your proposed change "from (-delta, delta) to (lat-delta, lat+delta)" is not needed. I have used other of your suggestions. Thanks for that.

In this testLowEarthOrbit2, I get TooManyEvaluationsException. I even got that sometimes averaging over 1 single day.

After all, I have some comments:
*I believe that the eccentric argument of latitude (AoL) should be computed in a different way because of the discontinuity 0-360deg. For instance, using state.getLM() or fitting the evolution of the eccentric AoL to a linear evolution and retrieving the value at the required epoch.

*Why in the value function (that belonging to the QuadratureComputation class) returns "elmt * (1 - h * sinF - k * cosF)" and not simply "elmt". Is that mistake or am I missing something? I tried to find an answer in the SST theory and in the paper you referred before but failed miserably :(

*Thinking about the SimpsonIntegration and the TooManyEvaluationsException. Do you think it would be OK to do the arithmetic mean (except for the AoL) over a given period of time and given time step (new input?)

Many thanks,
Daniel

#3 Updated by Thomas Neidhart about 5 years ago

Daniel Aguilar Taboada wrote:

Hi Thomas,

I have used the attached OsculatingToMeanElementsConverter.java. Suggests you compare it against that in git repository to see my changes. Note that I believe that your proposed change "from (-delta, delta) to (lat-delta, lat+delta)" is not needed. I have used other of your suggestions. Thanks for that.

Hi Daniel,

thanks for the testcase. You are right with the latitude, it is not needed after all as the propagator is being initialized with the current orbit state imho.

In this testLowEarthOrbit2, I get TooManyEvaluationsException. I even got that sometimes averaging over 1 single day.

After all, I have some comments:
*I believe that the eccentric argument of latitude (AoL) should be computed in a different way because of the discontinuity 0-360deg. For instance, using state.getLM() or fitting the evolution of the eccentric AoL to a linear evolution and retrieving the value at the required epoch.

*Why in the value function (that belonging to the QuadratureComputation class) returns "elmt * (1 - h * sinF - k * cosF)" and not simply "elmt". Is that mistake or am I missing something? I tried to find an answer in the SST theory and in the paper you referred before but failed miserably :(

*Thinking about the SimpsonIntegration and the TooManyEvaluationsException. Do you think it would be OK to do the arithmetic mean (except for the AoL) over a given period of time and given time step (new input?)

I will further look into these exceptions. For the question regarding the value function, I do not have an answer right now.

Thomas

#4 Updated by Thomas Neidhart about 5 years ago

I looked further into the TooManyEvaluationExceptions:

There are several places which use a SimpsonIntegrator for calculating the integral over the defined interval [-delta, delta]. Now that the delta is much larger than before, the integrator has to use many more evaluation points to reach the required tolerance level. For each iteration, the number of points is doubled. In the provided test case with 27 averaging revolutions, 1024 points are used for example. So I would suggest to adjust the MAX_EVALUATION parameter according to the number of revolutions that are requested.

#5 Updated by Romain Di Costanzo about 5 years ago

Thomas Neidhart wrote:

I looked further into the TooManyEvaluationExceptions:

There are several places which use a SimpsonIntegrator for calculating the integral over the defined interval [-delta, delta]. Now that the delta is much larger than before, the integrator has to use many more evaluation points to reach the required tolerance level. For each iteration, the number of points is doubled. In the provided test case with 27 averaging revolutions, 1024 points are used for example. So I would suggest to adjust the MAX_EVALUATION parameter according to the number of revolutions that are requested.

Hi,

Here are some informations about the method used in the OsculatingToMeanElementsConverter : The method is coming from the System Description And User's Guide for the GTDS R&D Averaged Orbit Generator, kindly furnished by P. Cefola. Unfortunately this paper is not in public domain and I can't post it. I'll try anyway to give some more informations about the general method.
First of all, you are all right, there was a misunderstanding about times and angles, used for integration.

About the averaging process :
The numerical osculating to mean elements computes initial mean equinoctial orbital elements by averaging osculating equinoctial elements over a time interval equal to one or more satellite revolutions. This integration is done numerically over a time centred on the epoch time, and shifted by the MEAN period of the satellite computed at epoch, times the number of satellite revolutions (t +- nT/2) .
This is why the first operation consists to compute the mean semimajor axis, needed to get the mean satellite's period, with the osculating semimajor axis at epoch used as the initial approximation to the mean semimajorAxis. This computation will trigger osculating parameters computation over the whole period, and generate ephemeris that will be used for any further computation (by interpolating data).
Once we get the mean semimajor axis computed(and the associated mean period), it's then possible to evaluate the averaging integral. The initial integral, defined with temporal bounds, is converted to a new integral with eccentric longitude bounds (F +- n*Pi, where F is the eccentric longitude, and n the . Those new variables lead to the (1 - h sinF - k cosF) factor.
I do not know actually why the temporal integral is replaced by an integration over the eccentric longitude (integration stability ?).
Anyway, reading back the article and comparing it with the code seems to turn out some incoherences. I need to have a deeper look inside the code, I will try to this next week,
Romain

#6 Updated by Romain Di Costanzo about 5 years ago

Romain Di Costanzo wrote:

Thomas Neidhart wrote:

I looked further into the TooManyEvaluationExceptions:

There are several places which use a SimpsonIntegrator for calculating the integral over the defined interval [-delta, delta]. Now that the delta is much larger than before, the integrator has to use many more evaluation points to reach the required tolerance level. For each iteration, the number of points is doubled. In the provided test case with 27 averaging revolutions, 1024 points are used for example. So I would suggest to adjust the MAX_EVALUATION parameter according to the number of revolutions that are requested.

Hi,

Here are some informations about the method used in the OsculatingToMeanElementsConverter : The method is coming from the System Description And User's Guide for the GTDS R&D Averaged Orbit Generator, kindly furnished by P. Cefola. Unfortunately this paper is not in public domain and I can't post it. I'll try anyway to give some more informations about the general method.
First of all, you are all right, there was a misunderstanding about times and angles, used for integration.

About the averaging process :
The numerical osculating to mean elements computes initial mean equinoctial orbital elements by averaging osculating equinoctial elements over a time interval equal to one or more satellite revolutions. This integration is done numerically over a time centred on the epoch time, and shifted by the MEAN period of the satellite computed at epoch, times the number of satellite revolutions (t +- nT/2) .
This is why the first operation consists to compute the mean semimajor axis, needed to get the mean satellite's period, with the osculating semimajor axis at epoch used as the initial approximation to the mean semimajorAxis. This computation will trigger osculating parameters computation over the whole period, and generate ephemeris that will be used for any further computation (by interpolating data).
Once we get the mean semimajor axis computed(and the associated mean period), it's then possible to evaluate the averaging integral. The initial integral, defined with temporal bounds, is converted to a new integral with eccentric longitude bounds (F +- n*Pi, where F is the eccentric longitude, and n the . Those new variables lead to the (1 - h sinF - k cosF) factor.
I do not know actually why the temporal integral is replaced by an integration over the eccentric longitude (integration stability ?).
Anyway, reading back the article and comparing it with the code seems to turn out some incoherences. I need to have a deeper look inside the code, I will try to this next week,
Romain

I've finally the answers I was looking for. But I'm afraid this method has no future :
The substitution of time by eccentric anomaly in the integral is done because in case of eccentric orbits, orbital parameters rate of change is more important at perigee than at apogee. So in case of using a fixed step integration method, we need to have more evaluations at perigee than at apogee. Using an eccentric longitude is then well adapted. This method then was suited for computation in the 80's (The paper used to implement the method is dated from 1978...).
As we now dispose of variable step integrators, there is no point of using this method.
In fact, we could replace it by doing a simple arithmetic average over orbital parameters, over the number of satellite revolutions given.
A general thinking has been started on this topic. I would suggest not to use this method in current state. It will be deprecated in next version.
Regards,

Romain

#7 Updated by Daniel Aguilar Taboada about 5 years ago

Hi Romain,

Thanks for your response.
For the time being I am using the attached file for my computations (simple regression), in case you find it useful for you.

Thanks again,
Daniel

PS. Attached also the corresponding test java file with some modifications to check that results are OK (not so bad).

#8 Updated by Pascal Parraud about 5 years ago

Hi Daniel,

Thanks for the alternative for osculating to mean elements conversion, your version works better than the original one but IMHO the averaging of the true longitude argument doesn't make sense.

By the way, the associated JUnit tests didn't test well: the first one (testTwoBodyPropagator) should give practically no difference between the initial orbit and the converted one as a numerical propagator without perturbations/forces is just like a keplerian propagator and in this case osculating are mean elements.

I'm working on a new conversion package that will allow to convert from any propagator to any one else with the side effect of converting from osculating to mean elements, it needs more testing but looks promising.

Thanks again,

Pascal

#9 Updated by Daniel Aguilar Taboada about 5 years ago

Hi Pascal,

Thanks for the info!

Regarding the true longitude argument, you may be right. In the end, I am only using the mean inclination.

Regarding the JUnit tests, you are right. I found out that it was due to truncation coming from the SimpleRegression class. I did not have time to look for a better solution :( What I had was enough good for me.

Many thanks again! Looking forward to see the new conversion package :)
Dani

#10 Updated by Pascal Parraud over 4 years ago

Closing issue as the OsculatingToMeanElementsConverter class has been completely revamped thanks to the new conversion package.

If you found an issue with this new class, please open an other issue in this tracker.

Also available in: Atom PDF