2-body LPCs @ perihelion

[Robert Jedicke]

Aloha.

I'm running the 2-body  ephemeris generation below trying to calculate the heliocentric distance of an LPC object at the time or perihelion.  I expected it to return a heliocentric distance = perihelion distance but it does not.  Am I doing or interpreting something incorrectly?  Thanks in advance, Rob.

import pyoorb as oo

import numpy as np

oo.pyoorb.oorb_init()

tperi=54795.821223409126

q = 0.0008171432801687972

orbit = np.array([[  1100582,      # id

                           q,                      # q

                           0.9999788457712048,                # e

    np.deg2rad(           69.63497042327923 ),      # incl

    np.deg2rad(          206.2620981971092 ),     # longnode

    np.deg2rad(           48.29342159640455 ),  # argper

                       tperi,                  # tperi

                           2,                              # orbit type

                       56368.0,                  # epoch

                           1,                              # timescale type

                           0.0,                     # absolute magnitude

                           0.15]],                         # slope parameter

                     dtype=np.double, 

                     order='F' )

print( 'orbit = ', orbit )

mjd = np.arange( tperi, tperi+1, 1, dtype=np.double )

print( 'mjd = ', mjd )

epoch = np.array( list(zip(mjd, [1]*len(mjd))), dtype=np.double, order='F')

print( 'epoch = ', epoch )

eph, err = oo.pyoorb.oorb_ephemeris_full( in_orbits = orbit,

                                         in_obscode = '500',

                                     in_date_ephems = epoch,

                                        in_dynmodel = '2')

print( 'err = ', err ) 

r = eph[0][0][7]

print( "oorb ephemeris tperi, r = ", tperi, r, "\n" )

print( "      expected tperi, q = ", tperi, q, "\n" )

[Bryce Bolin]

Hello Rob,

Here is what I get when I run your code snippet:

orbit = [[1.10058200e+06 8.17143280e-04 9.99978846e-01 1.21535951e+00
3.59995274e+00 8.42879214e-01 5.47958212e+04 2.00000000e+00
5.63680000e+04 1.00000000e+00 0.00000000e+00 1.50000000e-01]]
mjd = [54795.82122341]
epoch = [[5.47958212e+04 1.00000000e+00]]
err = 0
oorb ephemeris tperi, r = 54795.821223409126 0.002625574751974971

expected tperi, q = 54795.821223409126 0.0008171432801687972

Bryce

> --
> You received this message because you are subscribed to the Google Groups "OOrb" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to oorb+uns...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/oorb/bff2627c-1e80-4ca9-afce-50cb81dfcd17n%40googlegroups.com.

[Rob Jedicke]

Thanks Bryce.

Yes, that's exactly what I get.

Shouldn't r=q at t=tperi?

[Marc Buie]

are you sure it isn't a helio- vs. bary-center issue?

To view this discussion on the web visit https://groups.google.com/d/msgid/oorb/60800127-dcc8-4d89-a005-1c00d96de599n%40googlegroups.com.

[Rob Jedicke]

Thanks for thinking about it Marc.

I thought the same thing and double-checked the documentation before making my post.

Both the pyoorb docs at https://github.com/oorb/oorb/tree/master/python

and the oorb docs at https://github.com/oorb/oorb

specify the output as heliocentric.

I have not dug into the code...

I have checked that   r-q  at the time of perihelion gets worse as q decreases.

[Bryce Bolin]

Hello Rob,

You have already brought this up, but the difference in q and r is
less than 5 digits of precision for q>0.1 and gets larger at smaller
q. Have you tried adjusting the standard gravitational parameter in
the force model to see how it affects the size of the discrepancy?

Cheers,

Bryce

> To view this discussion on the web visit https://groups.google.com/d/msgid/oorb/d003715d-02bd-4512-9317-54dd04715cffn%40googlegroups.com.

[Man-To Hui]

Perhaps due to different time systems? Like, one is TDB and the other is UT. I didn't dive into the code so don't know whether or not this is the cause. But I've definitely seen people mistaken with time systems before.

-Man-To

[Rob Jedicke]

Thanks Man-To.

I also tried using all the different time systems and, as I recall, got essentially the same answers.

[Rob Jedicke]

Thanks Bryce.

Fiddling with G would make me feel like I'm playing God :) or, at the very least, pretending to be a cosmologist.

I thought it might have something to do with solar moments (i.e. Jn) but, as I understand it, that should not affect the 2-body calculation.

Rob

[Marc Buie]

ok, another stupid question, perhaps.  Is GR handled properly in oorb?

To view this discussion on the web visit https://groups.google.com/d/msgid/oorb/8251f0d5-5236-42e0-add4-f06cba300fe4n%40googlegroups.com.

[Mikael Granvik]

