OrbitCalculations.lua

-- Several portions of this code are ports from MuMechLib, and as such, this code is licensed under GPL.
-- Thanks to braeunig.us/space for much of the math.

--TODO: make library for calculating AN/DN, anomalies, times, etc.:
-- function DOSTUFF(subjectOrbitParameters[, other arguments])
-- OrbitParameters is a table of:
-- semi-major axis (meters) - SMA
-- periapsis (radius) (meters) - PE
-- apoapsis (radius) (meters) - AP
-- period (seconds) - P
-- eccentricity - E
-- inclination (degrees) - I
-- argument of periapsis (degrees) - AoP
-- Longitude of Ascending Node (degrees) - LAN
-- time to periapsis (seconds) - PeM
-- It should be possible to compute everything from these. Anything requiring further input should be in it's own script.
-- (must be GPL for including ports of MuMechLib C# functions)

-- DEBUG
--

-- TAKES: 3d vector A, 3d vector B
-- RETURNS: cross-product (3d vector)
function VectorCross(a, b)
    local c = {0,0,0}
    c[1] = a[2]*b[3] - a[3]*b[2]
    c[2] = a[3]*b[1] - a[1]*b[3]
    c[3] = a[1]*b[2] - a[2]*b[1]
    return c
end

-- TAKES: orbital position vector (3d vector)
-- RETURNS: geocentric latitude and longitude? (2d vector) at given orbital element vector
function LatLongofVector(vinput)
    -- get the geocentric latitude and longitude of a orbital element vector
    local rads = math.pi/180
    local c1 = vinput[1]
    local c2 = vinput[2]
    local c3 = vinput[3]
    local lat = math.atan(c3 / math.sqrt(c1*c1 + c2*c2))/rads
    local long = math.atan2(c2, c1)/rads
    if (c1 < 0) then
        long = long + 90
    elseif (c1 > 0) then
        long = long + 270
    end
    local coord = {lat, long}
    return coord;
end

-- TAKES: Semi-major axis, eccentricity
-- RETURNS: Periapsis radius
function GetPe(SMA, E)
    return SMA * (1 - E)
end

-- TAKES: Semi-major axis, eccentricity
-- RETURNS: Apoapsis radius
function GetAp(SMA, E)
    return SMA * (1 + E)
end

-- TAKES: true anomaly (radians)
-- RETURNS: eccentric anomaly (radians) at given true anomaly
function GetEccentricAnomalyAtTruAnomaly(orbitParams, truAnomaly)
    local anom = math.atan2(math.sqrt(1 - orbitParams.E*orbitParams.E) * math.sin(truAnomaly), orbitParams.E + math.cos(truAnomaly))
    if anom < 0 then
        return anom + 2*math.pi
    else
        return anom
    end
end

-- TAKES: eccentric anomaly (radians)
-- RETURNS: mean anomaly (radians) at given eccentric anomaly
function GetMeanAnomalyAtEccentricAnomaly(orbitParams, eccentricAnomaly)
    return eccentricAnomaly - orbitParams.E*math.sin(eccentricAnomaly)
end

-- TAKES: mean anomaly (radians)
-- RETURNS: time (PE minus) at given anomaly
function GetPeMinusAtMeanAnomaly(orbitParams, meanAnomaly)
    local SpR = (orbitParams.P)/(2*math.pi) -- seconds per radian
    return orbitParams.P - meanAnomaly * SpR
end

-- TAKES: true anomaly (radians)
-- RETURNS: time (PE minus) at given anomaly
function GetPeMinusAtTruAnomaly(orbitParams, truAnomaly)
    return GetPeMinusAtMeanAnomaly(orbitParams, GetMeanAnomalyAtEccentricAnomaly(orbitParams, GetEccentricAnomalyAtTruAnomaly(orbitParams, truAnomaly)))
end

-- TAKES: true anomaly (radians)
-- RETURNS: time (seconds) until orbit reaches true anomaly
function GetTimeToTruAnomaly(orbitParams, truAnomaly)
    local t = GetPeMinusAtTruAnomaly(orbitParams, truAnomaly)
    local now = orbitParams.PeM
    if (now > t) then
        return now - t
    else
        return orbitParams.P - t + now
    end
end

-- TAKES: geocentric coordinates (latitude, longitude)
-- RETURNS: true anomaly (radians) when orbit is over coordinates
function GetTrueAnomalyAtGeoCoordinates(orbitParams, coords)
    local temp2 = math.deg(math.atan(math.tan(math.rad(coords[2]))/math.cos(math.rad(orbitParams.I))))
--  print("temp2: " .. temp2)
    local truAnomaly = 180 + (temp2 - orbitParams.AoP)
    if (truAnomaly > 360) then truAnomaly = truAnomaly - 360 end
    return math.rad(truAnomaly)
end

-- TAKES: true anomaly (radians)
-- RETURNS: radius (meters)
function GetRadiusAtTruAnomaly(orbitParams, truAnomaly)
    return orbitParams.SMA * (1-orbitParams.E*orbitParams.E)/(1+orbitParams.E*math.cos(truAnomaly))
end

-- TAKES: target orbital parameters
-- RETURNS: (radians) true anomaly of ascending node
-- note: does not work yet.
function FindAN(orbitParams, targetParams)
end

-- TAKES: target orbital parameters
-- RETURNS: (radians) true anomaly of ascending node
-- note: does not work yet.
--function FindAN(orbitParams, targetParams)
--  print("finding AN")
--  local vA = math.sin(math.rad(orbitParams.I))*math.cos(math.rad(orbitParams.LAN))
--  local vB = math.sin(math.rad(orbitParams.I))*math.sin(math.rad(orbitParams.LAN))
--  local vC = math.cos(math.rad(orbitParams.I))
--  local tA = math.sin(math.rad(targetParams.I))*math.cos(math.rad(targetParams.LAN))
--  local tB = math.sin(math.rad(targetParams.I))*math.sin(math.rad(targetParams.LAN))
--  local tC = math.cos(math.rad(targetParams.I))
    -- vector3D.cross(v, t) inline here
--  local c = {}
--  c[1] = vB*tC - vC*tB
--  c[2] = vC*tA - vA*tC
--  c[3] = vA*tB - vB*tA
--  print("end vector stuff C")

--  local coordsA = LatLongofVector(c)
--  coordsA[2] = coordsA[2] - orbitParams.LAN --don't ask me.
--  if (coordsA[2] < 0) then
--      coordsA[2] = coordsA[2] + 360
--  end
--  print("AN coords: ".. coordsA[1], coordsA[2])
    -- Got latitude/longitude of node A
--  local latB = coordsA[1] * (-1)
--  local longB = coordsA[2] + 180
--  if (longB > 360) then
--      longB = longB - 360
--  end
--  local coordsB = {latB, longB}
    -- Got latitude/longitude of node B
--  print("end coordinate stuff C")

--  local nodeA = GetTrueAnomalyAtGeoCoordinates(orbitParams, coordsA)
--  local nodeB = GetTrueAnomalyAtGeoCoordinates(orbitParams, coordsB)
--  print("end node stuff A")
    -- TODO: figure out a way to distinguish ascending and descending nodes.
--  if (GetRadiusAtTruAnomaly(orbitParams, nodeA) > GetRadiusAtTruAnomaly(orbitParams, nodeB)) then
--      return nodeA
--  else
--      return nodeB
--  end
--end

-- TAKES: target orbital parameters
-- RETURNS: (radians) angle change required to match planes
function AngleChangeRequired(orbitParams, targetParams)
    local vA = math.sin(math.rad(orbitParams.I))*math.cos(math.rad(orbitParams.LAN))
    local vB = math.sin(math.rad(orbitParams.I))*math.sin(math.rad(orbitParams.LAN))
    local vC = math.cos(math.rad(orbitParams.I))
    local tA = math.sin(math.rad(targetParams.I))*math.cos(math.rad(targetParams.LAN))
    local tB = math.sin(math.rad(targetParams.I))*math.sin(math.rad(targetParams.LAN))
    local tC = math.cos(math.rad(targetParams.I))
    local dAngle = math.acos(vA*tA + vB*tB + vC*tC)
    return dAngle
end