Question additional forces about reboundx

[gacarita]

Yes, it's a little bit complicated.

Actually, I have done a few tests with a colleague (who has written the force function @safwanaljbaae) comparing both accelerations values and everything looks good.

I believe that could be a reference frame problem because even the non-perturbed orbit of REBOUND has 200m of error in comparison with the GMAT.

I agree with you, I don't believe that the problem is on REBOUND or even REBOUNDX, I know REBOUND and I'm using it for a while in astronomy problems but never in astrodynamics. I thought it was something that I was doing wrong with a code line.

I'll keep searching for the problem ...

If it's possible I strongly suggest the implementation of the feature of the gravity potential perturbation (not only in the Jn0 but all harmonics [Cnm and Snm]). This force function could be a starting point for the implementation. That could lead a larger number of users in astrodynamics.

I'll close this at the moment if I have any news I'll write here.

[hannorein]

It's hard to help you without seeing the function egm2008_kuga which seems to do all the hard work. Have you tried just using REBOUND with an additional force but without REBOUNDx? That would reduce one extra layer of complexity from your problem.

[gacarita]

Here are the function and the file containing the necessary data to do this calculation.

I tried to use REBOUND additional forces but the same error has shown.

I have tested the algorithm using the odeint routine in python and got +- 20 meters of error in comparison with an orbit given by GMAT tool.

I'm aware that the difference between numerical integrators could lead to different orbits, but I'm expecting some similar orbit.

EGM2008_ZeroTide_Truncated100.dat.zip

import pandas as pd
import numpy as np

def egm2008_kuga(req, gm, c, s, nm, gst, x_in):
    """
    :param req: Equatorial radius
    :param gm: The mass od the central body
    :param c: c-harmonics
    :param s: s-harmonics
    :param nm: Desired degree and degree
    :param gst: Greenwich Sideral Time
    :param x_in: Cartesian coordinates
    :return: Accelerations due to the geo-potential

    Kuga, Helio & Carrara, Valdemir. (2013). Fortran-and C-codes for higher order and degree geopotential and
    derivatives computation.

    https://mathworld.wolfram.com/LegendrePolynomial.html 44 -> eq 17
    """

    nm += 1
    cang = np.cos(gst)
    sang = np.sin(gst)

    x = [x_in[0] * cang + x_in[1] * sang, -x_in[0] * sang + x_in[1] * cang, x_in[2]]
    r = np.sqrt(x[0] ** 2 + x[1] ** 2 + x[2] ** 2)
    q = req / r
    ct = x[2] / r  # cos(theta), theta = colatitude
    st = np.sqrt(1 - ct ** 2)  # sin(theta)
    lamb = np.arctan2(x[1], x[0])  # Longitude
    slambda = np.sin(lamb)
    clambda = np.cos(lamb)
    Gmr = gm / r

    v_lambda, v_theta, Vr = 0, 0, 0

    pn, qn = [0] * nm, [0] * nm
    pn[0] = 1
    pn[1] = np.sqrt(3) * st
    qn[0] = 1
    qn[1] = q

    for m in range(2, nm):
        pn[m] = st * np.sqrt(1 + 0.5 / m) * pn[m - 1]  # (9)
        qn[m] = q * qn[m - 1]

    sm, cm = 0, 1

    for m in range(0, nm):
        Pnm = pn[m]
        dPnm = -m * Pnm * ct / st
        Pnm1m = Pnm
        Pnm2m = 0

        qc = qn[m] * c[m][m]
        qs = qn[m] * s[m][m]
        xmc = qc * Pnm  # (20)
        xms = qs * Pnm  # (20)
        xthetamc = qc * dPnm
        xthetams = qs * dPnm
        xrmc = (m + 1) * qc * Pnm  # (22)
        xrms = (m + 1) * qs * Pnm  # (22)

        for n in range(m + 1, nm):
            anm = np.sqrt(((n * 2 - 1) * (n * 2 + 1)) / ((n - m) * (n + m)))  # (7)
            bnm = np.sqrt(((n * 2 + 1) * (n + m - 1) * (n - m - 1)) / ((n - m) * (n + m) * (n * 2 - 3)))  # (8)

            Pnm = anm * ct * Pnm1m - bnm * Pnm2m  # (5) the Legendre polynomial

            fnm = np.sqrt(((n ** 2 - m ** 2) * (n * 2 + 1)) / (n * 2 - 1))  # (18)
            dPnm = -n * Pnm * ct / st + fnm * Pnm1m / st  # (17) first derivative of the Legendre polynomial

            Pnm2m = Pnm1m
            Pnm1m = Pnm

            if n >= 2:
                qc = qn[n] * c[n][m]
                qs = qn[n] * s[n][m]
                xmc = (xmc + qc * Pnm)  # (13)
                xms = (xms + qs * Pnm)  # (14)
                xthetamc = (xthetamc + qc * dPnm)  # (20)
                xthetams = (xthetams + qs * dPnm)  # (20)
                xrmc = (xrmc + (n + 1) * qc * Pnm)  # (22) ?????????
                xrms = (xrms + (n + 1) * qs * Pnm)  # (22) ?????????

        v_lambda = v_lambda + m * (-xms * cm + xmc * sm)  # (19) ?????????
        v_theta = v_theta + (xthetamc * cm + xthetams * sm)  # (21)
        Vr = Vr + (xrmc * cm + xrms * sm)  # (23)

        cm_temp = clambda * cm - sm * slambda
        sm_temp = clambda * sm + cm * slambda

        cm = cm_temp
        sm = sm_temp

    v_lambda = -Gmr * v_lambda
    v_theta = Gmr * v_theta
    Vr = -(Gmr / r) * Vr  # (23)

    ac = [st * clambda * Vr - ct * clambda * v_theta / r - slambda * v_lambda / (st * r),
          st * slambda * Vr - ct * slambda * v_theta / r + clambda * v_lambda / (st * r),
          ct * Vr + st * v_theta / r]  # (24)

    ac_out = [ac[0] * cang - ac[1] * sang, ac[0] * sang + ac[1] * cang, ac[2]]

    return ac_out




