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 + \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+\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 = m r(\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+\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)

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{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)

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)

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()