planetary obliquity, trace rotation axis along orbit with python

Question

I'm trying to make a simple simulation of a planet that is being orbited by a moon. So far I have a 2 body problem that solves the planet and moon orbit. Now I would like to add a fixed rotation axis to the planet and see how it is affected by the moon. Any idea how this can be done by using python?

The two body problem can be run with the code below:

import pylab
import math
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

# Set Constants
G = 6.67e-11
AU = 1.5e11
daysec = 24.0*60*60

Ma =5.972e24   # Planet mass in Kg
Mb = 7.348e22   # Moon mass in Kg


gravconst = G*Ma*Mb

# Set up starting conditions

# Planet
xa = 0.0
ya = 0.0
za = 0.0

xva = 0.0
yva = 0.0
zva = 0.0

# Moon
xb = 384400000
yb = 0.0
zb = 0.0

xvb = 0.0
yvb = 1000.0
zvb = 0.0

# Time steps
t = 0.0
dt = 0.01*daysec

# Coordinate lists
xalist = []
yalist = []

xblist = []
yblist = []

zalist = []
zblist = []

# Loop
while t < 100.0*daysec:
    # Compute Force
    rx = xb-xa
    ry = yb-ya
    rz = zb-za

    modr3 = (rx**2+ry**2+rz**2)**1.5

    fx = -gravconst*rx/modr3
    fy = -gravconst*ry/modr3
    fz = -gravconst*rz/modr3

    # Update quantities
    # Moon
    xvb += fx*dt/Mb
    yvb += fy*dt/Mb
    zvb += fz*dt/Mb

    xb += xvb*dt
    yb += yvb*dt
    zb += zvb*dt

    # Planet
    xva += -fx*dt/Ma
    yva += -fy*dt/Ma
    zva += -fz*dt/Ma

    xa += xva*dt
    ya += yva*dt
    za += zva*dt

    t += dt

    # Saving them in lists
    xalist.append(xa)
    yalist.append(ya)
    zalist.append(za)

    xblist.append(xb)
    yblist.append(yb)
    zblist.append(zb)


xalist[:] = [x / 1e6 for x in xalist]
yalist[:] = [x / 1e6 for x in yalist]
zalist[:] = [x / 1e6 for x in zalist]

xblist[:] = [x / 1e6 for x in xblist]
yblist[:] = [x / 1e6 for x in yblist]
zblist[:] = [x / 1e6 for x in zblist]

#Creating the point to represent the planet at the origin (not to       scale),
plt.scatter(0,0,s=200,color='blue')
plt.annotate('Planet', xy=(-45,-50))
plt.scatter(xblist[0],0,s=100,color='grey')
plt.annotate('Mond', xy=(xblist[0]-45,-50))

# Plotting
pylab.plot(xalist, yalist, "-g")
pylab.plot(xblist, yblist, "-r")
plt.axhline(0, color='black')
plt.axvline(0, color='black')
pylab.axis("equal")
pylab.xlabel("X (Mio. Meter)")
pylab.ylabel("Y (Mio. Meter)")
pylab.show()

Answer 1

Not an answer as I am no expert in the matter but just some hints instead (was unreadable in form of comment)

The stuff you want to add is quite complicated as you would need take into account:

1. moving masses of both bodies

   so you need "contact surfaces elevation" of body that has any moving masses (like oceans, magma, rotating core etc) so you can compute the real center of gravity for each time. Also you need to apply forces to the moving mass itself too (which is not driven just by gravity and rotation but also with resonance and mainly inertia) do not forget Earth has also core and magma around it not just oceans so you need to take into account at least 3 surfaces...

2. non homogenity of mass distribution of both bodies

   so you can compute center of gravity, and the quadratic mass inertia for rotations in respect to actual rotation axis

3. planet/moon is usually at least 3 body problem not just 2

   as there is the local star too affecting the moon orbit quite a lot...

Depending on the numbers some effects will be so small that can be discarded but some are not (especially with inertia+resonance in place).

The rotation equations are similar to the position/speed/acceleration like you already got. Its called Newton D'Alembert integration/physics but you would need to implement transform matrices to do this.

See few related QAs:

• Rigid Body Physics Rotations
• realistic n-body solar system

Pay attention to the last link for accuracy improvement of your integration as right now its very bad and the orbit will be deformed no matter how small the dt due to fact that the gravity vector changes with time but right now you affect it in wrong direction (that is correct only at the start of dt interval) for each integration iteration.

As you can see that is a lot of stuff to handle and most simulation programs I saw do not do it (mine included) ... they fake it by nutation and precession constants instead.