import julian
from skyfield.api import load
from skyfield import framelib
from skyfield.functions import mxv
import rebound
import reboundx


gm = [3986004.415e+8]

radius = [6378136.3e0]

nm = 10

if nm > 1:
    file_harmonics = 'EGM2008_ZeroTide_Truncated100.dat'
    harmonics = pd.read_csv(file_harmonics, delim_whitespace=True, skiprows=0,
                            names=['n', 'm', 'c', 's', 'x1', 'x2'])
    s = [[0] * (nm + 1) for _ in range(nm + 1)]
    c = [[0] * (nm + 1) for _ in range(nm + 1)]
    for n in range(nm + 1):
        for m in range(nm + 1):
            data = harmonics[(harmonics.n == n) & (harmonics.m == m)]
            if len(data) != 0:
                if 'D' in str(list(data.s)[0]).split()[0]:
                    s[n][m] = float(list(data.s)[0].replace('D', 'E'))
                    c[n][m] = float(list(data.c)[0].replace('D', 'E'))
                else:
                    s[n][m] = list(data.s)[0]
                    c[n][m] = list(data.c)[0]
else:
    nm = 2
    s = [[0] * (nm + 1) for _ in range(nm + 1)]
    c = [[0] * (nm + 1) for _ in range(nm + 1)]

day0 = 24
year0 =  2023 
month0 = 7
hour0 = 0.0
minute0 = int(0.0)
second0 = int(0.0)
ts = load.timescale()
M1 = framelib.itrs.rotation_at(ts.utc(year0, month0, day0,hour0, minute0, second0))

[hannorein]

Oh wow that is a complicated force function. I'm afraid I can't check that in detail for you. I don't think there is an issue with REBOUND. Also, the way you're interfacing to REBOUND seems good.

Here are some suggestions on what you might try next:

• Print out the accelerations and compare them to the code where you get the right results.
• Check that all your units are correct. Some of them such as sidereal time could be difficult to get right. Again, a direct comparison to the other code might give you important hints.

[hannorein]

You're welcome to provide an implementation for these forces. I'm sure @dtamayo would be happy to include them in REBOUNDx.