SimpleOrbitalMechanics0.001

-- SimpleOrbitalMechanics is a purely calculation based lua library, with some constants hardcoded for Kerbal Space Program
-- v0.001
-- to use, put
-- require "SimpleOrbitalMechanics"
-- in your lua script, or run it with
-- dofile("SimpleOrbitalMechanics")
-- Original Work Copyright 2012, Nadrek, insofar as any elements are protectable by copyright
-- licensed under your choice of GPLv2, GPLv3, or Creative Commons Attribution-ShareAlike 3.0 Unported (CC BY-SA 3.0)  http://www.gnu.org/licenses/gpl-2.0.html  or  http://www.gnu.org/licenses/gpl.html  or  http://creativecommons.org/licenses/by-sa/3.0/
-- KSP specific functions will have KSP in the name; other functions will not, and can be used for real life calculations as well.

-- This script is primarily reference material and pure math; there is no dependency on KSP.
-- This script is useful not only for playing with Kerbal Space Program (KSP), but also for calculating
--   actual real life orbital mechanics around a variety of celestial bodies.
-- While there are both RL and KSP data here, all data is entirely hardcoded, and no one kind of data
--   will interfere with other kinds of use - the Real Life (RL) data will not interfere with playing KSP,
--   and the KSP data will not interfere with RL usage.
-- Orbital references include http://what-when-how.com/remote-sensing-from-air-and-space/a-few-standard-orbits-orbital-mechanics-interlude-remote-sensing/
--     http://www.braeunig.us/space/orbmech.htm
--     http://www.faqs.org/faqs/space/math/

-- global variables used in this script will be prepended with "globalSom" to avoid conflicts with globals from other scripts.

function globalSomOverallInit()
  -- GM is the gravitational parameter (u) (used for orbiting objects of insignificant mass... like your spaceship or satellite)
  -- KSP... refers to objects in Kerbal Space Program
  -- RL... refers to objects in real life


  -- Initialize our arrays.  Yes, a two dimensional array would be better, but this is very simple to understand and write.
  globalSomCelestialBodyName = {} -- Universe_Galaxy_SolarSystem_BodyOrbitingSun_BodyOrbitingPriorBody
  globalSomCelestialBodyGM = {} -- in m^3/s^2
  globalSomCelestialBodySiderealRotationalPeriod = {} -- in seconds
  globalSomCelestialBodyMinStableAltitude = {} -- in meters
  globalSomCelestialBodyHillSphereOfInfluence = {} -- in meters
  globalSomCelestialBodyEquatorialSurfaceGravity = {} -- in m/s^2
  globalSomCelestialBodyEquatorialRadius = {} -- in meters

  -- WARNING: DO NOT DEPEND ON THE OFFSET REMAINING THE SAME!!!  USE THE NAME FIELD INSTEAD!!!

  -- KSP_Galaxy1_Kerbol_Kerbin per http://kerbalspaceprogram.com/~kerbalsp/wiki/index.php?title=Kerbin 20120813
  --  and per http://kspwiki.nexisonline.net/wiki/Kerbin 20120813
  globalSomCelestialBodyName[1] = "KSP_Galaxy1_Kerbol_Kerbin"
  globalSomCelestialBodyGM[1] = 3530461000000 -- in m^3/s^2
  globalSomCelestialBodySiderealRotationalPeriod[1] = 21600 -- in seconds (6 hours * 60 min/hr * 60 sec/min)
  globalSomCelestialBodyMinStableAltitude[1] = 70000 -- in meters
  globalSomCelestialBodyHillSphereOfInfluence[1] = 84275000 -- in meters
  globalSomCelestialBodyEquatorialSurfaceGravity[1] = 9.8068 -- in m/s^2
  globalSomCelestialBodyEquatorialRadius[1] = 600000 -- in meters

  -- KSP_Galaxy1_Kerbol_Kerbin_Mun per http://kerbalspaceprogram.com/~kerbalsp/wiki/index.php?title=Mun 20120813
  --  and per http://kspwiki.nexisonline.net/wiki/Mun 20120813
  globalSomCelestialBodyName[2] = "KSP_Galaxy1_Kerbol_Kerbin_Mun"
  globalSomCelestialBodyGM[2] = 65136000000 -- in m^3/s^2
  globalSomCelestialBodySiderealRotationalPeriod[2] = 147600 -- in seconds (41 hours * 60 min/hr * 60 sec/min)
  globalSomCelestialBodyMinStableAltitude[2] = 4700 -- in meters
  globalSomCelestialBodyHillSphereOfInfluence[2] = 2430000 -- in meters
  globalSomCelestialBodyEquatorialSurfaceGravity[2] = 1.628 -- in m/s^2
  globalSomCelestialBodyEquatorialRadius[2] = 200000 -- in meters

  -- KSP_Galaxy1_Kerbol_Kerbin_Minmus per http://kspwiki.nexisonline.net/wiki/Minmus 20120813
  globalSomCelestialBodyName[3] = "KSP_Galaxy1_Kerbol_Kerbin_Minmus"
  globalSomCelestialBodyGM[3] = 2825505000 -- in m^3/s^2
  globalSomCelestialBodySiderealRotationalPeriod[3] = 1077379 -- in seconds (299.272 hours * 60 min/hr * 60 sec/min, and equal to Orbital Period)
  globalSomCelestialBodyMinStableAltitude[3] = 5900 -- in meters
  globalSomCelestialBodyHillSphereOfInfluence[3] = 2713000  -- in meters
  globalSomCelestialBodyEquatorialSurfaceGravity[3] = 1.628 -- in m/s^2
  globalSomCelestialBodyEquatorialRadius[3] = 60000 -- in meters

  -- RL_MilkyWay_Sol_Earth per https://en.wikipedia.org/wiki/Earth 20120813
  --   and per http://www.earthobservatory.nasa.gov/Features/OrbitsCatalog/ 20120813
  globalSomCelestialBodyName[4] = "RL_MilkyWay_Sol_Earth"
  globalSomCelestialBodyGM[4] = 398600441800000 -- in m^3/s^2
  globalSomCelestialBodySiderealRotationalPeriod[4] = 86164.098903691 -- sidereal rotation, stellar day, in seconds ((23 hours * 60 min/hr + 56 min) * 60 sec/min)+4.1 sec
  globalSomCelestialBodyMinStableAltitude[4] = 180000 -- in meters
  globalSomCelestialBodyHillSphereOfInfluence[4] = 1500000000  -- in meters
  globalSomCelestialBodyEquatorialSurfaceGravity[4] = 9.780327 -- in m/s^2
  globalSomCelestialBodyEquatorialRadius[4] = 6378100 -- in meters

  -- math.pi will do instead of a global, but if you like, try: globalPi = 3.1415926535897932384 -- a couple more digits than double precision can handle, to near-maximize accuracy
  globale = 2.7182818284590452353 -- a couple more digits than double precision can handle, to near-maximize accuracy
