Warum umkreist mein Photon nicht den richtigen Radius?Python

Python-Programme
Guest
 Warum umkreist mein Photon nicht den richtigen Radius?

Post by Guest »

In Blender verfolge ich den Pfad einiger 2D-Photonen, die sich aufgrund eines nicht rotierenden Standard-Schwarzen Lochs, das durch die Schwarzschild-Metrik definiert ist, verbiegen. Ich stelle meine Anfangsgeschwindigkeiten in Form von r (Radius), Phi (Winkel) und t (Zeit) in Bezug auf den affinen Parameter Lambda ein und aktualisiere dann iterativ den Raum-Zeit-Vektor basierend auf den Christoffel-Symbolen am jeweiligen Punkt. Fügen Sie einfach den Code in ein Blender-Python-Skript ein und führen Sie es aus.
Bei der Projektion aus der Ferne kreist mein Photon nicht mit 3GM/(c^2). Als Referenz: Die innere Kugel ist das Schwarze Loch mit dem Radius 2GM/(c^2). Auf der rechten Seite der Schleife tritt das Photon innerhalb des Umlaufbahnradius ein und bewegt sich nicht spiralförmig zum Ereignishorizont!
Image

Bei der anfänglichen Projektion von x = 3GM/(c^2) funktioniert es jedoch tatsächlich in die Umlaufbahn Radius! (Irgendwann entkommt es, da es anfällig für leichte Ungenauigkeiten ist)
Image

Alles um 1e6 verkleinert, um es in Blender darzustellen

Code: Select all

import bpy
from math import sin, cos, pi, sqrt, atan2
from mathutils import Vector

# Constants
M = 9e32  # Mass of the black hole in kg
c = 299792458  # Speed of light in m/s
G = 6.6743e-11  # Gravitational constant in N m^2/kg^2

# Simulation parameters
curve_num = 4000  # Resolution of the curve
#test_point = Vector((0, 4, 0))  # Initial position
test_point = Vector((3*G*M/(1e6 * c**2), 0, 0))
alpha = 4 * pi/8 # initial projection angle (at x = 0, y = 4 try 2.55 * pi/8)
d_lambda = 1e-4  # Step size for integration

# Initialize variables
r = test_point.length * 1e6  # Convert to meters
p = atan2(test_point.y, test_point.x)  # p = phi = Initial planar angle

# initial radial and angular velocities
v_r = -c * cos(alpha)  # Radial component of velocity
v_p = c * sin(alpha) / r  # Angular component of velocity

# Schwarzschild metric components (2D)
g_tt = -(1 - 2 * G * M / (r * c**2)) # time
g_rr = 1 / (1 - 2 * G * M / (r * c**2)) # radial
g_pp = r**2 # phi

# Initialize velocity in the time direction
v_t = sqrt(-(g_rr * v_r**2 + g_pp * v_p**2) / g_tt)

# Debug: Initial null geodesic condition
print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2)

# Plot black hole sphere mesh
horizon = 2 * G * M / c**2  # Event horizon radius
bpy.ops.mesh.primitive_uv_sphere_add()
black_hole = bpy.context.object
black_hole.scale = (horizon / 1e6, horizon / 1e6, horizon / 1e6)  # Scale for Blender units

# Create a 2D trajectory curve
bpy.ops.mesh.primitive_circle_add(radius=1, vertices=curve_num + 1)
curve = bpy.context.object
verts = curve.data.vertices

# Delete unnecessary vertex to make the trajectory a 2D curve
bpy.ops.object.mode_set(mode='EDIT')
bpy.ops.mesh.select_all(action='DESELECT')
bpy.ops.object.mode_set(mode='OBJECT')
verts[-1].select = True
bpy.ops.object.mode_set(mode='EDIT')
bpy.ops.mesh.delete(type='VERT')
bpy.ops.object.mode_set(mode='OBJECT')

# Main simulation loop
for i, v in enumerate(verts):
if i == 0:
v.co = test_point
else:
if r > horizon + 1000:
# Relevant Christoffel symbols
Γ_r_tt = (G * M / (r**2 * c**2)) * (1 - 2 * G * M / (r * c**2))
Γ_r_rr = -(G * M / (r**2 * c**2)) / (1 - 2 * G * M / (r * c**2))
Γ_r_pp = -r * (1 - 2 * G * M / (r * c**2))
Γ_p_rp = 1 / r

# Get accelerations
dv_r = -Γ_r_tt * v_t**2 - Γ_r_rr * v_r**2 - Γ_r_pp * v_p**2
dv_p = -Γ_p_rp * v_r * v_p

# Update velocities
v_r += dv_r * d_lambda
v_p += dv_p * d_lambda

# Update positions
r += v_r * d_lambda
p += v_p * d_lambda

# Recalculate metric coefficients
g_tt = -(1 - 2 * G * M / (r * c**2))
g_rr = 1 / (1 - 2 * G * M / (r * c**2))
g_pp = r**2

# Debug: Print null geodesic condition
print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2)

# Convert spherical to Cartesian coordinates (2D in equatorial plane)
x = (r / 1e6) * cos(p)
y = (r / 1e6) * sin(p)

# Update vertex position
v.co = Vector((x, y, 0))

# Finalize the curve
bpy.context.view_layer.objects.active = curve
bpy.ops.object.convert(target='CURVE')
curve = bpy.context.object
curve.data.bevel_depth = 0.01
curve.data.fill_mode = 'FULL'
curve.data.bevel_resolution = 12
Ich erwarte, dass meine geodätische Nulldruckanweisung innerhalb der for-Schleife entweder 0 oder eine relativ kleine Ganzzahl ist. Zum Beispiel kann die erste geodätische Null-Druckanweisung 0, 16, -8 usw. sein, aufgrund einer kleinen Ungenauigkeit bei der großen Float-Addition (e17-Größen).
Das habe ich Ich habe versucht zu debuggen, indem ich meine „else“-Anweisung durch „elif j == 1“ ersetzt habe, und als ich mir die nächste Iteration ansah, konnte ich sehen, dass die Null-Geodäte einen viel größeren Float druckt.
Ich glaube, das herauszufinden Ein geodätischer Nullfehler zeigt, warum meine Flugbahnen falsch sind. Ich habe meine Christoffel-Symbole und die Implementierung mehrmals überprüft und kann nichts Auffälliges erkennen.

Quick Reply

Change Text Case: 
   
  • Similar Topics
    Replies
    Views
    Last post