irfpy.util.ballistic
¶
Ballistic trajectory related module
Associated link:
- class irfpy.util.ballistic.Trajectory(r0, th0, v0, 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.
v0 – 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.
- omega(th)[source]¶
Return \(omega\) as a function of th.
- Parameters
th – The phase angle in radians.
- v(th)[source]¶
Return the radial velocity as a function of
th
- Parameters
th – The phase angle in radians.
- 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 asinteg_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 inTrajectory
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.
- longitude(th)[source]¶
Return the longitude in degrees at the phase angle \(\theta\).
- Parameters
th – Phase angle \(\theta\)
- Returns
The longitude in degrees.