"""
Test 2: Velocity addition under the Euclidean unit-4-speed constraint.

Setup:
  Observer A at rest in the ambient frame.
  Observer B moves at velocity v (in x) relative to A.
  Particle P moves at velocity w (in x) relative to B.
  Question: what velocity does A measure for P?

In SR: v_AP = (v + w) / (1 + v*w/c^2)
In Galileo: v_AP = v + w

In our Euclidean model, both B and P are unit-4-speed worldlines.
B's 4-velocity in A's ambient frame: (v, u_B) with u_B = sqrt(c^2 - v^2).

The question is: how do we represent "P moves at velocity w relative to B"?

Option 1: Euclidean rotation in the (x, w4) plane.
  If B's 4-velocity makes angle theta_B = arctan(v / u_B) with the w4 axis,
  then "w relative to B" means P's 4-velocity makes angle theta_B + theta_w
  with the w4 axis, where tan(theta_w) = w / sqrt(c^2 - w^2).

  This is the direct Euclidean reading: rotations compose additively in angle.

Option 2: Some other composition law.
  If the geometry is really Euclidean, Option 1 is the natural reading, and
  any deviation from it requires additional structure.

Let's compute Option 1 and compare to SR.
"""

import numpy as np

c = 1.0

def angle_from_v(v):
    """Angle from w4 axis for velocity v (v < c)."""
    u = np.sqrt(c**2 - v**2)
    return np.arctan2(v, u)

def v_from_angle(theta):
    """Velocity corresponding to angle theta from w4 axis."""
    # v_3d / c = sin(theta), u/c = cos(theta)
    return c * np.sin(theta)

def sr_addition(v, w):
    return (v + w) / (1 + v*w/c**2)

def galilean_addition(v, w):
    return v + w

def euclidean_rotation_addition(v, w):
    """If rotations in (x, w4) plane compose additively in angle."""
    theta_v = angle_from_v(v)
    theta_w = angle_from_v(w)
    return v_from_angle(theta_v + theta_w)

print("=" * 80)
print("TEST 2: Velocity addition under Euclidean rotation composition")
print("=" * 80)
print()
print(f"{'v/c':>6} {'w/c':>6} {'Galilean':>12} {'Euclid. rot.':>14} {'SR':>12} {'Eucl.=SR?':>12}")
print(f"{'-'*6:>6} {'-'*6:>6} {'-'*12:>12} {'-'*14:>14} {'-'*12:>12} {'-'*12:>12}")

test_cases = [
    (0.1, 0.1),
    (0.3, 0.3),
    (0.5, 0.5),
    (0.7, 0.3),
    (0.8, 0.8),
    (0.9, 0.9),
    (0.99, 0.5),
    (0.5, 0.5),
]

for v, w in test_cases:
    gal = galilean_addition(v, w)
    eucl = euclidean_rotation_addition(v, w)
    sr = sr_addition(v, w)
    match = "YES" if abs(eucl - sr) < 1e-10 else f"NO (diff={eucl-sr:+.4f})"
    # If Euclidean gives v > c, flag it
    if abs(eucl) >= 1.0:
        eucl_str = f"{eucl:.4f} [>=c!]"
    else:
        eucl_str = f"{eucl:.4f}"
    print(f"{v:>6.2f} {w:>6.2f} {gal:>12.4f} {eucl_str:>14} {sr:>12.4f} {match:>12}")

print()
print("Expected behaviour:")
print("  If the Euclidean rotation rule is equivalent to SR:")
print("    tan(theta_v + theta_w) should relate to SR velocity addition.")
print("  Actually sin(theta_v + theta_w) = sin(theta_v)cos(theta_w) + cos(theta_v)sin(theta_w)")
print("         = (v/c)*sqrt(1-w^2/c^2) + sqrt(1-v^2/c^2)*(w/c)")
print("  This is NOT the SR addition formula.")
print()

# Let's check what happens for v = w = 0.8 specifically
v, w = 0.8, 0.8
theta_v = angle_from_v(v)
theta_w = angle_from_v(w)
print(f"Worked example: v = w = {v}")
print(f"  theta_v = {np.degrees(theta_v):.2f} deg")
print(f"  theta_w = {np.degrees(theta_w):.2f} deg")
print(f"  theta_v + theta_w = {np.degrees(theta_v + theta_w):.2f} deg")
print(f"  Euclidean result: v_3d = c*sin(sum) = {c * np.sin(theta_v + theta_w):.4f}")
print(f"  SR result:                         = {sr_addition(v, w):.4f}")
print(f"  Galilean:                          = {v + w:.4f}")
print()
print("Note: Euclidean rotation gives sin(theta_v + theta_w) which can exceed 1,")
print("      meaning the 'composed' 3D velocity would exceed c.")
print("      When theta_v + theta_w > 90 deg, we'd actually need to reflect")
print("      through the w4=0 plane, which corresponds to the forbidden region.")