First, r=q at t=tperi, so clearly something is wrong. The discrepancy is too
large for a timescale issue. GR doesn't matter for the 2-body approximation.
This may be a bug/feature in the code solving the Kepler problem. Similar issues
vanished when I switched to the new universal Kepler solver by Wisdom &
Hernandez, and I need to make sure that the old code isn't causing the problems
we see here.

Mikael

> https://groups.google.com/d/msgid/oorb/04c60c7b-bccc-954e-2397-e2ece16b72f0%40gmail.com.
>
>

[Mikael Granvik]

Okay, I have some answers... But first something I forgot to spell out
explicitly in my previous message: I'm *very* happy to see that the OpenOrb list
has active members that are willing to help out with issues such as the present
one. Thank you all for taking the time from your busy schedules to contribute to
the discussion!

Now, let's proceed to the answers. The short take home message is that there is
*no* error in the ephemeris that you compute.

The longer explanation is that you have to account for the fact that the
ephemeris provides you the situation as apparent from the location and date for
which you compute the ephemeris. Since the speed of light is finite and it takes
about 8 minutes for photons to get from the Sun to the Earth, an ephemeris
computed for the geocenter for a date equal to the time of perihelion will
actually show the coordinates some 8 minutes before the object reaches
perihelion. Although the 8-min difference may sound negligible (I was fooled by
it!), it turns out that the speed of the object at perihelion (q=0.0008...) is
so large, about 0.85 au/d = 1500 km/s based on elementary celestial mechanics,
that it explains the difference. (Note that this is a hypothetical situation
because with q<0.00465 an object would collide with the Sun.)

The silver lining with this thread is that I found another bug in the code and
it is the reason it took me a bit longer to figure out what was happening. I
started hunting for the bug by resorting to the command-line version of OpenOrb,
that is, the one with a Fortran front end. It has a configuration option that
allows you to determine which types of orbital elements are used in the
computations. Playing around with them I found out that by feeding in cometary
elements (as provided by RJ in the message that started this thread) and forcing
the code to calculate with cometary elements I could compute an ephemeris for
which r=q when t=tperi. So I figured that the use of the other element types
must lead to erroneous results, and started to hunt for it there. Only yesterday
did I realize that I had made a mistake in my thinking, and it's instead the use
of cometary elements that produce the wrong but correctly-looking result. I've
yet to locate and fix the bug, but it doesn't affect the Python version because
the incoming orbital elements are internally converted to cartesian elements.

Mikael

[Rob Jedicke]

Mahalo Mikael for looking into this issue and identifying a bug, even if it wasn't the bug we thought we reported.

Your explanation for the 'original' bug makes perfect sense, that it is the light travel time for the observer.  I am embarrassed that I didn't realize it myself.

So is there a way to use OpenOrb to determine the 'actual' location of an object like JPL Horizon's VECTOR ephemeris?

FWIW, I too am appreciative of the OpenOrb community of users. Bryce Bolin has already helped me address a few trivial use issues without having to resort to this users list.

Good luck on resolving the 'new' bug!

Rob

[Mikael Granvik]

Well, if you're embarrassed then I need to be doubly so because I should know
what the code is doing. But as I tried to convey, this example case is extreme
compared to the orbits that we're accustomed to and our gut feelings don't work.

To get something similar to the VECTOR type in JPL you just need to propagate
the orbits to your desired date and output them in the cartesian format. Not
sure if the Python bindings allow for this, but at least the oorb binary can do
it.

Mikael

> To view this discussion on the web visit https://groups.google.com/d/msgid/oorb/8874daf7-1ae7-40bf-a78c-96a0f3f3e525n%40googlegroups.com.
>
>

[Lynne Jones]

The python bindings should indeed do this -- oorb.oorb_propagation and oorb.oorb_element_transformation 

To view this discussion on the web visit https://groups.google.com/d/msgid/oorb/alpine.OSX.2.23.453.2009252134480.70913%40lm5-flafo-03.

--

=========================
Dr R. Lynne Jones
Research Scientist @
University of Washington / Dirac Institute / Rubin Observatory
she/hers
=========================

[Lynne Jones]

Hi Mikael, 

As a follow-on question here: is there an assumption somewhere inside OpenOrb that the timescales for orbital elements are actually TT, instead of UTC as Rob has here? I notice that doing this problem in python, I had two issues: 

* one was that the position (calculated using oorb_element_transformation and then oorb_propagation) matched the expected perihelion distance when using TT for the epoch timescales, but it did not match when using UTC for the epoch timescales

* the second was that doing the process in the reverse order (oorb_propagation and then oorb_element_transformation) led to an error (err=58) and very inaccurate results .. it seems like you can only feed cartesian orbits into oorb_propagation? I know this used to be the case, but then I thought that problem had gone away with some work that happened about the time Michael extended the pyoorb.f90 file .. and now it seems to be back. 

Any thoughts on these issues? 

Lynne