irfpy.util.ballistic

Ballistic trajectory related module

Associated link:

• irfpy.util.keplernumeric

class irfpy.util.ballistic.Trajectory(r0, th0, vr0, om0, mass, g=6.674e-11, m=1)

Bases: object

Initialize the trajectory.

Parameters:

• r0 – Initial position, $r_0$.

• th0 – Initial phase, ${\theta }_0$. Radians.

• vr0 – Initial radial velocity, $v_0=dr/dt(t=0)$.

• om0 – Initial angular veloicty, ${\omega }_0=d\theta /dt(t=0)$.

• mass – Mass of the central body (or effective mass of the center of gravity).

• g – Gravity of constant, $G$. Default value is for MKSA unit value.

• m – Mass of the moving object. It DOES NOT concern the trajectory, but some physical quantities do.

Theory: From the initial conditions, the trajectory is expressed in the parametric equation $r=LL/(1+e\cos (\theta +\mathrm{\Phi })$. Here, r is the distance from the central body, e is eccentricity, Phi represents the phase shift.

print_values()

angular_momentum()

Return the angular momentum.

The angular momentum is constant. It is defined by the initial condition: L=m r0^2 omega0.

$$L=m{r}_0^2{\omega }_0$$

eccentricity()

r(th)

Return the distance r as a function of th

Parameters:

th – The phase angle in radians.

Returns:

The radius, r.

Theory: From LL, e, and Phi, the radial distance is

$LL/(1+e\cos (\theta +\mathrm{\Phi }))$.

omega(th)

Return $omega$ as a function of th.

Parameters:

th – The phase angle in radians.

Returns:

The angular speed, omega.

Theory: The angular momentum (L see angular_momentum()) is conserved, so that the angular speed omega can be retrieved via the equation

$$L=mr(\theta {)}^2\omega (\theta ).$$

The distance r is retrieved from r()

v(th)

Return the radial velocity as a function of th.

Deprecated. Use vr().

vr(th)

Return the radial velocity as a function of th

Parameters:

th – The phase angle in radians.

Returns:

The radial speed.

Theory: We derived the vr as follows.

$$r=\frac{LL}{1+e\cos (\theta +\mathrm{\Phi })}.$$

By taking time derivative, we obtain

$$v_r=\frac{dr}{dt}=e\sin (\theta +\mathrm{\Phi })\frac{LL}{(1+e\cos (\theta +\mathrm{\Phi }){)}^2}\frac{d\theta }{dt}$$

By simplifying the equation, we obtain

$$v_r=\frac{el}{mLL}\sin \theta +\mathrm{\Phi }.$$

integ_distance_num(end_phase, ndiv=36000)

Integrate the distance flown from initial phase until the given end_phase.

The distance along the trajectory is integrated along the orbit. Simple daikei-kousiki (?what is the English?) is used.

Parameters:

• end_phase – The end of the phase (${\theta }_1$). Radians.

• ndiv – Division used for the separation.

Returns:

The distance.

Sample

Let’s see a lunar orbiter. The 100 km circular orbit (distance of 1838 km) is considered. It is easier to use the helper function speed_from_distances() to get the velocity for such orbiter.

>>> vper, vapo = speeds_from_distances(1838e3, 1838e3, 7.342e22)
>>> print('{:.1f} m/s'.format(vper))
1632.8 m/s

Create the Trajectory object.

>>> orbiter = Trajectory(1838e3, 0, 0, 1632.8 / 1838e3, 7.342e22)
>>> print('{:.2f}'.format(orbiter._e))   # Eccentrycity. This is (almost) circle.
0.00

Calculate the full length of the trajectory. It is $2\pi R$ where R is the radius of the orbiter (1838 km). Thus, 11548.5 km.

>>> trajectory_length = orbiter.integ_distance_num(2 * np.pi)
>>> print('{:.3e} m'.format(trajectory_length))
1.155e+07 m

integ_time_num(end_phase, ndiv=36000)

Calcuate the time from the initial time to the end at the end_phase.

The calculation uses numerical calculation

$$t={\int }_{{\theta }_0}^{{\theta }_1}\frac{d\theta }{\omega }$$

Parameters:

• end_phase – The end of the phase (${\theta }_1$). Radians.

• ndiv – Division used for the separation.

Returns:

The distance.

Sample

Let’s see a lunar orbiter. The 100 km circular orbit (distance of 1838 km) is considered. Create the Trajectory object, in the same manner as integ_distance_num().

>>> orbiter = Trajectory(1838e3, 0, 0, 1632.8 / 1838e3, 7.342e22)
>>> print('{:.2f}'.format(orbiter._e))   # Eccentrycity. This is (almost) circle.
0.00

Calculate the circular period. In this case, the total distance flying over is 11548.5 km (see integ_distance_num()). The speed of the orbiter is 1632.8 m/s. This resuts in 7072.8 s.

>>> period = orbiter.integ_time_num(2 * np.pi)
>>> print('{:.0f} s'.format(period))
7073 s

irfpy.util.ballistic.speeds_from_distances(dperi, dapo, mass, g=6.674e-11)

Calculate the speeds at peri- and apo-centers from the distances

It is common request to calulate the trajectory for, say, spacecraft with specification of the heights.

Parameters:

• dperi – Pericenter distance from the center

• dapo – Apocenter distance from the center

• mass – Mass of the central body

• g – Constant of gravity, $G$

Returns:

(vper, vapo), where vper is the pericenter spped and vapo for apocenter.

The formulation is simple: Using the energy and angular momentum conversations. Derivation is based on MKSA system.

$$\begin{array}{r}{v}_p=\frac{1}{d_p}\sqrt{\frac{2GM}{\frac{1}{d_a}+\frac{1}{d_p}}}\\ v_a=\frac{1}{d_a}\sqrt{\frac{2GM}{\frac{1}{d_a}+\frac{1}{d_p}}}\end{array}$$

Example below is an orbiter around the Moon. The pericenter altitude 30 km and apocenter altitude of 270 km are considered. Moon radius of 1738 km is taken.

>>> dp = (30 + 1738) * 1000    # in meter
>>> da = (270 + 1738) * 1000   # in meter
>>> from irfpy.util.planets import moon
>>> v_p, v_a = speeds_from_distances(dp, da, moon['mass'])
>>> print('{:.3f}'.format(v_p))
1717.547
>>> print('{:.3f}'.format(v_a))
1512.262

The results say the pericenter velocity of 1.718 km/s and apocenter velocity of 1.512 km/s

class irfpy.util.ballistic.Trajectory3D(r0vec, v0vec, mass, g=6.674e-11, m=1)

Bases: object

3-D trajectory of ballistic motion.

Theory

The class represents the ballistic motion of a mass m under a huge gravity star of mass M. (M >> m is assumed, but theoretically, if you use “equivalent mass” then one can extend this to more generic problems)

Consider the mass M at origin, and stationary. The mass m moves around under Kepler motion. Let’s take the initial position $\overrightarrow{r_0}$ with the velocity $\overrightarrow{v_0}$. $r_0=|\overrightarrow{r_0}|$ and $v_0=|\overrightarrow{v_0}|$ should not be zero. Also they should not be parallel, i.e. $\overrightarrow{r_0}\times \overrightarrow{v_0}\ne 0$.

This is a 2-D problem, which can be represented by Trajectory. This, a conversion between the 3-D coordinates to the 2-D in Trajectory is the key to solve.

Conversion:

Let’s take the new axis as follows.

• $\hat{s}$: Along the positon vector, i.e., $\hat{s}=\frac{\overrightarrow{r_0}}{r_0}$.

• $\hat{t}$: $\hat{s}-\hat{t}$ plane to include the velocity vector, with $v_t>0$. This means that $\hat{t}//(\overrightarrow{r_0}\times \overrightarrow{v_0})\times \overrightarrow{r_0}$.

Usage Let’s try the following example.

• An orbiter around the Moon (mass=7.342e22 kg). In inertia frame

• The position of orbiter at t=0 is 1850 km at the equator, longitude at 45.

>>> r0 = 1850e3 * np.array([np.cos(np.pi / 4), np.sin(np.pi / 4), 0])
>>> print(r0)
[1308147.54519511  1308147.54519511        0.        ]

• The initial velocity of the orbiter at t=0 is (0, 0.01, 1.65) km/s.

>>> v0 = 1000 * np.array([0, 0.01, 1.65])
>>> print(v0)
[    0.    10.  1650.]

Now it is time to make the Trajectory3D object.

>>> traj = Trajectory3D(r0, v0, 7.342e22)
>>> print(traj)
Trajectory3D: s=(0.707,0.707,0) t=(-0.00303,0.00303,1) r0=1.85e+06 vs0=7.07 w0=0.000892

The initial position, r0, should be converted to the original position in xyz.

>>> print(traj.st2xyz(traj.s(0), traj.t(0)))
[1308147.54519511  1308147.54519511        0.        ]

The initial velocity, v0, should be converted to the original velocity in xyz. >>> v0recal = traj.st2xyz(traj.vs(0), traj.vt(0)) >>> print(‘{v[0]:.3f} {v[1]:.3f} {v[2]:.3f}’.format(v=v0recal)) 0.000 10.000 1650.000

Later, the velocity can be retrieved by vxyz method. >>> vxyz1 = traj.vxyz(1) # Velocity at phase = 1 radian >>> print(‘{v[0]:.3f} {v[1]:.3f} {v[2]:.3f}’.format(v=vxyz1)) # It is in m/s -952.900 -947.372 912.080

>>> traj.print_values()
--- 2d ---
r0: 1850000.0
th0: 0.0
v0: 7.07106781187
om0: 0.00086487331077
...

The period of such orbit is ~7450 sec. >>> print(‘{:.3f}’.format(traj.integ_time_num(np.pi * 2))) 7452.076

Notes

All the quantities are derived via the parameter $\theta$. This is frequently a difficulty, because it is more natural to use the time t as a parameter. However, the theory of two-body collision is built such that. I do not want to dig into more in details now. Of course practically, it could be worthwhile to get a function such like $\theta =\theta (t)$ by inverting $t=t(\theta )$.

Initialize the trajectory.

Parameters:

• r0vec – Initial position, $\overrightarrow{r_0}$.

• v0vec – Initial radial velocity, $\overrightarrow{v_0}$.

• mass – Mass of the central body (or effective mass of the center of gravity).

• g – Gravity of constant, $G$. Default value is for MKSA unit value.

• m – Mass of the moving object. It DOES NOT concern the trajectory, but some deriving physical quantities do.

st2xyz(s, t)

r(th)

Return r from the phase angle as a parameter.

omega(th)

Return the tangential angular velocity

s(th)

Return the s value

t(th)

Return the t value

v(th)

Return the radial velocity

vr(th)

Return the radial velociyt. Same as v()

vth(th)

Return the phase velocity.

vs(th)

Return the s-direction velocity

vt(th)

Return the tangential velocity

xyz(th)

Return the (x, y, z) coordinates

longitude(th)

Return the longitude in degrees at the phase angle $\theta$.

Parameters:

th – Phase angle $\theta$

Returns:

The longitude in degrees.

latitude(th)

Return the Latitude in degrees at the phase angle $\theta$.

Parameters:

th – The phase angle $\theta$

Returns:

The latitude in degrees

vxyz(th)

Return (vx, vy, vz)

integ_time_num(end_phase, ndiv=36000)

integ_distance_num(end_phase, ndiv=36000)

print_values()