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
Use shamrock documentation style for matplotlib
14 shamrock.matplotlib.set_shamrock_mpl_style()
Define the unit system
19 si = shamrock.UnitSystem()
20 sicte = shamrock.Constants(si)
21 codeu = shamrock.UnitSystem(
22 unit_time=sicte.year(),
23 unit_length=sicte.au(),
24 unit_mass=sicte.sol_mass(),
25 )
26 ucte = shamrock.Constants(codeu)
27 G = ucte.G()
Utility to plot the resulting orbits
32 def plot_orbits(m1, m2, a, e, roll, pitch, yaw):
33 """
34 Plot the orbit of a binary system by varying the true anomaly
35 """
36 import matplotlib.pyplot as plt
37
38 ax = plt.figure().add_subplot(projection="3d")
39 ax.minorticks_off()
40
41 x1, x2, y1, y2, z1, z2 = [], [], [], [], [], []
42
43 vx1, vx2, vy1, vy2, vz1, vz2 = [], [], [], [], [], []
44
45 max_nu = np.pi
46 min_nu = -np.pi
47 if e >= 1: # if parabolic do not exceed pi
48 max_nu = 0.75 * np.pi
49 min_nu = -0.75 * np.pi
50
51 for nu in np.linspace(min_nu, max_nu, 200, endpoint=False):
52 # To see the orbit start
53 if nu > 1.8 * np.pi:
54 break
55
56 _x1, _x2, _v1, _v2 = shamrock.phys.get_binary_rotated(
57 m1=m1, m2=m2, a=a, e=e, nu=nu, G=G, roll=roll, pitch=pitch, yaw=yaw
58 )
59
60 # print(_x1, _x2, _v1, _v2)
61 x1.append(_x1[0])
62 x2.append(_x2[0])
63 y1.append(_x1[1])
64 y2.append(_x2[1])
65 z1.append(_x1[2])
66 z2.append(_x2[2])
67
68 vx1.append(_v1[0])
69 vx2.append(_v2[0])
70 vy1.append(_v1[1])
71 vy2.append(_v2[1])
72 vz1.append(_v1[2])
73 vz2.append(_v2[2])
74
75 ax.plot(x1, y1, z1, "-o", markevery=[0, 50, 100, 150])
76 ax.plot(x2, y2, z2, "-o", markevery=[0, 50, 100, 150])
77
78 for i in range(0, len(x1), 50):
79 vnorm = np.sqrt(vx1[i] ** 2 + vy1[i] ** 2 + vz1[i] ** 2) * 0.03
80 ax.quiver(x1[i], y1[i], z1[i], vx1[i], vy1[i], vz1[i], color="r", length=vnorm)
81 for i in range(0, len(x2), 50):
82 vnorm = np.sqrt(vx1[i] ** 2 + vy1[i] ** 2 + vz1[i] ** 2) * 0.03
83 ax.quiver(x2[i], y2[i], z2[i], vx2[i], vy2[i], vz2[i], color="b", length=vnorm)
84
85 ax.set_aspect("equal")
86 ax.set_xlabel("x")
87 ax.set_ylabel("y")
88 plt.show()
Orbit 1
93 plot_orbits(0.7, 0.3, 1.0, 0.3, 0.0, 0.0, 0.0)
Orbit 2
97 plot_orbits(0.5, 0.5, 1.0, 0.3, 1.0, 0.0, 0.0)
Orbit 3
101 plot_orbits(0.5, 0.5, 1.0, 0.0, 1.0, 0.0, 0.0)
Orbit 4
105 plot_orbits(0.9, 0.1, 1.0, 0.0, 0.0, 1.0, 0.0)
Orbit 5 (hyperbolic)
109 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.825 seconds)
Estimated memory usage: 159 MB