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 = 1.0 # Adjust this to control the "speed" of your simulation
const MIN_INFLUENCE_THRESHOLD = 0.00001
const ROCHE_LIMIT_MASS_MULTIPLIER = 0.5


func _physics_process(_delta: float) -> void:
	var star_system: StarSystem = GameManager.current_star_system
	if not star_system:
		return

	var star = star_system.get_star()
	var planetary_systems = star_system.get_planetary_systems()
	
	if not star:
		return
	
	var top_level_bodies: Array[OrbitalBody3D] = [star]
	top_level_bodies.append_array(planetary_systems)
	apply_n_body_forces(top_level_bodies)
	
	for system in planetary_systems:
		var system_attractors = system.get_internal_attractors()
		apply_n_body_forces(system_attractors)
		
	for star_orbiter in star_system.get_orbital_bodies():
		star_orbiter.apply_force_recursive(calculate_n_body_force(star_orbiter, top_level_bodies))

func calculate_gravitational_force(orbiter: OrbitalBody3D, primary: OrbitalBody3D) -> Vector3:
	if not is_instance_valid(orbiter) or not is_instance_valid(primary):
		# REFACTOR: Return Vector3.ZERO
		return Vector3.ZERO
		
	var distance_sq = orbiter.global_position.distance_squared_to(primary.global_position)
	
	if distance_sq < 1.0:
		return Vector3.ZERO
		
	var influence_a = primary.mass / distance_sq
	var influence_b = orbiter.mass / distance_sq
	
	if influence_a < MIN_INFLUENCE_THRESHOLD and influence_b < MIN_INFLUENCE_THRESHOLD:
		return Vector3.ZERO
	
	var force_magnitude = (G * primary.mass * orbiter.mass) / distance_sq
	# REFACTOR: direction_to is now 3D, this logic is fine
	var direction = orbiter.global_position.direction_to(primary.global_position)
	return direction * force_magnitude
	
func apply_n_body_forces(attractors: Array[OrbitalBody3D]):
	for i in range(attractors.size()):
		var body_a: OrbitalBody3D = attractors[i]
		if not is_instance_valid(body_a): continue

		for j in range(i + 1, attractors.size()):
			var body_b: OrbitalBody3D = attractors[j]
			if not is_instance_valid(body_b): continue
	
			# REFACTOR: force_vector is now Vector3
			var force_vector = calculate_gravitational_force(body_a, body_b)
			
			if force_vector != Vector3.ZERO:
				body_a.apply_force_recursive(force_vector)
				body_b.apply_force_recursive(-force_vector)

func calculate_n_body_force(body: OrbitalBody3D, attractors: Array[OrbitalBody3D]) -> Vector3:
	var total_pull: Vector3 = Vector3.ZERO
	for attractor in attractors:
		total_pull += calculate_gravitational_force(body, attractor)
	return total_pull

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

	var system_data = GameManager.get_system_data()
	if not system_data:
		return total_force

	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

func calculate_circular_orbit_velocity(orbiter: OrbitalBody3D, primary: OrbitalBody3D) -> Vector3:
	if not is_instance_valid(primary):
		return Vector3.ZERO

	var distance = orbiter.global_position.distance_to(primary.global_position)
	if distance == 0:
		return Vector3.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)
	
	# REFACTOR: This is the biggest 2D -> 3D logic change.
	# We can't just get a simple perpendicular vector. We need to define
	# an orbital plane. We'll assume a "flat" system on the XZ plane,
	# so the "up" vector is Vector3.UP.
	# We find the perpendicular by crossing "up" with the direction to the orbiter.
	var perpendicular_direction = Vector3.UP.cross(direction_to_orbiter).normalized()

	return perpendicular_direction * speed_magnitude
	
# REFACTOR: Returns PackedVector3Array
func _calculate_n_body_orbital_path(body_to_trace: OrbitalBody3D) -> PackedVector3Array:
	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 = PackedVector3Array()

	# for i in range(num_steps):
	# 	var ghost_body = OrbitalBody3D.new()
	# 	ghost_body.global_position = ghost_position
	# 	ghost_body.mass = body_to_trace.mass
		
	# 	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()
		
	return path_points

# REFACTOR: Returns PackedVector3Array
func _calculate_relative_orbital_path(body_to_trace: OrbitalBody3D) -> PackedVector3Array:
	if not is_instance_valid(body_to_trace) or not body_to_trace.has_method("get_primary") or not is_instance_valid(body_to_trace.get_primary()):
		return PackedVector3Array()

	var primary = body_to_trace.get_primary()
	var primary_mass = primary.mass
	var body_mass = body_to_trace.mass

	var ghost_relative_pos = body_to_trace.global_position - primary.global_position
	var ghost_relative_vel = body_to_trace.linear_velocity - primary.linear_velocity

	var r_magnitude = ghost_relative_pos.length()
	if r_magnitude == 0:
		return PackedVector3Array()

	var v_sq = ghost_relative_vel.length_squared()
	var mu = G * primary_mass
	var specific_energy = v_sq / 2.0 - mu / r_magnitude

	var num_steps = 200
	var time_step: float
	if specific_energy >= 0:
		time_step = 0.1
	else:
		var semi_major_axis = -mu / (2.0 * specific_energy)
		var orbital_period = 2.0 * PI * sqrt(pow(semi_major_axis, 3) / mu)
		time_step = orbital_period / float(num_steps)

	var path_points = PackedVector3Array()

	for i in range(num_steps):
		var distance_sq = ghost_relative_pos.length_squared()
		if distance_sq < 1.0:
			break
		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

		ghost_relative_vel += acceleration * time_step
		ghost_relative_pos += ghost_relative_vel * time_step
		path_points.append(ghost_relative_pos)

	return path_points
	
# --- These functions are scalar and need no changes ---

func calculate_hill_sphere(orbiter: OrbitalBody3D, primary: OrbitalBody3D) -> float:
	if not is_instance_valid(orbiter) or not is_instance_valid(primary) or primary.mass <= 0:
		return 0.0
	var distance = orbiter.global_position.distance_to(primary.global_position)
	var mass_ratio = orbiter.mass / (3.0 * primary.mass)
	if mass_ratio < 0: return 0.0
	return distance * pow(mass_ratio, 1.0/3.0)
	
func calculate_simplified_roche_limit(primary: OrbitalBody3D) -> float:
	if not is_instance_valid(primary) or primary.mass <= 0:
		return 100.0
	return sqrt(primary.mass) * ROCHE_LIMIT_MASS_MULTIPLIER

func get_orbital_time_in_seconds(orbiter: OrbitalBody3D, primary: OrbitalBody3D) -> float:
		var mu = OrbitalMechanics.G * primary.mass
		var r = orbiter.global_position.distance_to(primary.global_position)
		return TAU * sqrt(pow(r, 3) / mu)