# irfpy.util.ballistic¶

Ballistic trajectory related module

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

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 + \Phi)$$. Here, r is the distance from the central body, e is eccentricity, Phi represents the phase shift.

print_values()[source]
angular_momentum()[source]

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()[source]
r(th)[source]

Return the distance r as a function of th

Parameters:

th – The phase angle in radians.

Returns:

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

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

omega(th)[source]

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 = m r(\theta)^2 \omega(\theta).$

The distance r is retrieved from r()

v(th)[source]

Return the radial velocity as a function of th.

Deprecated. Use vr().

vr(th)[source]

Return the radial velocity as a function of th

Parameters:

th – The phase angle in radians.

Returns:

Theory: We derived the vr as follows.

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

By taking time derivative, we obtain

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

By simplifying the equation, we obtain

$v_r = \frac{e l}{m LL}\sin{\theta + \Phi}.$
integ_distance_num(end_phase, ndiv=36000)[source]

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)[source]

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)[source]

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{split}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{split}$

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)[source]

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 $$\vec{r_0}$$ with the velocity $$\vec{v_0}$$. $$r_0 = |\vec{r_0}|$$ and $$v_0 = |\vec{v_0}|$$ should not be zero. Also they should not be parallel, i.e. $$\vec{r_0}\times\vec{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{\vec{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} // (\vec{r_0}\times\vec{v_0})\times\vec{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, $$\vec{r_0}$$.

• v0vec – Initial radial velocity, $$\vec{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)[source]
r(th)[source]

Return r from the phase angle as a parameter.

omega(th)[source]

Return the tangential angular velocity

s(th)[source]

Return the s value

t(th)[source]

Return the t value

v(th)[source]

vr(th)[source]

Return the radial velociyt. Same as v()

vth(th)[source]

Return the phase velocity.

vs(th)[source]

Return the s-direction velocity

vt(th)[source]

Return the tangential velocity

xyz(th)[source]

Return the (x, y, z) coordinates

longitude(th)[source]

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

Parameters:

th – Phase angle $$\theta$$

Returns:

The longitude in degrees.

latitude(th)[source]

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

Parameters:

th – The phase angle $$\theta$$

Returns:

The latitude in degrees

vxyz(th)[source]

Return (vx, vy, vz)

integ_time_num(end_phase, ndiv=36000)[source]
integ_distance_num(end_phase, ndiv=36000)[source]
print_values()[source]