import numpy as np

# Real Apophis inbound state vector from paper around ~180k miles
r_apophis = np.array([-2.36202205e5, -1.50660503e5, -7.52367961e4])
v_apophis = np.array([4.33289446, 3.90501392, 1.69080403])

# Target distance and speed
target_dist_miles = 182000
target_speed_mph = 2660

target_dist_km = target_dist_miles * 1.609344
target_speed_km_s = target_speed_mph * (1.609344 / 3600)

# Normalize and scale position to exactly 182,000 miles
unit_r = r_apophis / np.linalg.norm(r_apophis)
new_r = unit_r * target_dist_km

# Normalize and scale velocity to exactly 2,660 mph keeping Apophis direction
unit_v = v_apophis / np.linalg.norm(v_apophis)
new_v = unit_v * target_speed_km_s

print(f"Target position km: {new_r}")
print(f"Target velocity km/s: {new_v}")
print(f"Position magnitude check miles: {np.linalg.norm(new_r)/1.609344}")
print(f"Velocity magnitude check mph: {np.linalg.norm(new_v)*3600/1.609344}")

import numpy as np
from scipy.integrate import solve_ivp

# Constants
G = 6.67430e-11 # m^3 kg^-1 s^-2
M = 5.972e24 # kg
mu = G * M / 1e9 # km^3 / s^2
R_e = 6378.1 # km

# Initial conditions from prompt logic
r_apophis_raw = np.array([-236202.205, -150660.503, -75236.7961])
v_apophis_raw = np.array([4.33289446, 3.90501392, 1.69080403])

target_dist_miles = 182000
target_speed_mph = 2660
km_per_mile = 1.609344
target_dist_km = target_dist_miles * km_per_mile
target_speed_km_s = target_speed_mph * (km_per_mile / 3600)

unit_position = r_apophis_raw / np.linalg.norm(r_apophis_raw)
unit_velocity = v_apophis_raw / np.linalg.norm(v_apophis_raw)

r0 = unit_position * target_dist_km
v0 = unit_velocity * target_speed_km_s

print("r0:", r0)
print("v0:", v0)

# Derivative function for 2-body problem
def equations(t, y):
    r = y[:3]
    v = y[3:]
    r_mag = np.linalg.norm(r)
    a = -mu * r / (r_mag**3)
    return np.concatenate([v, a])

# Event function to detect impact with Earth's surface
def impact_event(t, y):
    r = y[:3]
    return np.linalg.norm(r) - R_e

impact_event.terminal = True
impact_event.direction = -1

# Integrate
t_span = (0, 1000000) # large time span
y0 = np.concatenate([r0, v0])

sol = solve_ivp(equations, t_span, y0, events=impact_event, rtol=1e-9, atol=1e-9)

if sol.t_events[0].size > 0:
    t_impact = sol.t_events[0][0]
    y_impact = sol.y_events[0][0]
    r_impact = y_impact[:3]
    v_impact = y_impact[3:]
    print(f"Impact detected at t = {t_impact:.2f} seconds")
    print("r_impact:", r_impact)
    print("v_impact:", v_impact)
    
    # Calculate Lat, Long in Earth-Centered Inertial (ECI) framework at moment of impact
    # Note: To translate to actual geographic coordinates, we need Earth's rotation, 
    # but the prompt asks for "degrees east or west of the Trailing meridian".
    # Trailing meridian in orbit mechanics usually means the meridian 90 degrees behind the direction of motion, 
    # or let's look at the Earth's apex/trailing hemisphere.
    # The trailing meridian (or trailing longitudes) of Earth's orbit is at -90 degrees from the Apex of motion.
    # Earth moves along the +Y axis usually in sun-centered coordinates, but here it's Earth-centered.
    # Let's see how the position vector translates to lat/long.
    x, y, z = r_impact
    r_mag = np.linalg.norm(r_impact)
    lat = np.degrees(np.arcsin(z / r_mag))
    lon = np.degrees(np.arctan2(y, x))
    print(f"Inertial Latitude: {lat:.4f}")
    print(f"Inertial Longitude: {lon:.4f}")
    
    # Impact angle relative to the local horizontal surface
    # Angle between surface normal (-r_impact) and velocity vector (v_impact)
    dot_prod = np.dot(-r_impact, v_impact)
    v_mag = np.linalg.norm(v_impact)
    cos_theta = dot_prod / (r_mag * v_mag)
    impact_angle = np.degrees(np.arccos(cos_theta))
    # Local horizontal angle is 90 - angle with normal, or 90 - acccos(cos_theta) which is arcsin(cos_theta)
    impact_angle_horizontal = np.degrees(np.arcsin(dot_prod / (r_mag * v_mag)))
    print(f"Impact angle (relative to horizontal): {impact_angle_horizontal:.4f} degrees")
    
    # Let's find trailing meridian relation. In Earth-centered coordinates derived from the original vector,
    # let's determine the apex of Earth's motion if the original coordinates are in Ecliptic or ECI.
    # Apophis state vector from JPL is usually in ICRF/ECI or Heliocentric? 
    # The prompt gave: r_apophis_raw = [-236202.205, -150660.503, -75236.7961] which is roughly 292,000 km, 
    # which matches the close approach distance to Earth in Earth-centered coordinates.
    # In geocentric coordinates, the trailing meridian is the hemisphere facing away from Earth's orbital velocity vector.
    # Earth's orbital velocity points toward the Apex of Way (+Y in J2000 heliocentric roughly at equinox, but let's check).
else:
    print("No impact detected (missed Earth).")

import numpy as np
from scipy.integrate import solve_ivp

# Constants
G = 6.67430e-11   # m^3 kg^-1 s^-2
M = 5.972e24      # Earth mass (kg)
mu = G * M / 1e9  # km^3/s^2
R_e = 6378.1      # Earth Mean Radius (km)

START_DIST_MILES = 182000
START_SPEED_MPH  = 2660

dist_km = START_DIST_MILES * 1.609344
speed_km_s = START_SPEED_MPH * (1.609344 / 3600)

# True directions of Apophis 2029 approach
r_apophis_dir = np.array([-236202.205, -150660.503, -75236.7961])
v_apophis_dir = np.array([4.33289446, 3.90501392, 1.69080403])

# Initial State
r0 = (r_apophis_dir / np.linalg.norm(r_apophis_dir)) * dist_km
v0 = (v_apophis_dir / np.linalg.norm(v_apophis_dir)) * speed_km_s

def equations(t, state):
    r, v = state[:3], state[3:]
    a = -mu * r / (np.linalg.norm(r)**3)
    return np.concatenate([v, a])

def detect_impact(t, state):
    return np.linalg.norm(state[:3]) - R_e
detect_impact.terminal = True
detect_impact.direction = -1

sol = solve_ivp(equations, (0, 1e6), np.concatenate([r0, v0]), events=detect_impact, rtol=1e-9, atol=1e-9)

# Initial State Vectors
print("--- START VECTORS ---")
print(f"r0: {r0[0]:.4f}, {r0[1]:.4f}, {r0[2]:.4f}")
print(f"v0: {v0[0]:.6f}, {v0[1]:.6f}, {v0[2]:.6f}")

# Final State Vectors (at impact event)
if sol.y_events and len(sol.y_events[0]) > 0:
    impact_state = sol.y_events[0][0]
    rf = impact_state[:3]
    vf = impact_state[3:]
    print("\n--- IMPACT VECTORS ---")
    print(f"rf: {rf[0]:.4f}, {rf[1]:.4f}, {rf[2]:.4f}")
    print(f"vf: {vf[0]:.6f}, {vf[1]:.6f}, {vf[2]:.6f}")
    print(f"Final Magnitude check (Radius): {np.linalg.norm(rf):.2f} km")
    print(f"Final Magnitude check (Velocity): {np.linalg.norm(vf):.4f} km/s")
else:
    print("No impact detected in integration.")

# Let's perform a precise calculation.
# At any given UTC time, the Greenwich Meridian (0 degrees) has a local time equal to the UTC time.
# The Trailing Meridian is defined as the meridian where the local solar time is 18:00 (6:00 PM).
# Local Time = UTC_Time + (Longitude / 15)
# Therefore, 18.0 = UTC_Time_in_hours + (Longitude_of_Trailing_Meridian / 15)
# Longitude_of_Trailing_Meridian = (18.0 - UTC_Time_in_hours) * 15

# Let's check the UTC time at impact:
# 40.2 hours after 08:51 AM means:
# 08:51 AM + 40 hours 12 minutes = 01:03 AM (two days later)
utc_hour = 1.05 # 1 hour and 3 minutes

lon_trailing = (18.0 - utc_hour) * 15
# If longitude > 180, normalize to -180 to 180 or standard East/West.
if lon_trailing > 180:
    lon_trailing -= 360

# The target longitude is 31.27 degrees East of this Trailing Meridian
target_lon = lon_trailing + 31.27

print(f"Trailing Meridian Longitude: {lon_trailing:.4f}°")
print(f"Target Longitude: {target_lon:.4f}°")

# Re-calculating with the updated direction (31.27 degrees WEST of Trailing Meridian)
# Initial time: 08:51 AM
# Elapsed time: 40.2 hours (40 hours and 12 minutes)
# Final local time at Prime Meridian (0 deg longitude):
# 08:51 + 40 hours 12 minutes = 49 hours 03 minutes.
# 49 hours 03 minutes from the starting day's 00:00 means:
# Day 1: 08:51
# Day 2: 08:51 (24 hours)
# Day 2: 24:00 (39 hours 09 mins)
# Day 3: 01:03 AM (40 hours 12 mins from 08:51)

# So the Prime Meridian is at 01:03 AM (01.05 hours decimal).
# The Trailing Meridian is always at 6:00 PM (18:00 hours) local solar time.
# Time difference between 01:03 AM and 6:00 PM (previous day/behind):
# 18:00 is 7 hours and 3 minutes behind 01:03 AM (since 01:03 - 7h03m = 18:00).
# Therefore, Trailing Meridian longitude = -7.05 hours * 15 deg/hour = -105.75 degrees (105.75° W).

trailing_meridian_lon = -105.75
offset_west = 31.27

# Moving WEST means subtracting from the longitude
final_lon = trailing_meridian_lon - offset_west

# Wrap longitude to (-180, 180]
if final_lon <= -180:
    final_lon += 360

print(f"Final Longitude: {final_lon:.2f}")