Note
Go to the end to download the full example code.
Binary orbit functions#
This example shows how to use binary orbit functions
8 import numpy as np
9
10 import shamrock
Define the unit system
14 si = shamrock.UnitSystem()
15 sicte = shamrock.Constants(si)
16 codeu = shamrock.UnitSystem(
17 unit_time=sicte.year(),
18 unit_length=sicte.au(),
19 unit_mass=sicte.sol_mass(),
20 )
21 ucte = shamrock.Constants(codeu)
22 G = ucte.G()
Utility to plot the resulting orbits
27 def plot_orbits(m1, m2, a, e, roll, pitch, yaw):
28 """
29 Plot the orbit of a binary system by varying the true anomaly
30 """
31 import matplotlib.pyplot as plt
32
33 ax = plt.figure().add_subplot(projection="3d")
34
35 x1, x2, y1, y2, z1, z2 = [], [], [], [], [], []
36
37 vx1, vx2, vy1, vy2, vz1, vz2 = [], [], [], [], [], []
38
39 max_nu = np.pi
40 min_nu = -np.pi
41 if e >= 1: # if parabolic do not exceed pi
42 max_nu = 0.75 * np.pi
43 min_nu = -0.75 * np.pi
44
45 for nu in np.linspace(min_nu, max_nu, 200, endpoint=False):
46 # To see the orbit start
47 if nu > 1.8 * np.pi:
48 break
49
50 _x1, _x2, _v1, _v2 = shamrock.phys.get_binary_rotated(
51 m1=m1, m2=m2, a=a, e=e, nu=nu, G=G, roll=roll, pitch=pitch, yaw=yaw
52 )
53
54 # print(_x1, _x2, _v1, _v2)
55 x1.append(_x1[0])
56 x2.append(_x2[0])
57 y1.append(_x1[1])
58 y2.append(_x2[1])
59 z1.append(_x1[2])
60 z2.append(_x2[2])
61
62 vx1.append(_v1[0])
63 vx2.append(_v2[0])
64 vy1.append(_v1[1])
65 vy2.append(_v2[1])
66 vz1.append(_v1[2])
67 vz2.append(_v2[2])
68
69 ax.plot(x1, y1, z1, "-o", markevery=[0, 50, 100, 150])
70 ax.plot(x2, y2, z2, "-o", markevery=[0, 50, 100, 150])
71
72 for i in range(0, len(x1), 50):
73 vnorm = np.sqrt(vx1[i] ** 2 + vy1[i] ** 2 + vz1[i] ** 2) * 0.03
74 ax.quiver(x1[i], y1[i], z1[i], vx1[i], vy1[i], vz1[i], color="r", length=vnorm)
75 for i in range(0, len(x2), 50):
76 vnorm = np.sqrt(vx1[i] ** 2 + vy1[i] ** 2 + vz1[i] ** 2) * 0.03
77 ax.quiver(x2[i], y2[i], z2[i], vx2[i], vy2[i], vz2[i], color="b", length=vnorm)
78
79 ax.set_aspect("equal")
80 ax.set_xlabel("x")
81 ax.set_ylabel("y")
82 plt.show()
Orbit 1
87 plot_orbits(0.7, 0.3, 1.0, 0.3, 0.0, 0.0, 0.0)
Orbit 2
91 plot_orbits(0.5, 0.5, 1.0, 0.3, 1.0, 0.0, 0.0)
Orbit 3
95 plot_orbits(0.5, 0.5, 1.0, 0.0, 1.0, 0.0, 0.0)
Orbit 4
99 plot_orbits(0.9, 0.1, 1.0, 0.0, 0.0, 1.0, 0.0)
Orbit 5 (hyperbolic)
103 plot_orbits(0.9, 0.1, 1.0, 1.2, 0.0, 1.0, 0.0)
Total running time of the script: (0 minutes 1.687 seconds)
Estimated memory usage: 159 MB