Mars Lagrangian Point (L2)

[didimelli]

Hello,

first of all, thank you for developing this library, the api is really great!

I am using Skyfield for a big project at work, and I need to define the L2 for Mars.
What I need to achieve with this project is to get visibilities of Mars and its L2 from some locations on Earth during a given time window.

I believe that the ideal solution would be to have L2 defined as a VectorSum, so that I would be able to use .at(t) methods
(same as Mars is defined when you do something like

planets = load("de421.bsp")
mars: VectorSum = planets["mars"]

)

but i really do not understand how I could achieve this!

Do you have any advice for this kind of problem?
Is my approach (creating L2 as a VectorSum) correct and doable?
Is there any approach that may be more suitable for my problem?

Thank you very much!

[davidmikolas]

One problem is that the whole concept of Lagrange points is based on the circular restricted three body problem (CR3BP or CRTBP) which in this case means Mars and the Sun move around their common center of mass in perfect circular orbits, and there is no other body with any mass anywhere else in the universe.

https://en.wikipedia.org/wiki/Three-body_problem#Restricted_three-body_problem

In this case the mass of the other planets (especially Jupiter) is the secondary problem; Mars' very large eccentricity is the primary one. It's almost 0.1 It doesn't move in a circular orbit so there is no real Lagrange point at all.

However, if you want a very rough approximation you can use the math in this answer https://space.stackexchange.com/a/26037/12102

[didimelli]

Hello,

Thank you for your answer.
I get what you are saying and I understand the issue. However, I have the tracking of Mars' L2 as an input coming from outside my scope, so I think that the space.stackexchange answer whould be fine for me.

My problem was to translate the mathematical solution from that link to Skyfield Objects, just to have L2 behave as the bodies coming from load("de421.bsp").

[brandon-rhodes]

My to-do list tells me that I should soon expand the Skyfield documentation so that it explains how to build and return your own position objects — so this is a good chance for me to try putting together some sample code! Happily, you probably don’t need a vector sum, since you can probably compute the position as a single vector; a sum is only needed when combining multiple vectors, as is necessary when an ephemeris gives, say, a Barycenter→EarthMoon vector that then needs to be combined with its EarthMoon→Earth vector.

To keep this example simple, I won’t do any fancy Lagrangian-point math. Instead‚ I’ll simply take the Earth’s position and add a minus sign, turning it into the planet Gor which, I’ve heard, was in a science fiction book as a planet riding secretly on the other side of the Sun from us.

from skyfield import api
from skyfield.positionlib import Barycentric
from skyfield.vectorlib import VectorFunction

class Gor(VectorFunction):
    def __init__(self, ephemeris):
        self.earth = ephemeris['earth']

    def at(self, t):
        p = self.earth.at(t)
        x, y, z = p.position.au
        xdot, ydot, zdot = p.velocity.au_per_d
        b = Barycentric(
            [-x, -y, -z],           # other side of Sun
            [-xdot, -ydot, -zdot],  # opposite direction
            t=t,                    # time
            target=self,            # Gor is the target body
        )
        return b

ts = api.load.timescale()
t = ts.utc(2021, 4, 22)
eph = api.load('de421.bsp')

e = eph['earth'].at(t)
print('Earth is at', e.position.au)

g = Gor(eph)
p = g.at(t)
print('The planet Gor is at', p.position.au)

The essence here is that you need some kind of callable that given a time t returns a position object, which in this case is a Barycentric position since it is measured from the center of the Solar System. There are several ways to implement a callable in Python — a function; a closure; or a class. I’ve chosen to use a class for this example, since it provides a simple way to keep extra information around for use when at() is called — in this case I need to get to the earth object again each time at() is called, just as you will probably need to store away a mars object that you’ll use each time your at() is called.

One complication is that if t is a time object with an array of times inside, then the position returned needs to be a whole array of positions. This example makes it look easy, because -x applies equally well to either a scalar x or to a whole NumPy array x. But if you are using formulae for your L2 point, then you might run into the snag that your formula works fine for an individual time, but then gives you a NumPy broadcasting error or something when it’s given a time array, because your formula does not quite follow all the NumPy rules for array formulae. I can imagine two possible solutions:

1. The fast efficient best-case solution is that you then debug your formulae so they work with arrays. We could try to point people to NumPy documentation or Skyfield could have its own guide to writing such formulae.
2. The slow inefficient solution would be to come up with a wrapper that let the user write their L2 formula (or whatever their vector function is) for a scalar time t, such that if an array time is passed then the wrapper simply calls the formula n times and builds the results into an array to return. That would be wildly inefficient, because the user would build n separate position objects that would then all be thrown away and replaced with a single position object that held the array x,y, and z, but it might be better than nothing for users unfamiliar with doing math with NumPy arrays.
3. Okay, there’s a third possibility. The wrapper could ask the user to just return an x,y,z value without wrapping it up in a position first, so that if the wrapper needed to call the underlying routine n times, it at least would merely get back n values instead of n full objects that would then need to be thrown away. In fact, Skyfield vector functions do something like this internally, merely to avoid repeating the code for building the position object, but I’ve been considering throwing that mechanism away because it’s indirect and makes the code a bit harder to follow.

Anyway, let me know if the example code helps you put together a solution. If you could keep a list of questions and frustrations you encounter as you try adapting it to your own case, that could considerably improve the documentation I’d like to add to Skyfield to offer this pattern to other users!

[didimelli]

Hello,
thank you very much for the super clear explanation. I tested that with a simpler case and it works properly.

I will document my process and report that in this discussion.
Would you like me to do something special with that (i.e. a pr to the docs) or is it enough to report that here?

[brandon-rhodes]

I'm glad the explanation helped! It's enough to report back here; I'll be happy to take care of moving text into the docs myself!