end

function globalSomMainBodyReset(MainBodyName)
  -- MainBodyName must be one of RL_MilkyWay_Sol_Earth, KSP_Galaxy1_Kerbol_Kerbin_Mun, KSP_Galaxy1_Kerbol_Kerbin, or KSP_Galaxy1_Kerbol_Kerbin_Minmus

  -- in case we don't find the name we're looking for, set to defaults.
  globalSomMainBodyNumber = nil
  globalSomMainBodyName = nil
  globalSomMainBodyGM = nil
  globalSomMainBodySiderealRotationalPeriod = nil
  globalSomMainBodyMinStableAltitude = nil
  globalSomMainBodyHillSphereOfInfluence = nil
  globalSomMainBodyEquatorialSurfaceGravity = nil
  globalSomMainBodyEquatorialRadius = nil


  local i
  -- Loop through however many known Celestial Body Names we have
  for i = 1, #globalSomCelestialBodyName do
    -- If the name of that celestial body matches the name we were given, let's use that entry's data!
    if MainBodyName == globalSomCelestialBodyName[i] then
      globalSomMainBodyNumber = i + 0 -- this looks strange, but simply assigning = i ends up with #globalSomCelestialBodyName - presumably it acts as pointer assignment instead
      globalSomMainBodyName = MainBodyName
      globalSomMainBodyGM = globalSomCelestialBodyGM[i] -- in km^3/s^2
      globalSomMainBodySiderealRotationalPeriod = globalSomCelestialBodySiderealRotationalPeriod[i] -- in seconds
      globalSomMainBodyMinStableAltitude = globalSomCelestialBodyMinStableAltitude[i] -- in meters
      globalSomMainBodyHillSphereOfInfluence = globalSomCelestialBodyHillSphereOfInfluence[i] -- in meters
      globalSomMainBodyEquatorialSurfaceGravity = globalSomCelestialBodyEquatorialSurfaceGravity[i] -- in m/s^2
      globalSomMainBodyEquatorialRadius = globalSomCelestialBodyEquatorialRadius[i] -- in meters
    end
  end

  if globalSomMainBodyNumber == nil then
    print "ERROR: Unknown MainBodyName in globalSomMainBodyReset"
    print "  expecting RL_MilkyWay_Sol_Earth, KSP_Galaxy1_Kerbol_Kerbin, KSP_Galaxy1_Kerbol_Kerbin_Mun, or KSP_Galaxy1_Kerbol_Kerbin_Minmus"
  end

end



