import math
pi = math.pi

print("=== Unit mapping from article ===")
print()
print("Article states: a = r_a * i * 2 * l_P")
print("where i = number of Planck mass points in center mass")
print("and r_a is the dimensionless simulation radius")
print()
print("So: a_physical = r_sim * M * 2 * l_P  (using M for center mass)")
print("But earlier we said: R = r_sim * M * l_P")
print("The factor of 2 comes from the Schwarzschild radius = 2*G*M/c^2 = 2*M*l_P")
print()
print("Let's use the article's mapping: a = r_sim * i * 2 * l_P")
print()

# GR formula: theta = 6*pi*G*M / (a*c^2*(1-e^2))
# With a = r_sim * M * 2 * l_P and G*M/c^2 = M*l_P (Schwarzschild radius / 2):
# theta = 6*pi * M*l_P / (r_sim * M * 2 * l_P * (1-e^2))
# theta = 6*pi / (2 * r_sim * (1-e^2))
# theta = 3*pi / (r_sim * (1-e^2))
print("GR theta = 6*pi*G*M / (a*c^2*(1-e^2))")
print("         = 6*pi * (M*l_P) / (r_sim*M*2*l_P*(1-e^2))")
print("         = 6*pi / (2*r_sim*(1-e^2))")
print("         = 3*pi / (r_sim*(1-e^2))")
print()
print("THIS IS EXACTLY WHAT THE CODE USES!")
print("The (M+1) factor cancels when using the article's mapping a = r_sim*M*2*l_P")
print()
print("But wait - the article uses G*M, not G*(M+m).")
print("The standard GR formula uses M (just the central mass) when M >> m.")
print("When M is not >> m, the correct formula uses (M+m).")
print()

# Let's check with (M+m):
# theta = 6*pi*G*(M+m) / (a*c^2*(1-e^2))
# With a = r_sim * M * 2 * l_P:
# theta = 6*pi*(M+m)*l_P / (r_sim*M*2*l_P*(1-e^2))
# theta = 3*pi*(M+m) / (r_sim*M*(1-e^2))
# theta = 3*pi*(1 + m/M) / (r_sim*(1-e^2))
print("With (M+m):")
print("  theta = 3*pi*(M+m) / (r_sim*M*(1-e^2))")
print("        = 3*pi*(1 + 1/M) / (r_sim*(1-e^2))")
print()

M = 24
e = 0.931230421
r = 79590.8651832
dt = 0.012767161

gr_code = 3*pi / (r*(1-e*e))
gr_with_mass = 3*pi*(1 + 1/M) / (r*(1-e*e))
print(f"For M={M}:")
print(f"  Code GR (per half orbit) = {gr_code:.9f}")
print(f"  With (1+1/M) correction  = {gr_with_mass:.9f}")
print(f"  Simulated precession     = {dt}")
print(f"  Amplification (code)     = {dt/gr_code:.5f}")
print(f"  Amplification (with 1+1/M) = {dt/gr_with_mass:.5f}")
print()

# The (1+1/M) correction is small: for M=24, it's 1.042
# For M=56, it's 1.018
# This doesn't change the amplification significantly
print("The (1+1/M) correction is small:")
for M_val in [24, 32, 40, 48, 56]:
    print(f"  M={M_val}: factor = {1+1/M_val:.4f}")
print()
print("So the code's GR formula is essentially correct")
print("(neglecting the small m/M correction).")
print()
print("The amplification values in the table ARE correct.")
print()

# Now let's also check the article's equation 4 vs equation in the new text
print("=" * 60)
print("CHECK: Article equation 4 vs new text equation")
print()
print("Article Eq 4 (line 391): theta = 6*pi*G*M / (a*(1-e^2)*c^2)")
print("  Uses just M (central mass)")
print()
print("New text equation: delta_phi = 6*pi*G_N*(M+m) / (c^2*a*(1-e^2))")
print("  Uses (M+m)")
print()
print("These are INCONSISTENT. The article's own Eq 4 uses M only,")
print("matching the code. The new text should also use M only,")
print("or explain the (M+m) version.")
print()

# Now verify the Kepler claim
print("=" * 60)
print("KEPLER PERIOD CLAIM REVIEW")
print()
print("Source code line 123:")
print("  calc_period = 16*pi*alpha^1.5 * jmax^3 / (ipoints^2.5 * sqrt(jpoints))")
print()
print("This is computed BEFORE the simulation runs (line 123, sim starts ~183)")
print("It is labeled 'Theoretical Metrics (circular orbit)' (line 135)")
print()

M = 24
alpha = 137.0359991770
kr = 12
ipoints = M
jpoints = M + 1
jmax = kr * ipoints + 1
r = 2 * alpha * 2 * jmax**2 / ipoints**2
T_code = 16 * pi * alpha**1.5 * jmax**3 / (ipoints**2.5 * jpoints**0.5)
T_kepler = math.sqrt(4 * pi**2 * r**3 * M / (M+1))

print(f"  calc_period = {T_code:.1f}")
print(f"  Kepler (dimensionless) = sqrt(4*pi^2*r^3*M/(M+1)) = {T_kepler:.1f}")
print(f"  Ratio = {T_code/T_kepler:.10f}")
print()
print("  These are ALGEBRAICALLY IDENTICAL.")
print("  calc_period IS the Kepler formula, not an independent check.")
print()

# What about the actual simulation output?
print("Actual simulation sidereal half-periods (from raw data):")
halves = [69562247, 69562238, 69562229, 69562221, 69562212, 69562202, 69562192]
avg_half = sum(halves)/len(halves)
print(f"  Average = {avg_half:.1f}")
print(f"  Kepler half = {T_code/2:.1f}")
print(f"  Error = {abs(avg_half - T_code/2)/(T_code/2)*100:.4f}%")
print()
print("  BUT: these sidereal half-periods are for the ELLIPTICAL orbit (e~0.93)")
print("  The Kepler period for an ellipse depends on semi-major axis a,")
print("  not the circular radius. For this highly eccentric orbit,")
print("  the comparison is not straightforward.")