import numpy as np
import matplotlib.pyplot as plt
G = 6.6743e-11
c = 299792458
R_ns = 12000
m1 = m2 = 3.5e31
init_dist = 1e4
dt = 1e-8
pos1 = np.array([-init_dist / 2, 0], dtype=float)
pos2 = np.array([init_dist / 2, 0], dtype=float)
# Orbital velocity magnitude for each star
v = np.sqrt(G * (m1 + m2) / init_dist)
# Initial velocities perpendicular to separation vector (x-axis)
v1 = np.array([0, v], dtype=float) # Star 1 moves +y
v2 = np.array([0, -v], dtype=float) # Star 2 moves -y
def compute_accelerations(pos1, pos2):
r_vec = pos2 - pos1
r_mag = np.linalg.norm(r_vec)
r_hat = r_vec / r_mag
force_mag = G * m1 * m2 / r_mag ** 2
a1 = force_mag / m1 * r_hat
a2 = -force_mag / m2 * r_hat
return a1, a2
a1, a2 = compute_accelerations(pos1, pos2)
def verlet_integration(pos1, pos2, a1, a2, v1, v2):
pos1_new = pos1 + v1 * dt + 0.5 * a1 * dt ** 2
pos2_new = pos2 + v2 * dt + 0.5 * a2 * dt ** 2
a1_new, a2_new = compute_accelerations(pos1_new, pos2_new)
v1_new = v1 + 0.5 * (a1 + a1_new) * dt
v2_new = v2 + 0.5 * (a2 + a2_new) * dt
return pos1_new, pos2_new, v1_new, v2_new, a1_new, a2_new
def peters_mathews(pos1, pos2, v1, v2):
r_vec = pos2 - pos1
r = np.linalg.norm(r_vec)
dr_dt = - (64 / 5) * (G ** 3 * m1 * m2 * (m1 + m2)) / (c ** 5 * r ** 3)
delta_r = dr_dt * dt
r_new = r + delta_r # dr_dt negative, so r_new < r
if r_new < R_ns / 2:
print("Merger condition reached. Exiting.")
exit()
r_hat = r_vec / r
pos1_new = pos1 + 0.5 * (r_new - r) * r_hat
pos2_new = pos2 - 0.5 * (r_new - r) * r_hat
# Calculate new orbital speed magnitude
v_mag = np.sqrt(G * (m1 + m2) / r_new)
# Calculate tangential direction by rotating r_hat by 90 degrees CCW
tangential_dir = np.array([-r_hat[1], r_hat[0]])
v1_new = v_mag * tangential_dir
v2_new = -v1_new # opposite direction for star 2
return pos1_new, pos2_new, v1_new, v2_new
# For plotting
pos1_list = []
pos2_list = []
for i in range(100000):
pos1, pos2, v1, v2, a1, a2 = verlet_integration(pos1, pos2, a1, a2, v1, v2)
if i % 100 == 0:
pos1, pos2, v1, v2 = peters_mathews(pos1, pos2, v1, v2)
pos1_list.append(pos1.copy())
pos2_list.append(pos2.copy())
pos1_arr = np.array(pos1_list)
pos2_arr = np.array(pos2_list)
plt.figure(figsize=(8, 8))
plt.plot(pos1_arr[:, 0], pos1_arr[:, 1], label="Star 1")
plt.plot(pos2_arr[:, 0], pos2_arr[:, 1], label="Star 2")
plt.scatter([0], [0], color='black', marker='x', label="Center of Mass")
plt.xlabel("x position (m)")
plt.ylabel("y position (m)")
plt.title("Binary system orbit simulation with Peters-Mathews decay")
plt.legend()
plt.axis('equal')
plt.grid()
plt.show()
< /code>
Hier ist mein Code, den ich gemacht habe, um zwei Neutronensterne zu simulieren, was zu einer Kilonova führt. Ich möchte, dass diese beiden Neutronensterne aufgrund der Funktion peters_mathews, die für den Orbitalverfall verantwortlich ist, näher zueinander heranwachsen. Zusätzlich berechnete die Geschwindigkeitsverwandte I ihre Orbitalposition. Anstatt dass die Umlaufbahn verfälscht und sie näher kommen, wachsen sie weiter auseinander. Wie soll ich das beheben? Auch die Beschleunigung und Geschwindigkeit wurden modifiziert, um sie näher zusammenzubringen. Ich habe auch künstlichen Zerfall implementiert, was auch nicht zu funktionieren scheint. Ich möchte, dass sie näher kommen, nicht weiter voneinander entfernt. Auch versucht, das Problem zu debuggen, kann aber nicht identifizieren.>
Was bringt diese beiden Neutronensterne nicht näher zueinander? (Kilonova -Simulation mit Python) ⇐ Python
-
- Similar Topics
- Replies
- Views
- Last post