Hi am trying to make a space game with realistic orbital mechanics similar to that of KSP. I have a simple gravity system I am replicating it using line forces which works pretty well. What Im trying to do now is orbital prediction using Keplers laws of plantery motion and some physics we can calculate orbital parameters using just the “Satellites” initial Position and velocity, Semimajor Axis, Eccentricity give us the shape of the orbit. Then we need three other parameters Longitude of Ascending node, Inclination and Argument Of Periapsis these are used in the render to apply a rotation matrix. The problem is one of the parameters Argument of Periapsis isn’t calculated correctly and I have no Idea why. For more background information we are simulating this bodies motion
Its a tiny cube meant to represent our “Satellite” I have assigned it a random linear velocity and from there we try predict its movement around the glowly star object you see in the background. We ignore the force the Satellite exerts on the star. From there the orbits an ellipse and we can calculate the two parameters that tell us about the shape of the orbit Semimajor axis refers to half the ellipse along its longer side. Eccentricity really is just the relationship between the Short side and the long side of the Ellipse. After that we try to Render the orbit in 3D using the 3 other parameters I mentioned we define a Rotation matrix. Two of the three rotations work as you can see here Teal is the predicted orbit and green is the actual orbit
it lines up perfectly along 2 dimensions. Unfortunutely the last rotation using the parameter “Argument of periapsis” doesn’t really work
and I have the belief that its calculated wrong some how because I can guess the value from 0-360 untill its correct as you can see here
This here is the Code that I use to calculate it
local function CalculateNodeDot(n, e)
return n:Dot(e)
end
local nMag = n.magnitude
local eMag = e.magnitude
local DotEN = CalculateNodeDot(n, e)
local smallOmega = math.acos(DotEN / (eMag * nMag))
local omega = smallOmega
-- hours wasted counter on this = 4
This here is the formula so not sure what Im doing wrong
and yes I am 99% sure my inputs are correct.
Im thinking it could be the rotation matrix here is the full code. It is very messy and not optimized right now.
local G = 6.67430e-11 -- Gravitational constant in m^3 kg^-1 s^-2
local M = 4.225696418e15 -- Mass of the star in kg
local STUDS_TO_METERS = 3.7 -- Conversion factor from studs to meters
local part = game.Workspace.moon
local star = game.Workspace.star
-- Function to convert Roblox position to meters
local function studsToMeters(positionStuds)
return positionStuds / STUDS_TO_METERS
end
-- Function to calculate specific angular momentum vector
local function calculateSpecificAngularMomentum(r, v)
return r:Cross(v)
end
-- Function to calculate eccentricity vector
local function calculateEccentricityVector(r, v, h, mu)
return (v:Cross(h) / mu) - r.unit
end
-- Function to calculate eccentricity
local function calculateEccentricity(e_vector)
return e_vector.magnitude
end
-- Function to calculate semi-major axis using vis-viva equation
local function calculateSemiMajorAxis(r_mag, v_mag, mu)
return 1 / ((2 / r_mag) - (v_mag^2 / mu))
end
-- Function to calculate periapsis distance
local function calculatePeriapsisDistance(a, e)
return a * (1 - e)
end
local partPosition = part.Position
local partVelocity = part.Velocity
local starPosition = star.Position
-- Gravitational parameter
local mu = G * M -- in m^3/s^2
-- Calculate r and v in meters
local r = studsToMeters(partPosition - starPosition)
local v = studsToMeters(partVelocity)
-- Calculate specific angular momentum vector
local h = calculateSpecificAngularMomentum(r, v)
-- Calculate eccentricity vector
local e_vector = calculateEccentricityVector(r, v, h, mu)
local h_vector = calculateSpecificAngularMomentum
-- Calculate eccentricity
local eccentricity = calculateEccentricity(e_vector)
-- Calculate the magnitudes of r and v
local r_mag = r.magnitude
local v_mag = v.magnitude
-- Calculate semi-major axis
local semiMajorAxis = calculateSemiMajorAxis(r_mag, v_mag, mu)
-- Calculate periapsis distance
local periapsisDistance = calculatePeriapsisDistance(semiMajorAxis, eccentricity)
-- Calculate perihelion position
local perihelionPosition = starPosition + studsToMeters(periapsisDistance * e_vector.unit * STUDS_TO_METERS)
-- Print results
print("Specific Angular Momentum Vector (h):", h)
print("Eccentricity Vector (e):", e_vector)
print("Eccentricity:", eccentricity)
print("Semi-Major Axis :", semiMajorAxis*STUDS_TO_METERS)
print("Periapsis Distance (meters):", periapsisDistance*STUDS_TO_METERS)
print("Perihelion Position:", perihelionPosition*STUDS_TO_METERS)
-- Create a big red part at the calculated perihelion position
local function createPartAtPerihelion()
local bigRedPart = Instance.new("Part")
bigRedPart.Size = Vector3.new(100, 100, 100) -- Example size
bigRedPart.Position = perihelionPosition -- Set position to perihelion position
bigRedPart.Color = Color3.fromRGB(255, 0, 0) -- Red color
bigRedPart.Parent = game.Workspace -- Parent to the workspace or appropriate parent
bigRedPart.CanCollide = false
end
-- Call the function to create the part at the calculated perihelion position
createPartAtPerihelion()
-- Function to calculate semi-major axis in meters
local function calculateSemiMajorAxis(partPosition, partVelocity, starPosition, mu)
-- Calculate the relative position vector in meters
local r = studsToMeters(partPosition - starPosition)
-- Calculate the relative velocity vector in meters per second
local v = studsToMeters(partVelocity)
-- Calculate specific orbital angular momaentum vector
local h = r:Cross(v)
-- Calculate the magnitude of r and v
local r_mag = r.magnitude
local v_mag = v.magnitude
-- Calculate the specific orbital energy
local E = (v_mag^2)/2 - mu/r_mag
-- Calculate the semi-major axis in meters
local a = -mu/(2*E)
-- Convert semi-major axis back to Roblox studs and return
return a * STUDS_TO_METERS
end
local mu = G * M -- in m^3/s^2
-- Calculate semi-major axis in Roblox studs
local semiMajorAxisStuds = calculateSemiMajorAxis(partPosition, partVelocity, starPosition, mu)
local a = semiMajorAxisStuds
function calculateRa(a, perihelionPosition)
local Ra = (2 * a)-perihelionPosition
return Ra
end
local e = eccentricity -- Replace with your value for eccentricity
-- Calculate parameter p
local p = a * (1 - e^2)
-- Calculate aphelion distance Ra
local Ra = p / (1 - e)
local aphelionDistance = Ra
-- Print the results
local a = semiMajorAxisStuds
-- Function to calculate the aphelion position
local aphelionDistance = p/(1-e)
print("Parameter (p):", p)
print("Aphelion distance", aphelionDistance)
-- Define a function to calculate orbital elements
function calculateOrbitalElements(position, velocity, mu)
-- Calculate specific angular momentum vector
local h = position:Cross(velocity)
-- Calculate node line vector (n)
local k = Vector3.new(0, 1, 0) -- Reference direction (y-axis)
local n = k:Cross(h)
-- Calculate eccentricity vector (e)
local r = position.Magnitude
local vCrossH = velocity:Cross(h)
local e = (vCrossH / mu) - (position / r)
-- Calculate inclination (i)
local hMag = h.Magnitude
local i = math.acos(h.Y / hMag)
-- Calculate longitude of ascending node (Omega)
local nMag = n.Magnitude
local Omega = 0
if nMag ~= 0 then
Omega = math.acos(n.X / nMag)
if n.Z < 0 then
Omega = 2 * math.pi - Omega
end
end
local function CalculateNodeDot(n, e)
return n:Dot(e)
end
local nMag = n.magnitude
local eMag = e.magnitude
local DotEN = CalculateNodeDot(n, e)
local smallOmega = math.acos(DotEN / (eMag * nMag))
local omega = smallOmega
-- hours wasted counter on this = 4
return math.deg(i), math.deg(Omega), math.deg(omega) -- Return angles in degrees
end
-- Example values
local position = partPosition
local velocity = partVelocity
-- Calculate orbital elements
local inclination, longitudeOfAscendingNode, argumentOfPerihelion = calculateOrbitalElements(position, velocity, mu)
local argumentOfPerihelion = 125
-- Print the results
print("Inclination (i):", inclination)
print("Longitude of Ascending Node (Ω):", longitudeOfAscendingNode)
print("Argument of Perihelion (ω):", argumentOfPerihelion)
local function RobloxToPeriFocal(x, y, z)
return x, -z, y
end
local function PeriFocalToRoblox(i, j, k)
return i, k, -j
end
-- Calculate the transformation matrix from perifocal to ECI coordinates
local function CalculatePerifocalToECIMatrix(inclination, longitudeOfAscendingNode, argumentOfPeriapsis)
local i = math.rad(inclination)
local Omega = math.rad(longitudeOfAscendingNode)
local omega = math.rad(argumentOfPeriapsis)
local cosI = math.cos(i)
local sinI = math.sin(i)
local cosOmega = math.cos(Omega)
local sinOmega = math.sin(Omega)
local cosOmegaSmall = math.cos(omega)
local sinOmegaSmall = math.sin(omega)
-- Calculate the rotation matrix
local m11 = cosOmega * cosOmegaSmall - sinOmega * cosI * sinOmegaSmall
local m12 = -cosOmega * sinOmegaSmall - sinOmega * cosI * cosOmegaSmall
local m13 = sinOmega * sinI
local m21 = sinOmega * cosOmegaSmall + cosOmega * cosI * sinOmegaSmall
local m22 = -sinOmega * sinOmegaSmall + cosOmega * cosI * cosOmegaSmall
local m23 = -cosOmega * sinI
local m31 = sinI * sinOmegaSmall
local m32 = sinI * cosOmegaSmall
local m33 = cosI
return m11, m12, m13, m21, m22, m23, m31, m32, m33
end
local function drawEllipticalOrbit(semiMajorAxis, eccentricity, inclination, longitudeOfAscendingNode, argumentOfPerihelion)
local numPoints = 100
-- Calculate the transformation matrix
local m11, m12, m13, m21, m22, m23, m31, m32, m33 =
CalculatePerifocalToECIMatrix(inclination, longitudeOfAscendingNode, argumentOfPerihelion)
for j = 1, numPoints do
local angle = 2 * math.pi * j / numPoints -- Calculate angle around the orbit
-- Calculate distance from the center to the point on the orbit
local distance = semiMajorAxis * (1 - eccentricity^2) / (1 + eccentricity * math.cos(angle))
-- Calculate position in the orbital plane (x-y plane)
local x = distance * math.cos(angle)
local y = distance * math.sin(angle)
local z = 0
-- Apply the rotation matrix to the point (x, y, z) in perifocal coordinates
local x3 = m11 * x + m12 * y + m13 * z
local y3 = m21 * x + m22 * y + m23 * z
local z3 = m31 * x + m32 * y + m33 * z
local finalX, finalY, finalZ = PeriFocalToRoblox(-x3, y3, z3)
local finalPosition = Vector3.new(finalX, finalY, finalZ)
local orbitPart = Instance.new("Part")
orbitPart.Size = Vector3.new(30, 30, 30)
orbitPart.Position = finalPosition
orbitPart.Color = Color3.fromRGB(0, 255, 255)
orbitPart.Anchored = true
orbitPart.Parent = game.Workspace
orbitPart.CanCollide = false
end
end
-- Draw the orbit
drawEllipticalOrbit(semiMajorAxisStuds, eccentricity, inclination, longitudeOfAscendingNode, argumentOfPerihelion)
If you have any ideas on how to approximate the Argument of Periapsis I would like to hear them too!!! If it helps the 3 parameters I mentioned have a 1:1 correspondence with Euler’s angles
