{-|
Module      : OrbitalMechanics
Description : This module contains basic orbital mechanic functions.
License     : CC0
Maintainer  : frosch03@frosch03.de
Stability   : experimental
-}
module System.KSP.OrbitalMechanics
where

import System.KSP.DataConstructors


-- | 'Radius' of 'Orbit' at the current Position
type Radius = Double 

-- | 'Speed' is an alias for 'Double'
type Speed  = Double


-- | 'var_G' [Nm^2/kg^2] is the Gravitation constant in newton meter
-- squared over kilo gramms squared
var_G :: GravConst
var_G = 6.674e-11


-- | 'semiMajor' calculates the semi major axis of an orbit.
semiMajor   :: Orbit Body -> Double
semiMajor o
    = 0.5 * ((apoapsis o) + (periapsis o) + (2 * (r . celestial . centerBody $ o)))


-- | 'v' takes an orbit and a radius (from the center of the
-- centerBody) and calculates the orbital speed at that position.
v :: Orbit Body -> Radius -> Speed
v o radius =  sqrt $ mue * ( (2/radius) - (1/(semiMajor o)) ) 
    where
      mue = var_G * (m . celestial . centerBody $ o)


-- | 'v_e' calculates the escape velocity of that body.
v_e :: Body -> Speed
v_e b =  sqrt $ 2 * mue / (r . celestial $ b) 
    where
      mue = var_G * (m . celestial $ b)

      
-- | 'burnProgradeFromCircOrb' calculates the transformed orbit, after
-- the supplied delta V is applied to the initial orbit.
-- 
-- A prograde burn is done through a positive 'Speed' parameter, a
-- retrograde burn respective through a negative 'Speed'.
burnFromCircOrb :: Orbit Body -> Speed -> Orbit Body
burnFromCircOrb o dV
    = updateOrbit o (periapsis o) ((2 * a_new) - (2 * r_body) - (periapsis o))
    where
      a_new  = 1/((2/(semiMajor o))-(((v_circ+dV)^2)/mue))
      r_body = (r . celestial . centerBody $ o)
      m_body = (m . celestial . centerBody $ o)
      v_circ = v o (semiMajor o)
      mue    = var_G * m_body

-- | 'updateOrbit' is a function that takes an orbit and two
-- heights. It updates the apoapsis with the bigger height and the
-- periapsis with the smaller.
updateOrbit :: Orbit Body -> Double -> Double -> Orbit Body
updateOrbit o h1 h2
    | h1 > h2
    = o { apoapsis = h1, periapsis = h2 }

    | otherwise
    = o { apoapsis = h2, periapsis = h1 }

-- | 'burnAt' calculates the new orbit after a burn of 'Speed' delta V
-- is applied to the given orbit.
-- 
--     * 'f' is the function that calculates returns the distance of
--     the body within orbit, to the center body. Typical this is one
--     of 'apoapsis' or 'periapsis'
-- 
--     * 'o' is the initial orbit
-- 
--     * 'dV' is the amount of delta V to apply to the orbit. 
burnAt :: (Orbit Body -> Double) ->  Orbit Body -> Speed -> Orbit Body               
burnAt f o dV
    = updateOrbit o (f o) ((2 * a_new) - (2 * r_body) - (f o))
    where
      a_new     = 1/( (2/dist_body) - (((v_atPos+dV)^2)/mue) )
      dist_body = (f o) + r_body
      r_body    = (r . celestial . centerBody $ o)
      m_body    = (m . celestial . centerBody $ o)
      v_atPos   = v o dist_body
      mue       = var_G * m_body


-- | 'burnAtPeriapsis' calculates the new orbit after a burn of
-- 'Speed' delta V is applied to the given orbit at the periapsis.
burnAtPeriapsis :: Orbit Body -> Speed -> Orbit Body
burnAtPeriapsis = burnAt periapsis

-- | 'burnAt' calculates the new orbit after a burn of 'Speed' delta V
-- is applied to the given orbit at the apoapsis.
burnAtApoapsis :: Orbit Body -> Speed -> Orbit Body
burnAtApoapsis  = burnAt apoapsis


-- | 'hohmann' takes two orbits around the same centerBody. It
-- calculates both (v1 and v2 ) delta V changes for a hohmann
-- transfair.
hohmann :: Orbit Body -> Orbit Body -> (Double, Double)
hohmann o1 o2
    |   centerBody o1 /= centerBody o2
      || o1 == o2
    = (0,0)
hohmann o1 o2
    = ( sqrt(mue/r1) * (sqrt((2*r2)/(r1+r2)) - 1) 
      , sqrt(mue/r2) * (1 - (sqrt((2*r1)/(r1+r2))))
      )
    where 
      mue = var_G * (m . celestial . centerBody $ o1)
      r1 = (r . celestial . centerBody $ o1) + (apoapsis o1)
      r2 = (r . celestial . centerBody $ o2) + (apoapsis o2)