Question additional forces

gacarita commented 1 year ago

Hello,

I'm trying to run some simulations by adding the entire spherical harmonics simulations, but I'm encountering some problems with the program using REBOUND and REBOUNDx. I have tested my code with a simple RK45, and it works fine. However, now I'm trying to implement this using REBOUND, and I'm sure that I'm doing something wrong so far.

times = np.arange(0.0,24*60*60,30)

sim = rebound.Simulation()
date = "2023-07-24 00:00"
sim.units = ('s', 'm', 'kg')

x0 = 6478.136000000000*1e3
y0 =-2.194237229530863e-13*1e3
z0 = 1.340073885192083e-15*1e3
vx0 =-1.833551089721986e-17*1e3
vy0 = 7.844114000000000*1e3
vz0 = 6.845279345644194e-12*1e3

sim.add("399",date=date)
#sim.particles[0].m = 5.972346324997922e+24
sim.status()
sim.move_to_com()
sim.add(m=0.0,x=x0,y=y0,z=z0,vx=vx0,vy=vy0,vz=vz0)
ps = sim.particles
sim.integrator = 'BS'
#sim.ri_bs.eps_rel = 1e-12
#sim.ri_bs.eps_abs = 1e-12
x = np.zeros((2,Noutputs))
y = np.zeros((2,Noutputs))
z = np.zeros((2,Noutputs))
t0 = []
w_earth = 7.292115e-5
t=sim.t

def harm(reb_sim, rebx_force, particles, N):
    sim = reb_sim.contents
    vec_r = [particles[1].x,particles[1].y,particles[1].z] # meters
    gst = w_earth * sim.t
    rotat_vect = mxv(M1, vec_r)

    ACG = [0, 0, 0]
    ACG = egm2008_kuga(radius[0], gm[0], c, s, nm, gst, rotat_vect)
    #print(ACG,sim.t)
    particles[1].ax += ACG[0]
    particles[1].ay += ACG[1]
    particles[1].az += ACG[2]

rebx = reboundx.Extras(sim)
myforce = rebx.create_force("harmonics")

myforce.force_type = "pos"
myforce.update_accelerations = harm
rebx.add_force(myforce)

#a_s = np.zeros(len(times))
for i,time in enumerate(times):
    sim.integrate(time)
    t0.append(sim.t)
    for j in range(0,2):
        x[j][i] = ps[j].x
        y[j][i] = ps[j].y
        z[j][i] = ps[j].z

I've integrated this orbit around Earth for 1 day. I suspect that the acceleration is updating in every single step of the integrator, but I only want it to update in the successful steps.

In the figure, the blue orbit is what I'm expecting, and the orange is the rebound one.

test

I appreciate any help.

hannorein commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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.

gacarita commented 1 year ago

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 commented 1 year ago

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

dtamayo / reboundx

Question additional forces #109