Preliminary Orbit Determination - Gauss Method

Greetings Orekit community,

I was thinking about working with some preliminary orbit determination and I had a quick look at what Orekit has. I might have missed it, but I could not see Gauss’ method among other preliminary orbit determination methods.

Is Gauss method included somewhere else or have you not implemented it on purpose for some other reason?

Thank you for your answer in advance :blush:

Hi @Echulion

The Gauss method is not available in Orekit.
For angular based initial orbit determination you can use Laplace (i.e., or Gooding (i.e., methods.

I recommend you to try the Laplace method because it is close to the Gauss method in terms of accuracy and use.

Please not that the Gooding method is also an angular based IOD method but it requires two first guess for the ranges values. However, it is usually not available.

Best regards,

1 Like


According to the literature (see for example this), Gauss is much more useful than Laplace for Earth-orbiting objects, so I think it would be actually be nice addition to Orekit! There even exists an iterative version of it by the way, available in D. Vallado’s book.
Note that the advantage of both Laplace and Gauss is that as mentionned by Bryan they do not require any guess, so they could even be used to initialize Gooding.


I think it could be useful to open issues on the Orekit’s Gitlab repository to implement these new IOD methods. Gauss method will be easy to implement and validate following Vallado’s book. I think adding the Double-R iteration method could be also interesting.

I found this thesis which is very interesting about IOD.

Hello @bcazabonne , @Serrof

In the meantime, I myself were searching the literature for the accuracy of different IOD methods. As @Serrof mentioned, Gauss method prooves useful in certain cases which can be a nice feature for Orekit. I was also checking Vallado’s book and I think I may adapt his example to implement a Gauss method. I also found this thesis, which may prove useful particularly in this subject.

I will share the Gauss method once I write it, though it will be in Python since I am not competent in Java :sweat_smile: And I will try to open issue in Orekit’s Gitlab Repo if I can create the time.

Thanks for the discussion!


That’s great if you can contribute it! Don’t worry about the Python to Java translation, we can do it :slight_smile:

Best regards,

1 Like

Here is the Jupyter Notebook for the initial orbit determination with Gauss method. I’ve also compared it with Lambert and Gibbs methods.

It may be poor coding but I believe it does the job :nerd_face:

IOD with Gauss Method.ipynb (24.0 KB)


Greate! Thank you!

Could you open an issue in our Gitlab repository with a title like “Add Gauss method for initial orbit determination”.
Also, could you put in the description the link of the current forum thread?

Thank you

1 Like

A small question, for Gauss method they are three coordinates for the observer (i.e. one per measurement). Can the observer be the same (at three different epochs) or Gauss method needs three different observers at three different locations?

Best regards

Hello @bcazabonne,

Okay I will open an issue whenever I can, thanks.

On the other hand,

The code I have shared implements 3 observations from one observing location. But it is also possible and easy to apply different observation locations as well. From the below figure given in the Orbital Mechanics for Engineering Students book of Curtis it can also be seen:

Here, having three observations from same or different locations would not create a problem, since capital R vectors will change either way with the rotation of earth.

Kind regards,

1 Like

Hi there,

I haven’t looked at the code yet but I don’t understand what you mean by comparing it to Lambert and Gibbs. Unlike those Gauss is an angle only method. If you want to compare like for like, you need to use Laplace and Gooding.


Yes you are right, and thank you for the feedback @Serrof . Let me elaborate and fix what I meant there. By using 3 observations and utilizing angles-only Gauss method, we can find the respective position vectors of the satellite.

These three position vectors can later be used in Lambert and Gibbs method to estimate an initial orbit. This was what I meant, and I am sorry if I caused any misunderstanding.

Warning about code: In the code, I tried to test my findings with the results of an example which Vallado shared (Example 7-2). There has been some time since this book was first published so Vallado uses earth parameters slightly different than what we use here in Orekit. Therefore, the capital R vectors we find via this code does not match with Vallado’s results. Although I’ve tested and verified that this is the case, any other user is welcome to test it again. Its a good practice (though time consuming :hourglass:).

After being certain about the source of this discrepancy I just overwritten these R vectors manually to get aligning results in the end.

1 Like

Hello @Echulion,

I am translating what you wrote to resolve issue 982 on the gitlab, i have some questions about your code.
I understand that you followed the Vallado algorithm until the search of v2, where you switched on the algorithm 5.5 of Curtis book. Why did you chose to do that and not follow Vallado suggestion to try to use Gibbs or Lambert methods to find v2 ? Inversely, what not use Curtis algorithm the all time ?
Second i see that you truncated the Taylor serie function of f and g at the first order. Is it approximated like this in Vallado and Curtis but because it seems easier by hand. We could expand the all Taylor serie, or as Vallado suggests “that f and g functions should be used instead of the series coefficients”. What do you think of this ?
Finally, I am translating your code with 3 different positions and not one for now :slight_smile:


Julien A.

Hey @jasquier,

I saw earlier that you’ve started translating the code, its wonderful news! I will try to do my best and address your questions now.

Actually, I implemented both Lambert’s and Gibbs methods at the end of the code as Vallado suggests. The v2 result we get after applying Gibbs method perfectly matches with what Vallado provided in his example.

I cannot rightly remember now but I think that the data Vallado provided in his example was really on point to follow (in terms of verification at every step). So I started with his work. At some point in his example, he mentioned that the Gauss method formally ends once we found slant ranges. And since I’ve also found his v2 result applying Gibbs method, I wanted to see the v2 results of Algorithm 5.5 provided by Curtis. Long story short, I did not have any particular reason to switch between books, but you have every right to question this odd behavior. I would do the same tbh :smile: It would also be interesting to see the results following Algorithm 5.5. from start to end.

Yes, well Curtis ended up with the first order truncation of the series expressions, as he pointed out that if consecutive time intervals of observations are small enough, just first two terms could be retained. He also states that “we want to exclude all terms in f and g beyond the first two so that only the unknown r2 appears in equations”. So this first order truncation greatly reduces the computational work, and you can check it in the section “5.10 Gauss method of preliminary orbit determination” of Curtis.

Lastly, I couldn’t understand what you mean by this though, I am sorry.

Maybe you could also end the code translation once you found the slant ranges, as Vallado pointed out in his example. I think it makes sense that users ending up with slant ranges (or taking one more step, position vectors) using IodGauss class, since continuing further with velocity requires other IOD methods in play such as Gibbs or Lamberts methods.

Finally, it is wonderful to see that you are also improving the code by including 3 different observation positions! I thank you sincerely for your contribution :blush:


Hey guys,

I posted a comment on the gitlab issue but I’m rewriting it here for everyone. It’d be nice to be able to choose between the original, historical method of the Prince of Mathematics a.k.a. Gauss and iterative versions such as the one described in Vallado’s book.

Looking again at the IOD package, I made a scary discovery. I don’t see any Field version (my apologies if I missed it). Is there any particular reason for that appart from lack of time? In that case I will open an issue (if not already there, sorry I didn’t check) as I think it is very important to keep the gap between code for doubles and code for Field to a minimum. Maybe @jasquier since you’ve started working on IOD classes you’d be keen on tackling that? I’m tempted to say that one has not made a proper contribution to Orekit if one has not implemented some Field stuff :rofl:


That’s true, there is no Field implementation of the IOD methods. The reason: nobody asked for it before :slight_smile:

I agree with this comment.

Hey everyone,

I am sorry but I have never used Field implementations (or I don’t know if I used without knowing). Is there any source where I can check it out? If it is related with my field I could give it some thought.

Though my initial guess is, it is something that developers often deal with?

Field is the dark side of Orekit: powerful and dangerous.
Joking aside, fields are an extension of primitive double computations. In 99% of use cases, they are used to compute automatically derivatives alongside values of algorithms.

Many classes in Orekit have their Field counterpart. You can for instance compare CartesianOrbit to FieldCartesianOrbit. Other coordinates have their Field implementation, as well as propagators, measurement models, etc. You’ll see that it’s quite unreadable as unfortunately the Java language does not authorize operator overloading so “+” becomes “add”…

As Luc said, they are meant to represent classes which have their own algebra, the usual case being partial derivatives up to a given order. Check out Gradient and DerivativeStructure in Hipparchus for respectively the first and arbitrary order respectively.

Note that for partial derivatives, the so called Taylor map inversion can provide a workaround to a complete implementation of a given solver. See here for more details.


Hey @luc and @Serrof ,

Thank you for your explanation, I now have a rough image in my mind regarding Field stuff. I will still do my research further whenever I find time.

However, I feel like it is not something for me to worry about :laughing: