millimeters-of-aluminum/scripts/singletons/orbital_mechanics.gd

# OrbitalMechanics.gd
extends Node

# This singleton holds all universal physics constants and functions.

# The scaled gravitational constant for the entire simulation.
const G = 10.0 # Adjust this to control the "speed" of your simulation

func calculate_n_body_gravity_forces(body_to_affect: Node2D) -> Vector2:
	var total_force = Vector2.ZERO
	if not is_instance_valid(body_to_affect):
		return total_force

	# Get the list of all major gravitational bodies from the GameManager.
	var system_data = GameManager.get_system_data()
	if not system_data:
		return total_force

	# We only consider planets and the star as major attractors for performance.
	var attractors = system_data.all_bodies()

	for attractor in attractors:
		if is_instance_valid(attractor) and attractor != body_to_affect:
			total_force += calculate_gravitational_force(body_to_affect, attractor)
	
	return total_force

# Calculates the force of gravity exerted by a primary on an orbiter.
func calculate_gravitational_force(orbiter: Node2D, primary: RigidBody2D) -> Vector2:
	if not is_instance_valid(orbiter) or not is_instance_valid(primary):
		return Vector2.ZERO
		
	var direction = orbiter.global_position.direction_to(primary.global_position)
	var distance_sq = orbiter.global_position.distance_squared_to(primary.global_position)
	
	if distance_sq < 1.0: # Avoid division by zero
		return Vector2.ZERO
		
	var force_magnitude = (G * primary.mass * orbiter.mass) / distance_sq
	return direction * force_magnitude

# Simplified n-body gravitational pull on orbiters. 
# Station orbiting a moon will for instance calculate forces to the moon, the planet, and the star.
func simple_n_body_grav(primary: RigidBody2D, orbiter : RigidBody2D) -> Vector2:
	var pull = calculate_gravitational_force(primary, orbiter)
	var inner_primary : CelestialBody = primary

	while (inner_primary.primary is CelestialBody):
		pull = pull + calculate_gravitational_force(inner_primary.primary, orbiter)
		inner_primary = inner_primary.primary
	
	return pull
		
# Calculates the perfect initial velocity for a stable circular orbit.
func calculate_circular_orbit_velocity(orbiter: Node2D, primary: RigidBody2D) -> Vector2:
	if not is_instance_valid(primary):
		return Vector2.ZERO

	var distance = orbiter.global_position.distance_to(primary.global_position)
	if distance == 0:
		return Vector2.ZERO

	# v = sqrt(G * M / r)
	var speed_magnitude = sqrt(G * primary.mass / distance)
	
	var direction_to_orbiter = primary.global_position.direction_to(orbiter.global_position)
	var perpendicular_direction = Vector2(direction_to_orbiter.y, -direction_to_orbiter.x)
	
	# The final velocity must include the primary's own velocity.
	return (perpendicular_direction * speed_magnitude) + primary.linear_velocity
	
func _calculate_n_body_orbital_path(body_to_trace: OrbitalBody2D) -> PackedVector2Array:
	var num_steps = 10
	var time_step = 60
	
	var ghost_position = body_to_trace.global_position
	var ghost_velocity = body_to_trace.linear_velocity
	
	var path_points = PackedVector2Array()

	for i in range(num_steps):
		# Create a temporary "ghost" body to calculate forces on.
		var ghost_body = OrbitalBody2D.new()
		ghost_body.global_position = ghost_position
		ghost_body.mass = body_to_trace.mass
		
		# Use our library to get the total gravitational force at the ghost's position.
		var total_force = calculate_n_body_gravity_forces(ghost_body)
		var acceleration = total_force / ghost_body.mass
		
		ghost_velocity += acceleration * time_step
		ghost_position += ghost_velocity * time_step
		path_points.append(ghost_position)
		
		ghost_body.free() # Clean up the temporary node
		
	return path_points

# Calculates an array of points for the orbit RELATIVE to the primary body.
func _calculate_relative_orbital_path(body_to_trace: CelestialBody) -> PackedVector2Array:
	if not is_instance_valid(body_to_trace) or not is_instance_valid(body_to_trace.primary):
		return PackedVector2Array()
	
	# --- Initial State ---
	var primary = body_to_trace.primary
	var primary_mass = primary.mass
	var body_mass = body_to_trace.mass
	
	# The position of the body relative to its parent.
	var ghost_relative_pos = body_to_trace.global_position - primary.global_position
	# The velocity of the body relative to its parent's velocity.
	var ghost_relative_vel = body_to_trace.linear_velocity - primary.linear_velocity
	
	# --- NEW: Dynamically Calculate Simulation Time from Orbital Period ---
	var r_magnitude = ghost_relative_pos.length()
	if r_magnitude == 0:
		return PackedVector2Array()

	var v_sq = ghost_relative_vel.length_squared()
	var mu = OrbitalMechanics.G * primary_mass # Standard Gravitational Parameter

	# 1. Calculate the specific orbital energy. Negative energy means it's a stable orbit.
	var specific_energy = v_sq / 2.0 - mu / r_magnitude

	# --- Simulation Parameters ---
	var num_steps = 200  # The desired number of points for the orbit line's smoothness.
	var time_step: float # The duration of each step, which we will now calculate.

	if specific_energy >= 0:
		# Escape trajectory (parabolic or hyperbolic). The period is infinite.
		# We'll just draw a segment of its path using a fixed time_step.
		time_step = 0.1
	else:
		# Stable elliptical orbit.
		# 2. Calculate the semi-major axis from the energy.
		var semi_major_axis = -mu / (2.0 * specific_energy)
		
		# 3. Calculate the orbital period using Kepler's Third Law: T = 2π * sqrt(a³/μ)
		var orbital_period = 2.0 * PI * sqrt(pow(semi_major_axis, 3) / mu)
		
		# 4. Calculate the time_step needed to complete one period in num_steps.
		time_step = orbital_period / float(num_steps)

	var path_points = PackedVector2Array()

	for i in range(num_steps):
		# --- Physics Calculation (Primary is now at the origin (0,0)) ---
		var distance_sq = ghost_relative_pos.length_squared()
		
		if distance_sq < 1.0:
			break
			
		# Direction is simply towards the origin.
		var direction = -ghost_relative_pos.normalized()
		var force_magnitude = (G * primary_mass * body_mass) / distance_sq
		var force_vector = direction * force_magnitude
		
		var acceleration = force_vector / body_mass
		
		# --- Integration (Update ghost's relative state) ---
		ghost_relative_vel += acceleration * time_step
		ghost_relative_pos += ghost_relative_vel * time_step
		
		path_points.append(ghost_relative_pos)
		
	return path_points