Comments (5)
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.
from reboundx.
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.
from reboundx.
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))
from reboundx.
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.
from reboundx.
You're welcome to provide an implementation for these forces. I'm sure @dtamayo would be happy to include them in REBOUNDx.
from reboundx.
Related Issues (20)
- Problem with the example: full version of yarkovsky_effect: Python 2.7 HOT 8
- Yarkovsky effect in Reboundx : question about the 'k' coefficient HOT 2
- Signed int HOT 2
- Radiation forces blocked by particles HOT 3
- Need to attach reboundx.Extras instance to simulation before setting params. HOT 8
- Failed building editable for reboundx HOT 6
- Cannot reproduce the plot in InnerDiskEdge.ipynb HOT 1
- Trouble implementing new effect
- Extras not attatching to sim for Rebound 4.3.2 HOT 7
- "Process finished with exit code 139" when running example code HOT 4
- Error with setting gr.params in GR notebook HOT 12
- ReboundX - installing problem with Python HOT 1
- REBOUNDx trouble with satellites HOT 2
- Confusing function name in "type_I_migration.c" HOT 14
- Python QuickStart Guide HOT 2
- reboundx install failed on my Macbook HOT 6
- "OSError: librebound.cpython-37m-x86_64-linux-gnu.so: cannot open shared object file: No such file or directory" when importing reboundx on Google Colab HOT 6
- REBOUNDx error: Force 'type_I_migration' not found in REBOUNDx library. HOT 5
- Setting beta via sim.particles[i].beta does not work? HOT 2
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from reboundx.