function calcSomSemiMajorAxisOrbitingSphereFromApPeMbr(Ap, Pe, MainBodyRadius)
  -- Ap is the Apoapsis in meters
  -- Pe is the Periapsis in meters
  -- MainBodyRadius is the equatorial radius of the body being orbited in meters
  -- The semi-major axis (a) of an ellipse (or a circle) is half the longest dimension of the
  --   total ellipse; in other words, for a small body of trivial mass (your spaceship or satellite)
  --   compared to the body it orbits (the mainBody), it's the Apoapsis plus the body it orbit's radius
  --   (since it's assumed to orbit the center of the body), plus the same for Periapsis, divided by two.
  -- This is deliberately coded for maximum readability.
  return ((Ap + MainBodyRadius) + (Pe + MainBodyRadius))/2
end


function calcSomPeorapFromaAporpeMbr(a, Aporpe, MainBodyRadius)
  -- a is the semi-major axis in meters
  -- Aporpe is either Apoapsis or Periapsis, in meters - this returns the other one.
  -- MainBodyRadius is the equatorial radius of the body being orbited in meters
  -- This returns the opposite apsis from whatever you entered, in meters
  --   i.e. if you feed it Apoapsis, it tells you Periapsis
  --   i.e. if you feed it Periapsis, it tells you Apoapsis
  return (2*a-(Aporpe + MainBodyRadius*2))
end

function calcSomSiderealOrbitalPeriodFromaGm(a, GM)
  -- a is the semi-major axis in meters
  -- GM is the standard gravitational parameter in m^3/s^2 (also known as u)
  -- returns the sidereal orbital period in seconds
  -- For the equation, please see http://en.wikipedia.org/wiki/Orbital_period
  return (2 * math.pi * (a^3/GM)^0.5)
end

function calcSomSemiMajorAxisFromPeriodGm(Period, GM)
  -- Period in seconds
  -- GM is the standard gravitational parameter in m^3/s^2 (also known as u)
  -- returns the semi-major axis (a) in meters
  -- For the equation, please see http://en.wikipedia.org/wiki/Orbital_period and derive a = ...
  return ((GM*(Period/(2*math.pi))^2)^(1/3))
end

function calcSomSiderealOrbitalPeriodFromApPeMbrGm(Ap, Pe, MainBodyRadius, GM)
  -- relies on globalSomOverallInit
  -- Ap is the Apoapsis in meters
  -- Pe is the Periapsis in meters
  -- MainBodyRadius is the equatorial radius of the body being orbited in meters
  -- GM is the standard gravitational parameter in m^3/s^2
  -- returns the sidereal orbital period in seconds
  local semiMajorAxis = calcSomSemiMajorAxisOrbitingSphereFromApPeMbr(Ap, Pe, MainBodyRadius)
  return calcSomSiderealOrbitalPeriodFromaGm(semiMajorAxis, GM)
end


function calcSomOrbitalVelocityFromaGmRfcb(a, GM, RadiusFromCentralBody)
  -- a is the semi-major axis in meters
  -- GM is the standard gravitational parameter in m^3/s^2 (also known as u)
  -- RadiusFromCentralBody is the distance in meters from the orbiting object to
  --   the center of gravity it's orbiting (i.e. the center of the planet)
  -- returns the orbital velocity in m/s
  return ((GM*(2/RadiusFromCentralBody - 1/a))^0.5)
end


function calcSomOrbitalVelocityFromaGmAltaslMbr(a, GM, AltitudeASL, MainBodyRadius)
  -- a is the semi-major axis in meters
  -- GM is the standard gravitational parameter in m^3/s^2 (also known as u)
  -- AltitudeASL is the altitude above sea level in meters
  -- MainBodyRadius is the equatorial radius of the body being orbited in meters
  -- returns the orbital velocity in m/s
  local RadiusFromCentralBody = AltitudeASL + MainBodyRadius
  return calcSomOrbitalVelocityFromaGmRfcb(a, GM, RadiusFromCentralBody)
end


function calcSomGMFromaOvRfcb(a, OrbitalVelocity, RadiusFromCentralBody)
  -- a is the semi-major axis in meters
  -- OrbitalVelocity is the orbital velocity in m/s
  -- RadiusFromCentralBody is the distance in meters from the orbiting object to
  --   the center of gravity it's orbiting (i.e. the center of the planet)
  -- returns GM, the standard gravitational parameter in m^3/s^2 (also known as u)
  return (OrbitalVelocity^2 / (2/RadiusFromCentralBody - 1/a))
end


function calcSomGMFromaOvAltaslMbr(a, OrbitalVelocity, AltitudeASL, MainBodyRadius)
  -- a is the semi-major axis in meters
  -- OrbitalVelocity is the orbital velocity in m/s
  -- AltitudeASL is the altitude above sea level in meters
  -- MainBodyRadius is the equatorial radius of the body being orbited in meters
  -- returns GM, the standard gravitational parameter in m^3/s^2 (also known as u)
  local RadiusFromCentralBody = AltitudeASL + MainBodyRadius
  return calcSomGMFromaOvRfcb(a, OrbitalVelocity,RadiusFromCentralBody)
end


function guessSomWhatBodyIsTheMainBodyFromaOvAltaslMoe(a, OrbitalVelocity, AltitudeASL, MarginOfError)
  -- a is the semi-major axis in meters
  -- OrbitalVelocity is the orbital velocity in m/s
  -- AltitudeASL is the altitude above sea level in meters
  -- returns the SimpleOrbitalMechanics name of the body we think we're orbiting.
  --   suitable to pass to globalSomMainBodyReset, or UNKNOWN if it cannot tell
  -- NOTE: This goes through both RL and KSP bodies!
  -- NOTE: You actually have to be in orbit first.
  -- NOTE: Depends on globalSomOverallInit() having parameters for the body in question, and
  --   having been called before now.
  -- This is _extremely_ simplistic - we're just going to go through all the planets
  --   we know about, and check each one.  First one we find that matches closely enough wins!

  local i
  -- Loop through however many known Celestial Body Names we have
  for i = 1, #globalSomCelestialBodyName do
    -- If the name of that celestial body matches the name we were given, let's use that entry's data!
    if math.abs((calcSomGMFromaOvAltaslMbr(a, OrbitalVelocity, AltitudeASL, globalSomCelestialBodyEquatorialRadius[i]))/globalSomCelestialBodyGM[i] - 1.0) < MarginOfError then
      return globalSomCelestialBodyName[i]
    end
  end

  -- if we got here, it's because we didn't find a match and return with it from inside the for loop
  return "UNKNOWN"
end



function calcSomPeorapFromAporpeMbrGmSop(Aporpe, MainBodyRadius, GM, SiderealOrbitalPeriod)
  -- Aporpe is either Apoapsis or Periapsis, in meters - you don't even need to know which
  -- MainBodyRadius is the equatorial radius of the body being orbited in meters
  -- GM is the standard gravitational parameter in m^3/s^2
  -- SiderealOrbitalPeriod is the sidereal orbital period in seconds you wish the orbit to have
  -- returns the Peorap (if given Apoapsis, returns Periapsis; if given Periapsis, returns Apoapsis
  --   - and if you don't know which you're giving it going in, whichever is bigger at the end is Apoapsis)
  -- This is deliberately coded for maximum readability; even a genius Kerbal should be able to follow
  --   these instructions, if they take then slow and do them one small step at a time.

  -- NOTE: Molniya orbits circle twice per day, so the SiderealOrbitalPeriod needs to be half the SiderealRotationalPeriod
  -- NOTE: Tundra orbits circle once per day, so the SiderealOrbitalPeriod needs to be equal to the SiderealRotationalPeriod

  -- a is the semi-major axis, which we'll use to calculate the Peorap shortly.
  local a = calcSomSemiMajorAxisFromPeriodGm(SiderealOrbitalPeriod, GM)

  -- now that we have the semi-major axis, we can calculate the Peorap
  -- i.e. if we were given the Ap, let's find the Pe.  If we were given the Pe, find the Ap.
  local Peorap = calcSomPeorapFromaAporpeMbr(a, Aporpe, MainBodyRadius)

  -- WARNING WARNING WARNING - there's every chance this orbit will kill you or send you out of the local orbit!!!!
  --  Use isSomOrbitStableFromApPeHsoiMsa to check it before using!!!

  return Peorap
end

function calcSomMolniyaPeorapFromAporpeMbrGmSrp(Aporpe, MainBodyRadius, GM, SiderealRotationalPeriod)
  -- Aporpe is either Apoapsis or Periapsis, in meters - you don't even need to know which
  -- MainBodyRadius is the equatorial radius of the body being orbited in meters
  -- GM is the standard gravitational parameter in m^3/s^2
  -- SiderealRotationalPeriod is the sidereal rotational period in seconds of the body being orbited
  -- returns the Peorap (if given Apoapsis, returns Periapsis; if given Periapsis, returns Apoapsis
  --   - and if you don't know which you're giving it going in, whichever is bigger at the end is Apoapsis)
  --   for a Molniya orbit, i.e. a semisynchronous orbit (two orbits per day)
  -- NOTE: See calcSomPeorapFromAporpeMbrGmSop for details

  -- NOTE: Molniya orbits circle twice per day, so the SiderealOrbitalPeriod needs to be half the SiderealRotationalPeriod
  local molniyaSiderealOrbitalPeriod = SiderealRotationalPeriod / 2
  return calcSomPeorapFromAporpeMbrGmSop(Aporpe, MainBodyRadius, GM, molniyaSiderealOrbitalPeriod)

  -- WARNING WARNING WARNING - there's every chance this orbit will kill you or send you out of the local orbit!!!!
  --  Use isSomOrbitStableFromApPeHsoiMsa to check it before using!!!

end

function calcSomApPeFromAporpePeorap(Aporpe, Peorap)
  -- Aporpe is either Apoapsis or Periapsis, in meters - you don't even need to know which
  -- Peorap is either Periapsis or Apoapsis, in meters - you don't even need to know which
  -- returns Ap, Pe
  --[[ Ex:
       do
       local localAp, localPe = calcSomApPeFromAporpePeorap(75000,3095000)
       print("Ap: " .. localAp .. " Pe: " .. localPe)
       end
       ]]

  -- We'll find out which is which because the bigger one is the Ap.  If they're equal (perfectly circular orbit), who cares!
  return math.max(Aporpe, Peorap), math.min(Aporpe, Peorap)
end



function calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe(HillSphereOfInfluence, MinStableAltitude, MainBodyRadius, GM, SiderealRotationalPeriod, MarginOfError)
  -- HillSphereOfInfluence is the range in meters at which the gravity of the body is effective
  -- MinStableAltitude is the altitude above sea level in meters at which orbiting is safe and stable
  --   i.e. it's above both the atmosphere's effects, and the tallest mountains.
  -- MainBodyRadius is the equatorial radius of the body being orbited in meters
  -- GM is the standard gravitational parameter in m^3/s^2
  -- SiderealRotationalPeriod is the sidereal rotational period in seconds of the body being orbited
  -- MarginOfError 0.001 means allow 0.1% for a "safety zone" above the minimum stable altitude, or below the Hill sphere of influence.
  -- Returns Ap and Pe of the best possible Molniya orbit for a given body
  --   Returns -1, -1 if no stable Molniya orbit is possible
  --[[ ex.
    do
    local localAp, localPe = calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe(1500000000, 180000, 6378100, 398600441800000, 86164.098903691, 1.777777777)
    print("Ap " .. localAp .. " Pe " .. localPe)
    end
    ]]
  -- For oblate celestial bodies, of course, the inclination of the actual orbit must be 63.4 or 116.6 degrees, so the effect of the oblate equitorial belt is nullified.

  -- First, let's see if we can get a stable Molniya orbit with MinStableAltitude (including MarginOfError) as the Periapsis
  --   this would be the ideal Molniya orbit, with maximum Apoapsis dwell time.
  local Ap = calcSomMolniyaPeorapFromAporpeMbrGmSrp(MinStableAltitude * (1 + MarginOfError), MainBodyRadius, GM, SiderealRotationalPeriod)
  if isSomOrbitStableFromApPeHsoiMsa(Ap, MinStableAltitude  * (1 + MarginOfError), HillSphereOfInfluence, MinStableAltitude) == 0 then
   -- isSomOrbitStableFromApPeHsoiMsa returns 0 for a stable orbit.
   return Ap, MinStableAltitude * (1 + MarginOfError)
  end

  -- Well, if we're here, it wasn't a stable orbit.  We've got one more shot; let's see if we can get a crappy Molniya orbit
  --   Perhaps if the Apoapsis is the Hill sphere of influence (minus a 1000 meter buffer), we can get a periapsis that works
  local Pe = calcSomMolniyaPeorapFromAporpeMbrGmSrp(HillSphereOfInfluence  * (1 - MarginOfError), MainBodyRadius, GM, SiderealRotationalPeriod)
  if isSomOrbitStableFromApPeHsoiMsa(HillSphereOfInfluence  * (1 - MarginOfError), Pe, HillSphereOfInfluence, MinStableAltitude) == 0 then
    return HillSphereOfInfluence  * (1 - MarginOfError), Pe
  end

  -- And if we made it past both tries, it's just not going to happen.
  return -1, -1

end


function isSomOrbitStableFromApPeHsoiMsa(Ap, Pe, HillSphereOfInfluence, MinStableAltitude)
  -- Ap is the Apoapsis in meters
  -- Pe is the Periapsis in meters
  -- HillSphereOfInfluence is the range in meters at which the gravity of the body is effective
  -- MinStableAltitude is the altitude above sea level in meters at which orbiting is safe and stable
  --   i.e. it's above both the atmosphere's effects, and the tallest mountains.


  -- technically, this is simple bit flags.
  -- 0 = stable
  -- 1 = burn up in atmosphere
  -- 2 = escape SOI
  -- 3 = both burn up and escape SOI (whichever comes first)
  -- 4 = your Apoapsis is greater than your Periapsis, which is just wrong.
  returnValueisSomOrbitStableFromApPeHsoiMsa = 0 -- this should be a local, but when set to local scope the additions inside the ifs fail completely.
  if Pe < MinStableAltitude then
    returnValueisSomOrbitStableFromApPeHsoiMsa = returnValueisSomOrbitStableFromApPeHsoiMsa + 1
  end
  if Ap > HillSphereOfInfluence then
    returnValueisSomOrbitStableFromApPeHsoiMsa = returnValueisSomOrbitStableFromApPeHsoiMsa + 2
  end
  if Pe > Ap then
    returnValueisSomOrbitStableFromApPeHsoiMsa = returnValueisSomOrbitStableFromApPeHsoiMsa + 4
  end
  return returnValueisSomOrbitStableFromApPeHsoiMsa

end


function testSomTestHarness()
  -- If you get messages, then something broke!

  assert(calcSomSemiMajorAxisOrbitingSphereFromApPeMbr(4501000, 72001, 600000)>=2886500,"calcSomSemiMajorAxisOrbitingSphereFromApPeMbr test 1 Eccentric Kerbin orbit failed")
  assert(calcSomSemiMajorAxisOrbitingSphereFromApPeMbr(4501000, 72001, 600000)<=2886501,"calcSomSemiMajorAxisOrbitingSphereFromApPeMbr test 2 Eccentric Kerbin orbit failed")

  assert(calcSomPeorapFromaAporpeMbr(750000, 75000, 600000)==225000,"calcSomPeorapFromaAporpeMbr test 1 Silly Kerbin orbit failed")

  assert(calcSomSiderealOrbitalPeriodFromaGm(26378000,398600441800000)>=42633,"calcSomSiderealOrbitalPeriodFromaGm test 1 GPS Earth orbit failed")
  assert(calcSomSiderealOrbitalPeriodFromaGm(26378000,398600441800000)<=42637,"calcSomSiderealOrbitalPeriodFromaGm test 2 GPS Earth orbit failed")

  assert(calcSomSemiMajorAxisFromPeriodGm(86164,398600441800000)>=42164138,"calcSomSemiMajorAxisFromPeriodGm test 1 Geosync Earth orbit failed")
  assert(calcSomSemiMajorAxisFromPeriodGm(86164,398600441800000)<=42164142,"calcSomSemiMajorAxisFromPeriodGm test 2 Geosync Earth orbit failed")

  assert(calcSomSiderealOrbitalPeriodFromApPeMbrGm(20314000, 20047000, 6378100, 398600441800000)>=43072,"calcSomSiderealOrbitalPeriodFromApPeMbrGm test 1 GPS Earth orbit failed")
  assert(calcSomSiderealOrbitalPeriodFromApPeMbrGm(20314000, 20047000, 6378100, 398600441800000)<=43076,"calcSomSiderealOrbitalPeriodFromApPeMbrGm test 2 GPS Earth orbit failed")

  assert(calcSomOrbitalVelocityFromaGmRfcb(750000,3530461000000,750000)>=2169,"calcSomOrbitalVelocityFromaGmRfcb test 1 150km Kerbin circular orbit failed.")
  assert(calcSomOrbitalVelocityFromaGmRfcb(750000,3530461000000,750000)<=2170,"calcSomOrbitalVelocityFromaGmRfcb test 2 150km Kerbin circular orbit failed.")

  assert(calcSomOrbitalVelocityFromaGmAltaslMbr(725000,3530461000000,125000,600000)>=2206,"calcSomOrbitalVelocityFromaGmAltaslMbr test 1 125km Kerbin circular orbit failed.")
  assert(calcSomOrbitalVelocityFromaGmAltaslMbr(725000,3530461000000,125000,600000)<=2207,"calcSomOrbitalVelocityFromaGmAltaslMbr test 2 125km Kerbin circular orbit failed.")

  assert(calcSomGMFromaOvRfcb(362941.275006862,383.306541880459,399135.253432537)>=65138373282,"calcSomGMFromaOvRfcb test 1 126610-199271m Kerbin Mun orbit failed.")
  assert(calcSomGMFromaOvRfcb(362941.275006862,383.306541880459,399135.253432537)<=65138373302,"calcSomGMFromaOvRfcb test 2 126610-199271m Kerbin Mun orbit failed.")
  assert(calcSomGMFromaOvAltaslMbr(468746.649196013,407.483188333337,227130.402050318,200000)>=65138680733,"calcSomGMFromaOvAltaslMbr test 1 197405-340087km Kerbin Mun orbit failed.")
  assert(calcSomGMFromaOvAltaslMbr(468746.649196013,407.483188333337,227130.402050318,200000)<=65138680753,"calcSomGMFromaOvAltaslMbr test 2 197405-340087km Kerbin Mun orbit failed.")

  assert(guessSomWhatBodyIsTheMainBodyFromaOvAltaslMoe(468780,394.1,242942,0.01)=="KSP_Galaxy1_Kerbol_Kerbin_Mun","guessSomWhatBodyIsTheMainBodyFromaOvAltaslMoe test 1 Kerbin Mun failed.")
  assert(guessSomWhatBodyIsTheMainBodyFromaOvAltaslMoe(65955,207.0,5950,0.01)=="KSP_Galaxy1_Kerbol_Kerbin_Minmus","guessSomWhatBodyIsTheMainBodyFromaOvAltaslMoe test 2 Kerbin Minmus failed.")
  assert(guessSomWhatBodyIsTheMainBodyFromaOvAltaslMoe(2187000,1130.4,1839000,0.01)=="KSP_Galaxy1_Kerbol_Kerbin","guessSomWhatBodyIsTheMainBodyFromaOvAltaslMoe test 3 Kerbin failed.")
  assert(guessSomWhatBodyIsTheMainBodyFromaOvAltaslMoe(42164138,3075,35786000,0.01)=="RL_MilkyWay_Sol_Earth","guessSomWhatBodyIsTheMainBodyFromaOvAltaslMoe test 4 Earth Geosync failed.")
  assert(guessSomWhatBodyIsTheMainBodyFromaOvAltaslMoe(5,5,5,0.0001)=="UNKNOWN","guessSomWhatBodyIsTheMainBodyFromaOvAltaslMoe test 5 Unknown orbit failed.")

  assert(calcSomPeorapFromAporpeMbrGmSop(703000, 6378100, 398600441800000, 5928)>=699653,"calcSomMolniyaPeorapFromAporpeMbrGmSrp test 1 Classic Earth Landsat failed")
  assert(calcSomPeorapFromAporpeMbrGmSop(703000, 6378100, 398600441800000, 5928)<=699657,"calcSomMolniyaPeorapFromAporpeMbrGmSrp test 2 Classic Earth Landsat failed")
  assert(calcSomMolniyaPeorapFromAporpeMbrGmSrp(500000, 6378100, 398600441800000, 86164.098903691)>=39867327,"calcSomMolniyaPeorapFromAporpeMbrGmSrp test 1 Classic Earth Molniya failed")
  assert(calcSomMolniyaPeorapFromAporpeMbrGmSrp(500000, 6378100, 398600441800000, 86164.098903691)<=39867330,"calcSomMolniyaPeorapFromAporpeMbrGmSrp test 2 Classic Earth Molniya failed")
  assert(calcSomMolniyaPeorapFromAporpeMbrGmSrp(39867328.3148597, 6378100, 398600441800000, 86164.098903691)>=499998,"calcSomMolniyaPeorapFromAporpeMbrGmSrp test 3 Classic Earth Molniya failed")
  assert(calcSomMolniyaPeorapFromAporpeMbrGmSrp(39867328.3148597, 6378100, 398600441800000, 86164.098903691)<=500002,"calcSomMolniyaPeorapFromAporpeMbrGmSrp test 4 Classic Earth Molniya failed")
  do
    local localAp, localPe = calcSomApPeFromAporpePeorap(75000,3095000)
    assert(localAp==3095000,"calcSomApPeFromAporpePeorap test 1 Pe first failed")
    assert(localPe==75000,"calcSomApPeFromAporpePeorap test 2 Pe first failed")
    local localAp, localPe = calcSomApPeFromAporpePeorap(2,1)
    assert(localAp==2,"calcSomApPeFromAporpePeorap test 3 Ap first failed")
    assert(localPe==1,"calcSomApPeFromAporpePeorap test 4 Ap first failed")
  end

  do
    local localAp, localPe = calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe(1500000000, 180000, 6378100, 398600441800000, 86164.098903691, 1.777777777)
    assert(localAp>=39867318,"calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe test 1 Earth classic Molniya failed.")
    assert(localAp<=39867338,"calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe test 2 Earth classic Molniya failed.")
    assert(localPe>=499998,"calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe test 3 Earth classic Molniya failed.")
    assert(localPe<=500002,"calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe test 4 Earth classic Molniya failed.")
    local localAp, localPe = calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe(2430000, 4700, 200000, 65136000000, 147600, 0.01)
    assert(localAp>=2405690,"calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe test 5 Kerbin Mun crappy Molniya failed.")
    assert(localAp<=2405710,"calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe test 6 Kerbin Mun crappy Molniya failed.")
    assert(localPe>=1352323,"calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe test 7 Kerbin Mun crappy Molniya failed.")
    assert(localPe<=1352343,"calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe test 8 Kerbin Mun crappy Molniya failed.")
    local localAp, localPe = calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe(2130000, 4700, 200000, 65136000000, 197600, 0.01)
    assert(localAp == (-1),"calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe test 9 Impossible Molniya failed.") -- parenthesis around the -1 is required to get the desired  error
    assert(localPe == (-1),"calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe test 10 Impossible Molniya failed.") -- parenthesis around the -1 is required to get the desired  error
  end


  assert(isSomOrbitStableFromApPeHsoiMsa(100000,10000,200000,5000)==0,"isSomOrbitStableFromApPeHsoiMsa test 1 (stable orbit) failed")
  assert(isSomOrbitStableFromApPeHsoiMsa(100000,10000,200000,10001)==1,"isSomOrbitStableFromApPeHsoiMsa test 2 (areobrake orbit) failed")
  assert(isSomOrbitStableFromApPeHsoiMsa(200001,10000,200000,5000)==2,"isSomOrbitStableFromApPeHsoiMsa test 3 (escape orbit) failed")
  assert(isSomOrbitStableFromApPeHsoiMsa(200001,10000,200000,10001)==3,"isSomOrbitStableFromApPeHsoiMsa test 4 (deathtrap orbit) failed")
  assert(isSomOrbitStableFromApPeHsoiMsa(20000,30000,200000,5000)==4,"isSomOrbitStableFromApPeHsoiMsa test 5 (stupid orbit) failed")
  print "Tests complete."
end



function helpSomSimpleOrbitalMechanics1()
  print "Usage: globalSomOverallInit()"
  print "   Initializes unchanging global variables"
  print "Usage: globalSomMainBodyReset(MainBodyName)"
  print "   Initialized global variables that change based on the body"
  print "   that is being orbited."
  print "   MainBodyName must be one of RL_MilkyWay_Sol_Earth, KSP_Galaxy1_Kerbol_Kerbin_Mun, KSP_Galaxy1_Kerbol_Kerbin, or KSP_Galaxy1_Kerbol_Kerbin_Minmus"
  print "Usage: calcSomSemiMajorAxisOrbitingSphereFromApPeMbr(Ap, Pe, MainBodyRadius)"
  print "   Returns the semi-major axis based on Ap the Apoapsis,"
  print "   Pe the Periapsis, and MainBodyRadius, the equatorial radius"
  print "   of the body being orbited (assumed to be a sphere; at this time"
  print "   the effects of a highly oblate body with a high inclination"
  print "   orbit have not been investigated."
  print "   All units must be the same.  I.e. for a 50x125km orbit around"
  print "   Kerbin, calcSomSemiMajorAxisOrbitingSphereFromApPeMbr(125000,50000,600000)"
  print "Usage: calcSomPeorapFromaAporpeMbr(a, Aporpe, MainBodyRadius)"
  print "   a is the semi-major axis in meters"
  print "   Aporpe is either Apoapsis or Periapsis, in meters"
  print "   MainBodyRadius is the equatorial radius of the body being orbited in meters"
  print "   This returns the opposite apsis from whatever you entered, in meters"
  print "   i.e. if you feed it Apoapsis, it tells you Periapsis"
  print "   i.e. if you feed it Periapsis, it tells you Apoapsis"
  print "Usage: calcSomSiderealOrbitalPeriodFromaGm(a, GM)"
  print "   Returns the period in seconds of the orbit, given a,"
  print "   the semi-major axis in meters, and GM, the gravitational"
  print "   parameter in m^3/s^2 of the body being orbited."
  print "Usage: calcSomSemiMajorAxisFromPeriodGm(Period, GM)"
  print "  Period in seconds"
  print "  GM is the standard gravitational parameter in m^3/s^2"
  print "  returns the semi-major axis (a) in meters"
end

function helpSomSimpleOrbitalMechanics2()
  print "Usage: calcSomSiderealOrbitalPeriodFromApPeMbrGm(Ap, Pe, MainBodyRadius, GM)"
  print "   Returns the period in seconds of the orbit, given Ap"
  print "   the Apoapsis, Pe the Periapsis, and GM the gravitational"
  print "   parameter of the body being orbited."
  print "Usage: calcSomOrbitalVelocityFromaGmRfcb(a, GM, RadiusFromCentralBody)"
  print "   a is the semi-major axis in meters"
  print "   GM is the standard gravitational parameter in m^3/s^2 (also known as u)"
  print "   RadiusFromCentralBody is the distance in meters from the orbiting object"
  print "     to the center of gravity it's orbiting (i.e. the center of the planet)"
  print "   returns the orbital velocity in m/s"
  print "Usage: calcSomOrbitalVelocityFromaGmAltaslMbr(a, GM, AltitudeASL, MainBodyRadius)"
  print "   a is the semi-major axis in meters"
  print "   GM is the standard gravitational parameter in m^3/s^2 (also known as u)"
  print "   AltitudeASL is the altitude above sea level in meters"
  print "   MainBodyRadius is the equatorial radius of the body being orbited in meters"
  print "Usage: calcSomGMFromaOvRfcb(a, OrbitalVelocity, RadiusFromCentralBody)"
  print "   a is the semi-major axis in meters"
  print "   OrbitalVelocity is the orbital velocity in m/s"
  print "   RadiusFromCentralBody is the distance in meters from the orbiting object to"
  print "     the center of gravity it's orbiting (i.e. the center of the planet)"
  print "   returns GM, the standard gravitational parameter in m^3/s^2 (also known as u)"
end

function helpSomSimpleOrbitalMechanics3()
  print "Usage: calcSomGMFromaOvAltaslMbr(a, OrbitalVelocity, AltitudeASL, MainBodyRadius)"
  print "   a is the semi-major axis in meters"
  print "   OrbitalVelocity is the orbital velocity in m/s"
  print "   AltitudeASL is the altitude above sea level in meters"
  print "   MainBodyRadius is the equatorial radius of the body being orbited in meters"
  print "   returns GM, the standard gravitational parameter in m^3/s^2 (also known as u)"
  print "Usage: guessSomWhatBodyIsTheMainBodyFromaOvAltaslMoe(a, OrbitalVelocity, AltitudeASL, MarginOfError)"
  print "   a is the semi-major axis in meters"
  print "   OrbitalVelocity is the orbital velocity in m/s"
  print "   AltitudeASL is the altitude above sea level in meters"
  print "   returns the SimpleOrbitalMechanics name of the body we think we're orbiting."
  print "     suitable to pass to globalSomMainBodyReset, or UNKNOWN if it cannot tell"
  print "Usage: calcSomPeorapFromAporpeMbrGmSop(Aporpe, MainBodyRadius, GM, SiderealOrbitalPeriod)"
  print "   Aporpe is either Apoapsis or Periapsis, in meters - you don't even need to know which"
  print "   MainBodyRadius is the equatorial radius of the body being orbited in meters"
  print "   GM is the standard gravitational parameter in m^3/s^2"
  print "   SiderealOrbitalPeriod is the sidereal orbital period in seconds you wish the orbit to have"
  print "   returns the Peorap (if given Apoapsis, returns Periapsis; if given Periapsis, returns Apoapsis"
  print "      and if you don't know which you're giving it going in, whichever is bigger at the end is Apoapsis)"
  print "   WARNING WARNING WARNING - there's every chance this orbit will kill you or send you out of the local orbit!!!!"
  print "      Use isSomOrbitStableFromApPeHsoiMsa to check your orbit before using!!!"
end

function helpSomSimpleOrbitalMechanics4()
  print "Usage: calcSomMolniyaPeorapFromAporpeMbrGmSrp(Aporpe, MainBodyRadius, GM, SiderealRotationalPeriod)"
  print "   Aporpe is either Apoapsis or Periapsis, in meters - you don't even need to know which"
  print "   MainBodyRadius is the equatorial radius of the body being orbited in meters"
  print "   GM is the standard gravitational parameter in m^3/s^2"
  print "   SiderealRotationalPeriod is the sidereal rotational period in seconds of the body being orbited"
  print "   returns the Peorap (if given Apoapsis, returns Periapsis; if given Periapsis, returns Apoapsis"
  print "      and if you don't know which you're giving it going in, whichever is bigger at the end is Apoapsis)"
  print "      for a Molniya orbit, i.e. a semisynchronous orbit (two orbits per day)"
  print "   WARNING WARNING WARNING - there's every chance this orbit will kill you or send you out of the local orbit!!!!"
  print "      Use isSomOrbitStableFromApPeHsoiMsa to check your orbit before using!!!"
  print "Usage: calcSomApPeFromAporpePeorap(Aporpe, Peorap)"
  print "   Aporpe is either Apoapsis or Periapsis, in meters - you don't even need to know which"
  print "   Peorap is either Periapsis or Apoapsis, in meters - you don't even need to know which"
  print "   returns Ap, Pe"
  print "   Ex: "
  print "       local localAp, localPe = calcSomApPeFromAporpePeorap(75000,3095000)"
  print "       print(\"Ap: \" .. localAp .. \" Pe: \" .. localPe)"
  print "Usage: testSomTestHarness()"
  print "   Just run it - if you get a message, something's broken.  Read this code for usage guidance!!!"
  print "Usage: isSomOrbitStableFromApPeHsoiMsa(Ap, Pe, HillSphereOfInfluence, MinStableAltitude)"
  print "   Returns 0 for stable, 1 for burn up in atmosphere, 2 for escape"
  print "   sphere of influence (SOI), and 3 for burn up AND escape SOI"
  print "   All values must be in the same units, including Ap Apoapsis,"
  print "   Pe Periapsis, HillSphereOfInfluence the Hill sphere of"
  print "   gravitational influence, and MinStableAltitude, the minimum"
  print "   altitude to be out of the atmosphere and above the highest"
  print "   mountain."
end

function helpSomSimpleOrbitalMechanics5()
  print "Usage: calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe(HillSphereOfInfluence, MinStableAltitude, MainBodyRadius, GM, SiderealRotationalPeriod, MarginOfError)"
  print "   HillSphereOfInfluence is the range in meters at which the gravity of the body is effective"
  print "   MinStableAltitude is the altitude above sea level in meters at which orbiting is safe and stable"
  print "     i.e. it's above both the atmosphere's effects, and the tallest mountains."
  print "   MainBodyRadius is the equatorial radius of the body being orbited in meters"
  print "   GM is the standard gravitational parameter in m^3/s^2"
  print "   SiderealRotationalPeriod is the sidereal rotational period in seconds of the body being orbited"
  print "   MarginOfError 0.001 means allow 0.1% for a \"safety zone\" above the minimum stable altitude, or below the Hill sphere of influence."
  print "   Returns Ap and Pe of the best possible Molniya orbit for a given body"
  print "     Returns -1, -1 if no stable Molniya orbit is possible"
  print "   ex."
  print "       do"
  print "       local localAp, localPe = calcSomMolniyaApPeFromHsoiMsaMbrGmSrpMoe(1500000000, 180000, 6378100, 398600441800000, 86164.098903691, 1.777777777)"
  print "       print(\"Ap \" .. localAp .. \" Pe \" .. localPe)"
  print "       end"
  print "   For oblate celestial bodies, of course, the inclination of the actual orbit must be 63.4 or 116.6 degrees, so the effect"
  print "     of the oblate equitorial belt is nullified."
end

globalSomOverallInit() -- initialize required variables up front

print "Usage: helpSomSimpleOrbitalMechanics1() through helpSomSimpleOrbitalMechanics5()"
print "   gives detailed help on all Som functions."
print "   Som functions relate to orbital mechanics as a whole, in real life and in KSP."