TransWikia.com

Rotation Reference Change Transformation CRTBP Python

Physics Asked by John Shepard on February 26, 2021

I have created the following script that will propagate the dynamics of the CRTBP set for the Earth-Moon system. I have the initial conditions set for a vertical Lyaponov orbit. The goal I am trying to achieve is to map this orbit in both relative and inertial reference frames.

My attempt below is using the following equations found on the lines prior to the comment "Plotting of the data":

x_in = x_ro  cos⁡(φ(t)) - y_ro sin⁡(φ(t))
y_in = x_ro  sin⁡(φ(t)) + y_ro cos⁡(φ(t))
z_in = z_ro

Ideally, I would also like to do the conversion for the velocity as well, but I am currently struggling with the position transformation.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import librationPoints
from mpl_toolkits import mplot3d


def crtbp(y, t, mu):
    G = 1
    r1 = np.sqrt((mu+y[0])**2 + (y[1]**2) + (y[2]**2))
    r2 = np.sqrt((1-mu-y[0])**2 + (y[1]**2) + (y[2]**2))

    m1 = 1 - mu
    m2 = mu

    ydot = [y[3],
            y[4],
            y[5],
            y[0]+2*y[4]+G*m1*(-mu-y[0])/(r1**3)+G*m2*(1-mu-y[0])/(r2**3),
            y[1]-2*y[3]-G*m1*(y[1])/(r1**3)-G*m2*y[1]/(r2**3),
            -G*m1*y[2]/(r1**3)-G*m2*y[2]/(r2**3)]
    return ydot


# mass of primary bodies
m1_i = 5.9722 * 10**24
m2_i = 7.34767309 * 10**22

M = m1_i + m2_i                       # total mass of the system
mu_i = m2_i/M                         # order of unity for dimensional form

# distance from the Earth to the Moon
d = 384400  # km # unit length

# simulation end time
tau = 6.2935

# initial state vector
x_i = 1.0619
y_i = 0.0
z_i = 0.0
vx_i = 0
vy_i = -1.8721
vz_i = 0.6525

# RK45 integration for the CRTBP Setup
y0 = np.array([x_i, y_i, z_i, vx_i, vy_i, vz_i])  # initial value
print(y0)
t0 = 0                                  # simulation start time
N = 1000000                             # minimal acceptable performance N = 1000000
step = 1/N                              # step size

# RK45 integration for the CRTBP Setup
numSteps = 2000
tspan = np.linspace(t0, tau, numSteps)
solution = odeint(crtbp, y0, tspan, args=(mu_i,), rtol=1e-12, atol=1e-12)

G = 1e-9 * 6.674e-11
t = np.sqrt((d*d*d) / (G * M))
t_values = tspan

x_values = solution[:, 0]
y_values = solution[:, 1]
z_values = solution[:, 2]

phi = 1/t
nx = [(x * np.cos(phi*t)) - (y * np.sin(phi*t)) for x, t, y in zip(x_values, t_values, y_values)]
ny = [(x * np.sin(phi*t)) + (y * np.cos(phi*t)) for y, t, x in zip(y_values, t_values, x_values)]
nz = z_values

vx_values = solution[:, 3]
vy_values = solution[:, 4]
vz_values = solution[:, 5]

# Plotting of the data
# 2D Results
plt.title("Transfer Maneuver")
plt.xlabel("x (non-dimensional)")
plt.ylabel("y (non-dimensional)")
plt.plot([1-mu_i], [0], marker='o', markersize=3, color='blue')  # P2
plt.plot([-mu_i], [0], marker='o', markersize=3, color='blue')  # P1
plt.plot(x_values, y_values)
plt.plot(nx, ny)
plt.show()

# 3D Results
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot3D(x_values, y_values, z_values, 'b')
ax.plot3D(nx, ny, nz, 'r')
plt.show()

Currently when I run the code the transformation does not seem to be changing the position. Two orbits are overlapping which is not correct. In the rotational reference frame I expect to see the figure eight or hour glass shape. However in the inertial frame I am expecting to see an orbit about Earth with some inclination with respect to the Earth-Moon plane. What is it in the transformation that I am missing or not doing correct? Any help would be appreciated.

*** Update ***

When I change the relationship from 1/t to t I get the following image.
enter image description here

However is the orbit when converting into inertial frame suppose to be open or closed, is there something wrong with the equations, how I am using them, or is there some understanding that I am missing?

*** I think I solved it ***

When adjusting the model run time tau = 0.0000011 the image is corrected ending and starting at the same location. This seems to be very sensitive to the initial inputs and requires that one rotation is completed. Making this change produces the following image.
enter image description here

Then I apply this to the velocity list
enter image description here

and I get a similar image. I believe this to be solved now for both position and velocity.

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP