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
47 # To see the orbit start
48 if nu > 1.8 * np.pi:
49 break
50
51 _x1, _x2, _v1, _v2 = shamrock.phys.get_binary_rotated(
52 m1=m1, m2=m2, a=a, e=e, nu=nu, G=G, roll=roll, pitch=pitch, yaw=yaw
53 )
54
55 # print(_x1, _x2, _v1, _v2)
56 x1.append(_x1[0])
57 x2.append(_x2[0])
58 y1.append(_x1[1])
59 y2.append(_x2[1])
60 z1.append(_x1[2])
61 z2.append(_x2[2])
62
63 vx1.append(_v1[0])
64 vx2.append(_v2[0])
65 vy1.append(_v1[1])
66 vy2.append(_v2[1])
67 vz1.append(_v1[2])
68 vz2.append(_v2[2])
69
70 ax.plot(x1, y1, z1, "-o", markevery=[0, 50, 100, 150])
71 ax.plot(x2, y2, z2, "-o", markevery=[0, 50, 100, 150])
72
73 for i in range(0, len(x1), 50):
74 vnorm = np.sqrt(vx1[i] ** 2 + vy1[i] ** 2 + vz1[i] ** 2) * 0.03
75 ax.quiver(x1[i], y1[i], z1[i], vx1[i], vy1[i], vz1[i], color="r", length=vnorm)
76 for i in range(0, len(x2), 50):
77 vnorm = np.sqrt(vx1[i] ** 2 + vy1[i] ** 2 + vz1[i] ** 2) * 0.03
78 ax.quiver(x2[i], y2[i], z2[i], vx2[i], vy2[i], vz2[i], color="b", length=vnorm)
79
80 ax.set_aspect("equal")
81 ax.set_xlabel("x")
82 ax.set_ylabel("y")
83 plt.show()
Orbit 1
88 plot_orbits(0.7, 0.3, 1.0, 0.3, 0.0, 0.0, 0.0)
Orbit 2
92 plot_orbits(0.5, 0.5, 1.0, 0.3, 1.0, 0.0, 0.0)
Orbit 3
96 plot_orbits(0.5, 0.5, 1.0, 0.0, 1.0, 0.0, 0.0)
Orbit 4
100 plot_orbits(0.9, 0.1, 1.0, 0.0, 0.0, 1.0, 0.0)
Orbit 5 (hyperbolic)
104 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.182 seconds)
Estimated memory usage: 65 MB