{-# language QuasiQuotes #-}

-- | Types and functions for dealing with Kepler orbits.
module Physics.Orbit
  ( -- * The Orbit data type and dependencies
    Orbit(..)
  , InclinationSpecifier(..)
  , PeriapsisSpecifier(..)
  , Classification(..)

    -- * Functions for dealing with orbits
    -- ** Utilities
  , isValid
  , classify
  , normalizeOrbit
    -- ** Orbital elements
  , apoapsis
  , meanMotion
  , period
  , arealVelocity
    -- *** Geometry
  , semiMajorAxis
  , semiMinorAxis
  , semiLatusRectum
  , hyperbolicApproachAngle
  , hyperbolicDepartureAngle
    -- ** Conversions

    -- *** To time since periapse
  , timeAtMeanAnomaly
  , timeAtEccentricAnomaly
  , timeAtHyperbolicAnomaly
  , timeAtTrueAnomaly

    -- *** To mean anomaly
  , meanAnomalyAtTime
  , meanAnomalyAtEccentricAnomaly
  , meanAnomalyAtHyperbolicAnomaly
  , meanAnomalyAtTrueAnomaly

    -- *** To eccentric anomaly
  , eccentricAnomalyAtTime
  , eccentricAnomalyAtMeanAnomaly
  , eccentricAnomalyAtMeanAnomalyFloat
  , eccentricAnomalyAtTrueAnomaly

    -- *** To hyperbolic anomaly
  , hyperbolicAnomalyAtTime
  , hyperbolicAnomalyAtMeanAnomaly
  , hyperbolicAnomalyAtMeanAnomalyDouble
  , hyperbolicAnomalyAtTrueAnomaly

    -- *** To true anomaly
  , trueAnomalyAtTime
  , trueAnomalyAtMeanAnomaly
  , trueAnomalyAtEccentricAnomaly
  , trueAnomalyAtHyperbolicAnomaly

    -- *** Properties of orbits
  , specificAngularMomentum
  , specificOrbitalEnergy
  , specificPotentialEnergyAtTrueAnomaly
  , specificKineticEnergyAtTrueAnomaly
  , speedAtTrueAnomaly
  , radiusAtTrueAnomaly

    -- *** Other utilities
  , escapeVelocityAtDistance

    -- * Unit synonyms
  , Quantity
  , Time
  , Distance
  , Speed
  , Mass
  , Angle
  , AngleH
  , RadianHyperbolic(..)
  , PlaneAngleHyperbolic(..)
  , Unitless

    -- * Reexported from 'Data.CReal'
  , Converge
  ) where

import           Control.Monad                  ( (<=<) )
import           Data.Bifunctor                 ( bimap
                                                , second
                                                )
import           Data.CReal.Converge            ( Converge
                                                , convergeErr
                                                )
import           Data.Constants.Mechanics.Extra
import           Data.Maybe                     ( fromJust )
import           Data.Metrology
import           Data.Metrology.Extra
import           Data.Metrology.Show            ( )
import           Data.Metrology.Unsafe          ( Qu(..)
                                                , UnsafeQu(..)
                                                )
import           Data.Units.SI.Parser
import           Numeric.AD                     ( Mode
                                                , Scalar
                                                , auto
                                                )
import           Numeric.AD.Halley              ( findZero
                                                , findZeroNoEq
                                                )
import           Numeric.AD.Internal.Identity   ( Id(..) )
import qualified Numeric.AD.Newton.Double      as Newton
import           Physics.Orbit.Metrology

--------------------------------------------------------------------------------
-- Types
--------------------------------------------------------------------------------

-- | Data type defining an orbit parameterized by the type used to
-- represent values
data Orbit a = Orbit { -- | The orbit's eccentricity, e.
                       --
                       -- 'eccentricity' must be non-negative.
                       --
                       -- An eccentricity of 0 describes a circular orbit.
                       --
                       -- An eccentricity of less than 1 describes an elliptic
                       -- orbit.
                       --
                       -- An eccentricity equal to 1 describes a parabolic orbit.
                       --
                       -- An eccentricity greater than 1 describes a hyperbolic
                       -- orbit.
                       Orbit a -> Unitless a
eccentricity                  :: !(Unitless a)
                       -- | The orbit's periapsis, q.
                       --
                       -- 'periapsis' must be positive.
                       --
                       -- The periapsis is the distance between the bodies at
                       -- their closest approach.
                     , Orbit a -> Distance a
periapsis                     :: !(Distance a)
                       -- | The 'inclinationSpecifier' describes the angle
                       -- between the obtital plane and the reference plane.
                     , Orbit a -> InclinationSpecifier a
inclinationSpecifier          :: !(InclinationSpecifier a)
                       -- | 'periapsisSpecifier' is 'Circular' iff
                       -- 'eccentricity' is 0
                       --
                       -- The periapsis specifier describes any rotation of
                       -- the orbit relative to the reference direction in the
                       -- orbital plane.
                     , Orbit a -> PeriapsisSpecifier a
periapsisSpecifier            :: !(PeriapsisSpecifier a)
                       -- | The gravitational parameter of the system's
                       -- primary, μ.
                       --
                       -- μ is equal to the mass of the primary times
                       -- <https://en.wikipedia.org/wiki/Gravitational_constant
                       -- G>.
                       --
                       -- 'primaryGravitationalParameter' must be positive.
                     , Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
primaryGravitationalParameter :: !(Quantity [si| m^3 s^-2 |] a)
                     }
  deriving (Int -> Orbit a -> ShowS
[Orbit a] -> ShowS
Orbit a -> String
(Int -> Orbit a -> ShowS)
-> (Orbit a -> String) -> ([Orbit a] -> ShowS) -> Show (Orbit a)
forall a. Show a => Int -> Orbit a -> ShowS
forall a. Show a => [Orbit a] -> ShowS
forall a. Show a => Orbit a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Orbit a] -> ShowS
$cshowList :: forall a. Show a => [Orbit a] -> ShowS
show :: Orbit a -> String
$cshow :: forall a. Show a => Orbit a -> String
showsPrec :: Int -> Orbit a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> Orbit a -> ShowS
Show, Orbit a -> Orbit a -> Bool
(Orbit a -> Orbit a -> Bool)
-> (Orbit a -> Orbit a -> Bool) -> Eq (Orbit a)
forall a. Eq a => Orbit a -> Orbit a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Orbit a -> Orbit a -> Bool
$c/= :: forall a. Eq a => Orbit a -> Orbit a -> Bool
== :: Orbit a -> Orbit a -> Bool
$c== :: forall a. Eq a => Orbit a -> Orbit a -> Bool
Eq)

-- | Along with 'PeriapsisSpecifier' the 'InclinationSpecifier' describes
-- orbital elements extra to its geometry.
data InclinationSpecifier a = -- | The orbit does not lie exactly in the
                              -- reference plane
                              Inclined { -- | The longitude of the ascending
                                         -- node, Ω.
                                         --
                                         -- The angle between the reference
                                         -- direction and the point where the
                                         -- orbiting body crosses the reference
                                         -- plane in the positive z direction.
                                         InclinationSpecifier a -> Angle a
longitudeOfAscendingNode :: !(Angle a)
                                         -- | The orbit's inclination, i.
                                         --
                                         -- The angle between the reference
                                         -- plane and the orbital plane
                                       , InclinationSpecifier a -> Angle a
inclination              :: !(Angle a)
                                       }
                              -- | The orbit lies in the reference plane
                            | NonInclined
  deriving (Int -> InclinationSpecifier a -> ShowS
[InclinationSpecifier a] -> ShowS
InclinationSpecifier a -> String
(Int -> InclinationSpecifier a -> ShowS)
-> (InclinationSpecifier a -> String)
-> ([InclinationSpecifier a] -> ShowS)
-> Show (InclinationSpecifier a)
forall a. Show a => Int -> InclinationSpecifier a -> ShowS
forall a. Show a => [InclinationSpecifier a] -> ShowS
forall a. Show a => InclinationSpecifier a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [InclinationSpecifier a] -> ShowS
$cshowList :: forall a. Show a => [InclinationSpecifier a] -> ShowS
show :: InclinationSpecifier a -> String
$cshow :: forall a. Show a => InclinationSpecifier a -> String
showsPrec :: Int -> InclinationSpecifier a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> InclinationSpecifier a -> ShowS
Show, InclinationSpecifier a -> InclinationSpecifier a -> Bool
(InclinationSpecifier a -> InclinationSpecifier a -> Bool)
-> (InclinationSpecifier a -> InclinationSpecifier a -> Bool)
-> Eq (InclinationSpecifier a)
forall a.
Eq a =>
InclinationSpecifier a -> InclinationSpecifier a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: InclinationSpecifier a -> InclinationSpecifier a -> Bool
$c/= :: forall a.
Eq a =>
InclinationSpecifier a -> InclinationSpecifier a -> Bool
== :: InclinationSpecifier a -> InclinationSpecifier a -> Bool
$c== :: forall a.
Eq a =>
InclinationSpecifier a -> InclinationSpecifier a -> Bool
Eq)

-- | Along with 'InclinationSpecifier' the 'PeriapsisSpecifier' describes
-- orbital elements extra to its geometry.
data PeriapsisSpecifier a = -- | The orbit is not circular
                            Eccentric { -- | The argument of periapsis, ω.
                                        --
                                        -- The 'argumentOfPeriapsis' is the
                                        -- angle of the periapsis relative to
                                        -- the reference direction in the
                                        -- orbital plane.
                                        PeriapsisSpecifier a -> Angle a
argumentOfPeriapsis :: !(Angle a)
                                      }
                            -- | The orbit has an eccentricity of 0 so the
                            -- 'argumentOfPeriapsis' is indeterminate.
                          | Circular
  deriving (Int -> PeriapsisSpecifier a -> ShowS
[PeriapsisSpecifier a] -> ShowS
PeriapsisSpecifier a -> String
(Int -> PeriapsisSpecifier a -> ShowS)
-> (PeriapsisSpecifier a -> String)
-> ([PeriapsisSpecifier a] -> ShowS)
-> Show (PeriapsisSpecifier a)
forall a. Show a => Int -> PeriapsisSpecifier a -> ShowS
forall a. Show a => [PeriapsisSpecifier a] -> ShowS
forall a. Show a => PeriapsisSpecifier a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [PeriapsisSpecifier a] -> ShowS
$cshowList :: forall a. Show a => [PeriapsisSpecifier a] -> ShowS
show :: PeriapsisSpecifier a -> String
$cshow :: forall a. Show a => PeriapsisSpecifier a -> String
showsPrec :: Int -> PeriapsisSpecifier a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> PeriapsisSpecifier a -> ShowS
Show, PeriapsisSpecifier a -> PeriapsisSpecifier a -> Bool
(PeriapsisSpecifier a -> PeriapsisSpecifier a -> Bool)
-> (PeriapsisSpecifier a -> PeriapsisSpecifier a -> Bool)
-> Eq (PeriapsisSpecifier a)
forall a.
Eq a =>
PeriapsisSpecifier a -> PeriapsisSpecifier a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: PeriapsisSpecifier a -> PeriapsisSpecifier a -> Bool
$c/= :: forall a.
Eq a =>
PeriapsisSpecifier a -> PeriapsisSpecifier a -> Bool
== :: PeriapsisSpecifier a -> PeriapsisSpecifier a -> Bool
$c== :: forall a.
Eq a =>
PeriapsisSpecifier a -> PeriapsisSpecifier a -> Bool
Eq)

-- | What form the orbit's geometry takes. This is dependant only on the
-- 'eccentricity', e >= 0, of the orbit.
data Classification = -- | 0 <= e < 1
                      --
                      -- This includes circular orbits.
                      Elliptic
                      -- | e == 1
                    | Parabolic
                      -- | e > 1
                    | Hyperbolic
  deriving (Int -> Classification -> ShowS
[Classification] -> ShowS
Classification -> String
(Int -> Classification -> ShowS)
-> (Classification -> String)
-> ([Classification] -> ShowS)
-> Show Classification
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Classification] -> ShowS
$cshowList :: [Classification] -> ShowS
show :: Classification -> String
$cshow :: Classification -> String
showsPrec :: Int -> Classification -> ShowS
$cshowsPrec :: Int -> Classification -> ShowS
Show, ReadPrec [Classification]
ReadPrec Classification
Int -> ReadS Classification
ReadS [Classification]
(Int -> ReadS Classification)
-> ReadS [Classification]
-> ReadPrec Classification
-> ReadPrec [Classification]
-> Read Classification
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [Classification]
$creadListPrec :: ReadPrec [Classification]
readPrec :: ReadPrec Classification
$creadPrec :: ReadPrec Classification
readList :: ReadS [Classification]
$creadList :: ReadS [Classification]
readsPrec :: Int -> ReadS Classification
$creadsPrec :: Int -> ReadS Classification
Read, Classification -> Classification -> Bool
(Classification -> Classification -> Bool)
-> (Classification -> Classification -> Bool) -> Eq Classification
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Classification -> Classification -> Bool
$c/= :: Classification -> Classification -> Bool
== :: Classification -> Classification -> Bool
$c== :: Classification -> Classification -> Bool
Eq)

-- TODO, use the neat "UnsafeQu" newtype for unsafe instances
unsafeMapUnit :: (a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit :: (a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit f :: a -> b
f = UnsafeQu u l b -> Qu u l b
forall (d :: [Factor *]) (l :: LCSU *) n.
UnsafeQu d l n -> Qu d l n
qu (UnsafeQu u l b -> Qu u l b)
-> (Qu u l a -> UnsafeQu u l b) -> Qu u l a -> Qu u l b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> b) -> UnsafeQu u l a -> UnsafeQu u l b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> b
f (UnsafeQu u l a -> UnsafeQu u l b)
-> (Qu u l a -> UnsafeQu u l a) -> Qu u l a -> UnsafeQu u l b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Qu u l a -> UnsafeQu u l a
forall (d :: [Factor *]) (l :: LCSU *) n.
Qu d l n -> UnsafeQu d l n
UnsafeQu

unsafeMapOrbit :: (a -> b) -> Orbit a -> Orbit b
unsafeMapOrbit :: (a -> b) -> Orbit a -> Orbit b
unsafeMapOrbit f :: a -> b
f (Orbit e :: Unitless a
e q :: Distance a
q i :: InclinationSpecifier a
i p :: PeriapsisSpecifier a
p μ :: Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ) = Unitless b
-> Distance b
-> InclinationSpecifier b
-> PeriapsisSpecifier b
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     b
-> Orbit b
forall a.
Unitless a
-> Distance a
-> InclinationSpecifier a
-> PeriapsisSpecifier a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
-> Orbit a
Orbit ((a -> b) -> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU b
forall a b (u :: [Factor *]) (l :: LCSU *).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit a -> b
f Qu '[] 'DefaultLCSU a
Unitless a
e)
                                           ((a -> b)
-> Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu '[ 'F Length One] 'DefaultLCSU b
forall a b (u :: [Factor *]) (l :: LCSU *).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit a -> b
f Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
q)
                                           ((a -> b) -> InclinationSpecifier a -> InclinationSpecifier b
forall a b.
(a -> b) -> InclinationSpecifier a -> InclinationSpecifier b
unsafeMapInclinationSpecifier a -> b
f InclinationSpecifier a
i)
                                           ((a -> b) -> PeriapsisSpecifier a -> PeriapsisSpecifier b
forall a b.
(a -> b) -> PeriapsisSpecifier a -> PeriapsisSpecifier b
unsafeMapPeriapsisSpecifier a -> b
f PeriapsisSpecifier a
p)
                                           ((a -> b)
-> Qu
     '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu
     '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU b
forall a b (u :: [Factor *]) (l :: LCSU *).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit a -> b
f Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ)

unsafeMapInclinationSpecifier :: (a -> b)
                              -> InclinationSpecifier a -> InclinationSpecifier b
unsafeMapInclinationSpecifier :: (a -> b) -> InclinationSpecifier a -> InclinationSpecifier b
unsafeMapInclinationSpecifier f :: a -> b
f s :: InclinationSpecifier a
s = case InclinationSpecifier a
s of
  Inclined _Ω :: Angle a
_Ω i :: Angle a
i -> Angle b -> Angle b -> InclinationSpecifier b
forall a. Angle a -> Angle a -> InclinationSpecifier a
Inclined ((a -> b)
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU b
forall a b (u :: [Factor *]) (l :: LCSU *).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit a -> b
f Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
_Ω) ((a -> b)
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU b
forall a b (u :: [Factor *]) (l :: LCSU *).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit a -> b
f Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
i)
  NonInclined   -> InclinationSpecifier b
forall a. InclinationSpecifier a
NonInclined

unsafeMapPeriapsisSpecifier :: (a -> b)
                            -> PeriapsisSpecifier a -> PeriapsisSpecifier b
unsafeMapPeriapsisSpecifier :: (a -> b) -> PeriapsisSpecifier a -> PeriapsisSpecifier b
unsafeMapPeriapsisSpecifier f :: a -> b
f p :: PeriapsisSpecifier a
p = case PeriapsisSpecifier a
p of
  Circular    -> PeriapsisSpecifier b
forall a. PeriapsisSpecifier a
Circular
  Eccentric a :: Angle a
a -> Angle b -> PeriapsisSpecifier b
forall a. Angle a -> PeriapsisSpecifier a
Eccentric ((a -> b)
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU b
forall a b (u :: [Factor *]) (l :: LCSU *).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit a -> b
f Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
a)

--------------------------------------------------------------------------------
-- Functions
--------------------------------------------------------------------------------

-- | Determines if the orbital elements are valid (@e >= 0@ etc...). The
-- behavior of all the other functions in this module is undefined when given
-- an invalid orbit.
isValid :: (Ord a, Num a) => Orbit a -> Bool
isValid :: Orbit a -> Bool
isValid o :: Orbit a
o = Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a -> Bool
forall a. Ord a => a -> a -> Bool
>= 0 Bool -> Bool -> Bool
&&
            ((Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a -> Bool
forall a. Eq a => a -> a -> Bool
== 0) Bool -> Bool -> Bool
`iff` (Orbit a -> PeriapsisSpecifier a
forall a. Orbit a -> PeriapsisSpecifier a
periapsisSpecifier Orbit a
o PeriapsisSpecifier a -> PeriapsisSpecifier a -> Bool
forall a. Eq a => a -> a -> Bool
== PeriapsisSpecifier a
forall a. PeriapsisSpecifier a
Circular)) Bool -> Bool -> Bool
&&
            Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
q Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu '[ 'F Length One] 'DefaultLCSU a -> Bool
forall a. Ord a => a -> a -> Bool
> Qu '[ 'F Length One] 'DefaultLCSU a
forall n (dimspec :: [Factor *]) (l :: LCSU *).
Num n =>
Qu dimspec l n
zero Bool -> Bool -> Bool
&&
            Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu
     '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Bool
forall a. Ord a => a -> a -> Bool
> Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
forall n (dimspec :: [Factor *]) (l :: LCSU *).
Num n =>
Qu dimspec l n
zero
  where
    iff :: Bool -> Bool -> Bool
iff = Bool -> Bool -> Bool
forall a. Eq a => a -> a -> Bool
(==) :: Bool -> Bool -> Bool
    e :: Unitless a
e = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
    q :: Distance a
q = Orbit a -> Distance a
forall a. Orbit a -> Distance a
periapsis Orbit a
o
    μ :: Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ = Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
forall a.
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
primaryGravitationalParameter Orbit a
o

-- | What shape is the orbit
classify :: (Num a, Ord a) => Orbit a -> Classification
classify :: Orbit a -> Classification
classify o :: Orbit a
o | Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a -> Bool
forall a. Ord a => a -> a -> Bool
< 1     = Classification
Elliptic
           | Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a -> Bool
forall a. Eq a => a -> a -> Bool
== 1    = Classification
Parabolic
           | Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a -> Bool
forall a. Ord a => a -> a -> Bool
> 1     = Classification
Hyperbolic
           | Bool
otherwise = String -> Classification
forall a. HasCallStack => String -> a
error "classify: NaN eccentricity"
  where e :: Unitless a
e = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o

-- | Return an equivalent orbit such that
--
-- - i ∈ [0..π)
-- - Ω ∈ [0..2π)
-- - ω ∈ [0..2π)
-- - inclinationSpecifier == NonInclined if i = 0
-- - periapsisSpecifier == Circular if e == 0 and ω == 0
normalizeOrbit :: (Floating a, Real a) => Orbit a -> Orbit a
normalizeOrbit :: Orbit a -> Orbit a
normalizeOrbit (Orbit e :: Unitless a
e q :: Distance a
q inc :: InclinationSpecifier a
inc per :: PeriapsisSpecifier a
per μ :: Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ) = Unitless a
-> Distance a
-> InclinationSpecifier a
-> PeriapsisSpecifier a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
-> Orbit a
forall a.
Unitless a
-> Distance a
-> InclinationSpecifier a
-> PeriapsisSpecifier a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
-> Orbit a
Orbit Unitless a
e Distance a
q InclinationSpecifier a
inc' PeriapsisSpecifier a
per' Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ
 where
  -- Were we actually given a descending node and have to flip things
  (inc' :: InclinationSpecifier a
inc', flipped :: Bool
flipped) = case InclinationSpecifier a
inc of
    NonInclined              -> (InclinationSpecifier a
forall a. InclinationSpecifier a
NonInclined, Bool
False)
    Inclined _ i :: Angle a
i | Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
i Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Bool
forall a. Eq a => a -> a -> Bool
== Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall n (dimspec :: [Factor *]) (l :: LCSU *).
Num n =>
Qu dimspec l n
zero -> (InclinationSpecifier a
forall a. InclinationSpecifier a
NonInclined, Bool
False)
    Inclined _Ω :: Angle a
_Ω i :: Angle a
i ->
      let iR :: Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
iR  = Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
i Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a (u :: [Factor *]) (l :: LCSU *).
Real a =>
Qu u l a -> Qu u l a -> Qu u l a
`mod'` Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
turn
          i' :: Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
i'  = if Bool
flipped then Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
turn Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|-| Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
iR else Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
iR
          _Ω' :: Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
_Ω' = (if Bool
flipped then Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
_Ω Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
halfTurn else Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
_Ω) Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a (u :: [Factor *]) (l :: LCSU *).
Real a =>
Qu u l a -> Qu u l a -> Qu u l a
`mod'` Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
turn
      in  (Angle a -> Angle a -> InclinationSpecifier a
forall a. Angle a -> Angle a -> InclinationSpecifier a
Inclined Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
_Ω' Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
i', Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
iR Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Bool
forall a. Ord a => a -> a -> Bool
>= Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
halfTurn)

  per' :: PeriapsisSpecifier a
per' = case PeriapsisSpecifier a
per of
    Circular | Bool
flipped   -> Angle a -> PeriapsisSpecifier a
forall a. Angle a -> PeriapsisSpecifier a
Eccentric Angle a
forall a. Floating a => PlaneAngle a
halfTurn
             | Bool
otherwise -> PeriapsisSpecifier a
forall a. PeriapsisSpecifier a
Circular
    Eccentric ω :: Angle a
ω
      | Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
ω Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Bool
forall a. Eq a => a -> a -> Bool
== Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall n (dimspec :: [Factor *]) (l :: LCSU *).
Num n =>
Qu dimspec l n
zero, Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a -> Bool
forall a. Eq a => a -> a -> Bool
== 0, Bool -> Bool
not Bool
flipped
      -> PeriapsisSpecifier a
forall a. PeriapsisSpecifier a
Circular
      | Bool
otherwise
      -> Angle a -> PeriapsisSpecifier a
forall a. Angle a -> PeriapsisSpecifier a
Eccentric (Angle a -> PeriapsisSpecifier a)
-> Angle a -> PeriapsisSpecifier a
forall a b. (a -> b) -> a -> b
$ (if Bool
flipped then Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
ω Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
halfTurn else Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
ω) Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a (u :: [Factor *]) (l :: LCSU *).
Real a =>
Qu u l a -> Qu u l a -> Qu u l a
`mod'` Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
turn

-- | Calculate the semi-major axis, a, of the 'Orbit'. Returns 'Nothing' when
-- given a parabolic orbit for which there is no semi-major axis. Note that the
-- semi-major axis of a hyperbolic orbit is negative.
semiMajorAxis :: (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)
semiMajorAxis :: Orbit a -> Maybe (Distance a)
semiMajorAxis o :: Orbit a
o =
  case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
    Parabolic -> Maybe (Distance a)
forall a. Maybe a
Nothing
    _         -> Qu ('F Length One : Normalize '[]) 'DefaultLCSU a
-> Maybe (Distance a)
forall a. a -> Maybe a
Just (Qu ('F Length One : Normalize '[]) 'DefaultLCSU a
 -> Maybe (Distance a))
-> Qu ('F Length One : Normalize '[]) 'DefaultLCSU a
-> Maybe (Distance a)
forall a b. (a -> b) -> a -> b
$ Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
q Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F Length One] @- '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| (1 Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|-| Qu '[] 'DefaultLCSU a
Unitless a
e)
  where
    q :: Distance a
q = Orbit a -> Distance a
forall a. Orbit a -> Distance a
periapsis Orbit a
o
    e :: Unitless a
e = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o

-- | Calculate the semi-minor axis, b, of the 'Orbit'. Like 'semiMajorAxis'
-- @\'semiMinorAxis\' o@ is negative when @o@ is a hyperbolic orbit. In the
-- case of a parabolic orbit 'semiMinorAxis' returns @0m@.
semiMinorAxis :: (Floating a, Ord a) => Orbit a -> Distance a
semiMinorAxis :: Orbit a -> Distance a
semiMinorAxis o :: Orbit a
o =
  case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
    Elliptic   -> Qu '[ 'F Length One] 'DefaultLCSU a
a Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F Length One] @+ '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[] 'DefaultLCSU a -> Qu ('[] @/ 'S One) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (1 Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|-| Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a -> Int -> Qu '[] 'DefaultLCSU a
forall a b. (Num a, Integral b) => a -> b -> a
^ (2::Int))
    Parabolic  -> Distance a
forall n (dimspec :: [Factor *]) (l :: LCSU *).
Num n =>
Qu dimspec l n
zero
    Hyperbolic -> Qu '[ 'F Length One] 'DefaultLCSU a
a Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F Length One] @+ '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[] 'DefaultLCSU a -> Qu ('[] @/ 'S One) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a -> Int -> Qu '[] 'DefaultLCSU a
forall a b. (Num a, Integral b) => a -> b -> a
^ (2::Int) Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|-| 1)
  where
    e :: Unitless a
e = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
    Just a :: Qu '[ 'F Length One] 'DefaultLCSU a
a = Orbit a -> Maybe (Distance a)
forall a. (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)
semiMajorAxis Orbit a
o

-- | Calculate the semiLatusRectum, l, of the 'Orbit'
semiLatusRectum :: (Num a) => Orbit a -> Distance a
semiLatusRectum :: Orbit a -> Distance a
semiLatusRectum orbit :: Orbit a
orbit = Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a
-> Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu (Normalize ('[] @+ '[ 'F Length One])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
q Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu '[ 'F Length One] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
q
  where q :: Distance a
q = Orbit a -> Distance a
forall a. Orbit a -> Distance a
periapsis Orbit a
orbit
        e :: Unitless a
e = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
orbit

-- | Calculate the distance between the bodies when they are at their most
-- distant. 'apoapsis' returns 'Nothing' when given a parabolic or hyperbolic
-- orbit.
apoapsis :: (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)
apoapsis :: Orbit a -> Maybe (Distance a)
apoapsis o :: Orbit a
o =
  case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
    Elliptic -> Qu ('F Length One : Normalize '[]) 'DefaultLCSU a
-> Maybe (Distance a)
forall a. a -> Maybe a
Just (Qu ('F Length One : Normalize '[]) 'DefaultLCSU a
 -> Maybe (Distance a))
-> Qu ('F Length One : Normalize '[]) 'DefaultLCSU a
-> Maybe (Distance a)
forall a b. (a -> b) -> a -> b
$ Qu '[ 'F Length One] 'DefaultLCSU a
a Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F Length One] @+ '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| (1 Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| Qu '[] 'DefaultLCSU a
Unitless a
e)
    _        -> Maybe (Distance a)
forall a. Maybe a
Nothing
  where
    Just a :: Qu '[ 'F Length One] 'DefaultLCSU a
a = Orbit a -> Maybe (Distance a)
forall a. (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)
semiMajorAxis Orbit a
o
    e :: Unitless a
e = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o

-- | Calculate the mean motion, n, of an orbit
--
-- This is the rate of change of the mean anomaly with respect to time.
meanMotion :: (Floating a, Ord a) => Orbit a -> Quantity [si|rad/s|] a
meanMotion :: Orbit a -> Quantity (Radian :/ Second) a
meanMotion o :: Orbit a
o =
  case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
    Elliptic   -> Qu '[ 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Quantity (Radian :/ Second) a
forall (b :: [Factor *]) a.
Qu b 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ b)) 'DefaultLCSU a
addRad (Qu '[ 'F Time ('P 'Zero)] 'DefaultLCSU a
 -> Quantity (Radian :/ Second) a)
-> Qu '[ 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Quantity (Radian :/ Second) a
forall a b. (a -> b) -> a -> b
$ Qu '[ 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu ('[ 'F Time ('P ('P 'Zero))] @/ 'S One) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu '[ 'F Length ('S ('S One))] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
         @- '[ 'F Length ('S ('S One))]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu
     (Normalize
        (Normalize ('[ 'F Length One] @+ '[ 'F Length One])
         @+ '[ 'F Length One]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *).
Num n =>
Qu a l n -> Qu (Normalize (Normalize (a @+ a) @+ a)) l n
qCube Qu '[ 'F Length One] 'DefaultLCSU a
a)
    Hyperbolic -> Qu '[ 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Quantity (Radian :/ Second) a
forall (b :: [Factor *]) a.
Qu b 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ b)) 'DefaultLCSU a
addRad (Qu '[ 'F Time ('P 'Zero)] 'DefaultLCSU a
 -> Quantity (Radian :/ Second) a)
-> Qu '[ 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Quantity (Radian :/ Second) a
forall a b. (a -> b) -> a -> b
$ Qu '[ 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu ('[ 'F Time ('P ('P 'Zero))] @/ 'S One) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu '[ 'F Length ('S ('S One))] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
         @- '[ 'F Length ('S ('S One))]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Qu '[ 'F Length ('S ('S One))] 'DefaultLCSU a
-> Qu '[ 'F Length ('S ('S One))] 'DefaultLCSU a
forall n (d :: [Factor *]) (l :: LCSU *).
Num n =>
Qu d l n -> Qu d l n
qNegate (Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu
     (Normalize
        (Normalize ('[ 'F Length One] @+ '[ 'F Length One])
         @+ '[ 'F Length One]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *).
Num n =>
Qu a l n -> Qu (Normalize (Normalize (a @+ a) @+ a)) l n
qCube Qu '[ 'F Length One] 'DefaultLCSU a
a))
    Parabolic  -> Qu '[ 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Quantity (Radian :/ Second) a
forall (b :: [Factor *]) a.
Qu b 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ b)) 'DefaultLCSU a
addRad (Qu '[ 'F Time ('P 'Zero)] 'DefaultLCSU a
 -> Quantity (Radian :/ Second) a)
-> Qu '[ 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Quantity (Radian :/ Second) a
forall a b. (a -> b) -> a -> b
$ 2 Qu '[] 'DefaultLCSU a
-> Qu '[ 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Qu (Normalize ('[] @+ '[ 'F Time ('P 'Zero)])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[ 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu ('[ 'F Time ('P ('P 'Zero))] @/ 'S One) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu '[ 'F Length ('S ('S One))] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
         @- '[ 'F Length ('S ('S One))]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu
     (Normalize
        (Normalize ('[ 'F Length One] @+ '[ 'F Length One])
         @+ '[ 'F Length One]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *).
Num n =>
Qu a l n -> Qu (Normalize (Normalize (a @+ a) @+ a)) l n
qCube Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
l)
  where
    Just a :: Qu '[ 'F Length One] 'DefaultLCSU a
a = Orbit a -> Maybe (Distance a)
forall a. (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)
semiMajorAxis Orbit a
o
    μ :: Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ = Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
forall a.
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
primaryGravitationalParameter Orbit a
o
    l :: Distance a
l = Orbit a -> Distance a
forall a. Num a => Orbit a -> Distance a
semiLatusRectum Orbit a
o

-- | Calculate the orbital period, p, of an elliptic orbit.
--
-- 'period' returns Nothing if given a parabolic or hyperbolic orbit.
period :: (Floating a, Ord a) => Orbit a -> Maybe (Time a)
period :: Orbit a -> Maybe (Time a)
period o :: Orbit a
o =
  case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
    Elliptic -> Qu ('F Time One : Normalize '[]) 'DefaultLCSU a
-> Maybe (Qu ('F Time One : Normalize '[]) 'DefaultLCSU a)
forall a. a -> Maybe a
Just Qu ('F Time One : Normalize '[]) 'DefaultLCSU a
Qu
  (Normalize
     ('[ 'F PlaneAngle One]
      @- '[ 'F PlaneAngle One, 'F Time ('P 'Zero)]))
  'DefaultLCSU
  a
p
    _ -> Maybe (Time a)
forall a. Maybe a
Nothing
  where
    n :: Quantity (Radian :/ Second) a
n = Orbit a -> Quantity (Radian :/ Second) a
forall a.
(Floating a, Ord a) =>
Orbit a -> Quantity (Radian :/ Second) a
meanMotion Orbit a
o
    p :: Qu
  (Normalize
     ('[ 'F PlaneAngle One]
      @- '[ 'F PlaneAngle One, 'F Time ('P 'Zero)]))
  'DefaultLCSU
  a
p = Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
turn Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One, 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F PlaneAngle One]
         @- '[ 'F PlaneAngle One, 'F Time ('P 'Zero)]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Qu '[ 'F PlaneAngle One, 'F Time ('P 'Zero)] 'DefaultLCSU a
Quantity (Radian :/ Second) a
n


-- | Calculate the areal velocity, A, of the orbit.
--
-- The areal velocity is the area <https://xkcd.com/21/ swept out> by the line
-- between the orbiting body and the primary per second.
arealVelocity :: (Ord a, Floating a) => Orbit a -> Quantity [si|m^2/s|] a
arealVelocity :: Orbit a -> Quantity ((Meter :^ Succ (Succ 'Zero)) :/ Second) a
arealVelocity o :: Orbit a
o = Qu
  '[ 'F Length ('S ('S ('S One))), 'F Time ('P ('P 'Zero))]
  'DefaultLCSU
  a
-> Qu
     ('[ 'F Length ('S ('S ('S One))), 'F Time ('P ('P 'Zero))]
      @/ 'S One)
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
l Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu
     '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length One]
         @+ '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ) Qu '[ 'F Length ('S One), 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu
     (Normalize ('[ 'F Length ('S One), 'F Time ('P 'Zero)] @- '[]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| 2
  where l :: Distance a
l = Orbit a -> Distance a
forall a. Num a => Orbit a -> Distance a
semiLatusRectum Orbit a
o
        μ :: Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ = Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
forall a.
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
primaryGravitationalParameter Orbit a
o

-- | Calculate the angle at which a body leaves the system when on an escape
-- trajectory relative to the argument of periapsis. This is the limit of the
-- true anomaly as time tends towards infinity minus the argument of periapsis.
-- The departure angle is in the closed range (π/2..π).
--
-- This is the negation of the approach angle.
--
-- 'hyperbolicDepartureAngle' returns Nothing when given an elliptic orbit and
-- π when given a parabolic orbit.
hyperbolicDepartureAngle :: (Floating a, Ord a) => Orbit a -> Maybe (Angle a)
hyperbolicDepartureAngle :: Orbit a -> Maybe (Angle a)
hyperbolicDepartureAngle o :: Orbit a
o =
  case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
    Hyperbolic ->
      let e :: Unitless a
e = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
          θ :: Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
θ = Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
forall (b :: [Factor *]) a.
Qu b 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ b)) 'DefaultLCSU a
addRad (Qu '[] 'DefaultLCSU a
 -> Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a)
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
forall a b. (a -> b) -> a -> b
$ Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Floating a => a -> a
acos (-1 Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Fractional a => a -> a -> a
/ Qu '[] 'DefaultLCSU a
Unitless a
e)
      in Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall a. a -> Maybe a
Just Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
θ
    Parabolic -> Qu ('F PlaneAngle One : Normalize '[]) 'DefaultLCSU a
-> Maybe (Qu ('F PlaneAngle One : Normalize '[]) 'DefaultLCSU a)
forall a. a -> Maybe a
Just (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
turn Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @- '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| 2)
    _ -> Maybe (Angle a)
forall a. Maybe a
Nothing

-- | Calculate the angle at which a body leaves the system when on a hyperbolic
-- trajectory relative to the argument of periapsis. This is the limit of the
-- true anomaly as time tends towards -infinity minus the argument of
-- periapsis. The approach angle is in the closed range (-π..π/2).
--
-- This is the negation of the departure angle.
--
-- 'hyperbolicApproachAngle' returns Nothing when given a non-hyperbolic orbit
-- and -π when given a parabolic orbit.
hyperbolicApproachAngle :: (Floating a, Ord a) => Orbit a -> Maybe (Angle a)
hyperbolicApproachAngle :: Orbit a -> Maybe (Angle a)
hyperbolicApproachAngle = (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
 -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall n (d :: [Factor *]) (l :: LCSU *).
Num n =>
Qu d l n -> Qu d l n
qNegate (Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
 -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> (Orbit a -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> Orbit a
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Orbit a -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall a. (Floating a, Ord a) => Orbit a -> Maybe (Angle a)
hyperbolicDepartureAngle

----------------------------------------------------------------
-- ## Conversions between time and anomolies
----------------------------------------------------------------

---------
-- To time
---------

-- | Calculate the time since periapse, t, when the body has the given
-- <https://en.wikipedia.org/wiki/Mean_anomaly mean anomaly>, M. M may be
-- negative, indicating that the orbiting body has yet to reach periapse.
--
-- The sign of the time at mean anomaly M is the same as the sign of M.
--
-- The returned time is unbounded.
timeAtMeanAnomaly :: (Floating a, Ord a) => Orbit a -> Angle a -> Time a
timeAtMeanAnomaly :: Orbit a -> Angle a -> Time a
timeAtMeanAnomaly o :: Orbit a
o _M :: Angle a
_M = Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
_M Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One, 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F PlaneAngle One]
         @- '[ 'F PlaneAngle One, 'F Time ('P 'Zero)]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Qu '[ 'F PlaneAngle One, 'F Time ('P 'Zero)] 'DefaultLCSU a
Quantity (Radian :/ Second) a
n
  where n :: Quantity (Radian :/ Second) a
n = Orbit a -> Quantity (Radian :/ Second) a
forall a.
(Floating a, Ord a) =>
Orbit a -> Quantity (Radian :/ Second) a
meanMotion Orbit a
o

-- | Calculate the time since periapse, t, of an elliptic orbit when at
-- eccentric anomaly E.
--
-- 'timeAtEccentricAnomaly' returns Nothing if given a parabolic or hyperbolic
-- orbit.
timeAtEccentricAnomaly :: (Floating a, Ord a) => Orbit a -> Angle a -> Maybe (Time a)
timeAtEccentricAnomaly :: Orbit a -> Angle a -> Maybe (Time a)
timeAtEccentricAnomaly o :: Orbit a
o = (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
 -> Qu '[ 'F Time One] 'DefaultLCSU a)
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> Maybe (Qu '[ 'F Time One] 'DefaultLCSU a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Orbit a -> Angle a -> Time a
forall a. (Floating a, Ord a) => Orbit a -> Angle a -> Time a
timeAtMeanAnomaly Orbit a
o) (Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
 -> Maybe (Qu '[ 'F Time One] 'DefaultLCSU a))
-> (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
    -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F Time One] 'DefaultLCSU a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Orbit a -> Angle a -> Maybe (Angle a)
forall a.
(Floating a, Ord a) =>
Orbit a -> Angle a -> Maybe (Angle a)
meanAnomalyAtEccentricAnomaly Orbit a
o

-- | Calculate the time since periapse, t, of a hyperbolic orbit when at
-- hyperbolic anomaly H.
--
-- Returns Nothing if given an elliptic or parabolic orbit.
timeAtHyperbolicAnomaly
  :: (Floating a, Ord a) => Orbit a -> AngleH a -> Maybe (Time a)
timeAtHyperbolicAnomaly :: Orbit a -> AngleH a -> Maybe (Time a)
timeAtHyperbolicAnomaly o :: Orbit a
o =
  (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
 -> Qu '[ 'F Time One] 'DefaultLCSU a)
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> Maybe (Qu '[ 'F Time One] 'DefaultLCSU a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Orbit a -> Angle a -> Time a
forall a. (Floating a, Ord a) => Orbit a -> Angle a -> Time a
timeAtMeanAnomaly Orbit a
o) (Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
 -> Maybe (Qu '[ 'F Time One] 'DefaultLCSU a))
-> (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
    -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F Time One] 'DefaultLCSU a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Orbit a -> AngleH a -> Maybe (Angle a)
forall a.
(Floating a, Ord a) =>
Orbit a -> AngleH a -> Maybe (Angle a)
meanAnomalyAtHyperbolicAnomaly Orbit a
o

-- | Calculate the time since periapse given the true anomaly, ν, of an
-- orbiting body.
--
-- Returns 'Nothing' if the body never passed through the specified true
-- anomaly.
timeAtTrueAnomaly
  :: (Real a, Floating a) => Orbit a -> Angle a -> Maybe (Time a)
timeAtTrueAnomaly :: Orbit a -> Angle a -> Maybe (Time a)
timeAtTrueAnomaly o :: Orbit a
o ν :: Angle a
ν = case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  _ | Just d :: Angle a
d <- Orbit a -> Maybe (Angle a)
forall a. (Floating a, Ord a) => Orbit a -> Maybe (Angle a)
hyperbolicDepartureAngle Orbit a
o, Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a (l :: LCSU *) (u :: [Factor *]).
Num a =>
Qu u l a -> Qu u l a
qAbs Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
ν Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Bool
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Ord n) =>
Qu d1 l n -> Qu d2 l n -> Bool
|>| Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
d -> Maybe (Time a)
forall a. Maybe a
Nothing
  Parabolic ->
    let _D :: Unitless a
_D = Angle a -> Unitless a
forall a. Floating a => Angle a -> Unitless a
qTan (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
ν Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @- '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| 2)
        t :: Qu (Normalize ('[ 'F Time One] @+ '[])) 'DefaultLCSU a
t  = 0.5 Qu '[] 'DefaultLCSU a
-> Qu '[ 'F Time One] 'DefaultLCSU a
-> Qu (Normalize ('[] @+ '[ 'F Time One])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[ 'F Time ('S One)] 'DefaultLCSU a
-> Qu ('[ 'F Time ('S One)] @/ 'S One) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu
     (Normalize
        (Normalize ('[ 'F Length One] @+ '[ 'F Length One])
         @+ '[ 'F Length One]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *).
Num n =>
Qu a l n -> Qu (Normalize (Normalize (a @+ a) @+ a)) l n
qCube Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
l Qu '[ 'F Length ('S ('S One))] 'DefaultLCSU a
-> Qu
     '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S ('S One))]
         @- '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ) Qu '[ 'F Time One] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F Time One] @+ '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| (Qu '[] 'DefaultLCSU a
Unitless a
_D Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| (Qu '[] 'DefaultLCSU a
-> Qu (Normalize (Normalize ('[] @+ '[]) @+ '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *).
Num n =>
Qu a l n -> Qu (Normalize (Normalize (a @+ a) @+ a)) l n
qCube Qu '[] 'DefaultLCSU a
Unitless a
_D Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[] @- '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| 3))
    in  Qu ('F Time One : Normalize '[]) 'DefaultLCSU a
-> Maybe (Qu ('F Time One : Normalize '[]) 'DefaultLCSU a)
forall a. a -> Maybe a
Just Qu ('F Time One : Normalize '[]) 'DefaultLCSU a
Qu (Normalize ('[ 'F Time One] @+ '[])) 'DefaultLCSU a
t
  _ -> (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
 -> Qu '[ 'F Time One] 'DefaultLCSU a)
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> Maybe (Qu '[ 'F Time One] 'DefaultLCSU a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Orbit a -> Angle a -> Time a
forall a. (Floating a, Ord a) => Orbit a -> Angle a -> Time a
timeAtMeanAnomaly Orbit a
o) (Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
 -> Maybe (Qu '[ 'F Time One] 'DefaultLCSU a))
-> (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
    -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F Time One] 'DefaultLCSU a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Orbit a -> Angle a -> Maybe (Angle a)
forall a.
(Real a, Floating a) =>
Orbit a -> Angle a -> Maybe (Angle a)
meanAnomalyAtTrueAnomaly Orbit a
o (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Maybe (Time a))
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Maybe (Time a)
forall a b. (a -> b) -> a -> b
$ Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
ν
 where
  μ :: Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ = Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
forall a.
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
primaryGravitationalParameter Orbit a
o
  l :: Distance a
l = Orbit a -> Distance a
forall a. Num a => Orbit a -> Distance a
semiLatusRectum Orbit a
o

---------
-- To mean anomaly
---------

-- | Calculate the <https://en.wikipedia.org/wiki/Mean_anomaly mean anomaly>,
-- M, at the given time since periapse, t. t may be negative, indicating that
-- the orbiting body has yet to reach periapse.
--
-- The sign of the mean anomaly at time t is the same as the sign of t.
--
-- The returned mean anomaly is unbounded.
meanAnomalyAtTime :: (Floating a, Ord a) => Orbit a -> Time a -> Angle a
meanAnomalyAtTime :: Orbit a -> Time a -> Angle a
meanAnomalyAtTime o :: Orbit a
o t :: Time a
t = Qu '[ 'F Time One] 'DefaultLCSU a
Time a
t Qu '[ 'F Time One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One, 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Time One] @+ '[ 'F PlaneAngle One, 'F Time ('P 'Zero)]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[ 'F PlaneAngle One, 'F Time ('P 'Zero)] 'DefaultLCSU a
Quantity (Radian :/ Second) a
n
  where n :: Quantity (Radian :/ Second) a
n = Orbit a -> Quantity (Radian :/ Second) a
forall a.
(Floating a, Ord a) =>
Orbit a -> Quantity (Radian :/ Second) a
meanMotion Orbit a
o

-- | Calculate the mean anomaly, M, of an elliptic orbit when at eccentric
-- anomaly E
--
-- 'meanAnomalyAtEccentricAnomaly' returns Nothing if given a parabolic or
-- hyperbolic orbit.
--
-- The number of orbits represented by the anomalies is preserved;
-- i.e. M `div` 2π = E `div` 2π
meanAnomalyAtEccentricAnomaly :: (Floating a, Ord a) => Orbit a -> Angle a -> Maybe (Angle a)
meanAnomalyAtEccentricAnomaly :: Orbit a -> Angle a -> Maybe (Angle a)
meanAnomalyAtEccentricAnomaly o :: Orbit a
o _E :: Angle a
_E = case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
                                       Elliptic -> Qu ('F PlaneAngle One : Normalize '[]) 'DefaultLCSU a
-> Maybe (Qu ('F PlaneAngle One : Normalize '[]) 'DefaultLCSU a)
forall a. a -> Maybe a
Just Qu ('F PlaneAngle One : Normalize '[]) 'DefaultLCSU a
Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
_M
                                       _ -> Maybe (Angle a)
forall a. Maybe a
Nothing
  where e :: Unitless a
e = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
        untypedE :: Qu
  (Normalize ('[ 'F PlaneAngle One] @- '[ 'F PlaneAngle One]))
  'DefaultLCSU
  a
untypedE = Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu
     (Normalize ('[ 'F PlaneAngle One] @- '[ 'F PlaneAngle One]))
     'DefaultLCSU
     a
forall (u :: [Factor *]) a.
Qu u 'DefaultLCSU a
-> Qu (Normalize (u @- '[ 'F PlaneAngle One])) 'DefaultLCSU a
delRad Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
_E
        _M :: Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
_M = Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
forall (b :: [Factor *]) a.
Qu b 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ b)) 'DefaultLCSU a
addRad (Qu '[] 'DefaultLCSU a
Qu
  (Normalize ('[ 'F PlaneAngle One] @- '[ 'F PlaneAngle One]))
  'DefaultLCSU
  a
untypedE Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|-| Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[] @+ '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Floating a => a -> a
sin Qu '[] 'DefaultLCSU a
Qu
  (Normalize ('[ 'F PlaneAngle One] @- '[ 'F PlaneAngle One]))
  'DefaultLCSU
  a
untypedE)

-- | Calculate the mean anomaly, M, of a hyperbolic orbit when at hyperbolic
-- anomaly H
meanAnomalyAtHyperbolicAnomaly
  :: (Floating a, Ord a) => Orbit a -> AngleH a -> Maybe (Angle a)
meanAnomalyAtHyperbolicAnomaly :: Orbit a -> AngleH a -> Maybe (Angle a)
meanAnomalyAtHyperbolicAnomaly o :: Orbit a
o _H :: AngleH a
_H = case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Hyperbolic -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall a. a -> Maybe a
Just Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
_M
  _          -> Maybe (Angle a)
forall a. Maybe a
Nothing
 where
  e :: Unitless a
e  = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
  _M :: Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
_M = Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
forall (b :: [Factor *]) a.
Qu b 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ b)) 'DefaultLCSU a
addRad (Qu '[] 'DefaultLCSU a
 -> Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a)
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
forall a b. (a -> b) -> a -> b
$ Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Num a => a -> a -> a
* AngleH a -> Unitless a
forall a. Floating a => AngleH a -> Unitless a
qSinh AngleH a
_H Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Num a => a -> a -> a
- a -> Qu '[] 'DefaultLCSU a
forall n (l :: LCSU *). n -> Qu '[] l n
quantity (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
AngleH a
_H Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
-> RadianHyperbolic -> a
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# RadianHyperbolic
RadianHyperbolic)

-- | Calculate the mean anomaly, M, of an orbiting body when at the given true
-- anomaly, ν.
--
-- The number of orbits represented by the anomalies is preserved;
-- i.e. M `div` 2π = ν `div` 2π
--
-- Returns 'Nothing' for parabolic orbits.
--
-- Returns 'Nothing' when the trajectory is not defined for the given true
-- anomaly.
meanAnomalyAtTrueAnomaly :: (Real a, Floating a)
                         => Orbit a -> Angle a -> Maybe (Angle a)
meanAnomalyAtTrueAnomaly :: Orbit a -> Angle a -> Maybe (Angle a)
meanAnomalyAtTrueAnomaly o :: Orbit a
o ν :: Angle a
ν = case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Parabolic -> Maybe (Angle a)
forall a. Maybe a
Nothing
  Elliptic -> Orbit a -> Angle a -> Maybe (Angle a)
forall a.
(Floating a, Ord a) =>
Orbit a -> Angle a -> Maybe (Angle a)
meanAnomalyAtEccentricAnomaly Orbit a
o (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
 -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
    -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall (m :: * -> *) b c a.
Monad m =>
(b -> m c) -> (a -> m b) -> a -> m c
<=<
              Orbit a -> Angle a -> Maybe (Angle a)
forall a.
(Floating a, Real a) =>
Orbit a -> Angle a -> Maybe (Angle a)
eccentricAnomalyAtTrueAnomaly Orbit a
o (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Maybe (Angle a))
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Maybe (Angle a)
forall a b. (a -> b) -> a -> b
$ Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
ν
  Hyperbolic -> Orbit a -> AngleH a -> Maybe (Angle a)
forall a.
(Floating a, Ord a) =>
Orbit a -> AngleH a -> Maybe (Angle a)
meanAnomalyAtHyperbolicAnomaly Orbit a
o (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
 -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
    -> Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a))
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall (m :: * -> *) b c a.
Monad m =>
(b -> m c) -> (a -> m b) -> a -> m c
<=<
                Orbit a -> Angle a -> Maybe (AngleH a)
forall a.
(Floating a, Ord a) =>
Orbit a -> Angle a -> Maybe (AngleH a)
hyperbolicAnomalyAtTrueAnomaly Orbit a
o (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Maybe (Angle a))
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Maybe (Angle a)
forall a b. (a -> b) -> a -> b
$ Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
ν

---------
-- To eccentric
---------

-- | Calculate the eccentric anomaly, E, of an elliptic orbit at time t.
--
-- 'eccentricAnomalyAtTime' returns Nothing when given a parabolic or
-- hyperbolic orbit.
--
-- The number of orbits represented by the time is preserved;
-- i.e. t `div` p = E `div` 2π
eccentricAnomalyAtTime :: (Converge [a], Floating a, Real a)
                       => Orbit a -> Time a -> Maybe (Angle a)
eccentricAnomalyAtTime :: Orbit a -> Time a -> Maybe (Angle a)
eccentricAnomalyAtTime o :: Orbit a
o t :: Time a
t = case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Elliptic -> Orbit a -> Angle a -> Maybe (Angle a)
forall a.
(Converge [a], Floating a, Real a) =>
Orbit a -> Angle a -> Maybe (Angle a)
eccentricAnomalyAtMeanAnomaly Orbit a
o (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
 -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> (Qu '[ 'F Time One] 'DefaultLCSU a
    -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> Qu '[ 'F Time One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Orbit a -> Time a -> Angle a
forall a. (Floating a, Ord a) => Orbit a -> Time a -> Angle a
meanAnomalyAtTime Orbit a
o (Qu '[ 'F Time One] 'DefaultLCSU a -> Maybe (Angle a))
-> Qu '[ 'F Time One] 'DefaultLCSU a -> Maybe (Angle a)
forall a b. (a -> b) -> a -> b
$ Qu '[ 'F Time One] 'DefaultLCSU a
Time a
t
  _ -> Maybe (Angle a)
forall a. Maybe a
Nothing

-- | Calculate the eccentric anomaly, E, of an elliptic orbit when at mean
-- anomaly M. This function is considerably slower than most other conversion
-- functions as it uses an iterative method as no closed form solution exists.
--
-- The number of orbits represented by the anomalies is preserved;
-- i.e. M `div` 2π = E `div` 2π
--
-- 'eccentricAnomalyAtMeanAnomaly' returns Nothing when given a parabolic or
-- hyperbolic orbit.
eccentricAnomalyAtMeanAnomaly :: forall a. (Converge [a], Floating a, Real a)
                              => Orbit a -> Angle a -> Maybe (Angle a)
eccentricAnomalyAtMeanAnomaly :: Orbit a -> Angle a -> Maybe (Angle a)
eccentricAnomalyAtMeanAnomaly o :: Orbit a
o _M :: Angle a
_M = case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
                                       Elliptic -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
Maybe (Angle a)
_E
                                       _ -> Maybe (Angle a)
forall a. Maybe a
Nothing
  where (n :: Qu '[] 'DefaultLCSU Integer
n, wrappedM :: a
wrappedM) = (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> a)
-> (Qu '[] 'DefaultLCSU Integer,
    Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> (Qu '[] 'DefaultLCSU Integer, a)
forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Radian -> a
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si|rad|]) (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
_M Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> (Qu '[] 'DefaultLCSU Integer,
    Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall a b (u :: [Factor *]) (l :: LCSU *).
(Real a, Integral b) =>
Qu u l a -> Qu u l a -> (Qu '[] l b, Qu u l a)
`divMod'` Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
turn)
        e :: a
e = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o Qu '[] 'DefaultLCSU a -> Number -> a
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si||]
        _MFloat :: Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
_MFloat = Float -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
forall a. Fractional a => a -> Angle a
rad (Float -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float)
-> (a -> Float) -> a -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac (a -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float)
-> a -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
forall a b. (a -> b) -> a -> b
$ a
wrappedM
        oFloat :: Orbit Float
oFloat = (a -> Float) -> Orbit a -> Orbit Float
forall a b. (a -> b) -> Orbit a -> Orbit b
unsafeMapOrbit a -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac Orbit a
o
        initialGuessFloat :: Angle Float
        Just initialGuessFloat = Orbit Float -> Angle Float -> Maybe (Angle Float)
eccentricAnomalyAtMeanAnomalyFloat Orbit Float
oFloat Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
Angle Float
_MFloat
        initialGuess :: a
initialGuess = Float -> a
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Float -> a)
-> (Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float -> Float)
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
-> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float -> Radian -> Float
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si|rad|]) (Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float -> a)
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float -> a
forall a b. (a -> b) -> a -> b
$ Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
Angle Float
initialGuessFloat
        err :: (Mode b, Floating b, Scalar b ~ a) => b -> b
        err :: b -> b
err _E :: b
_E = Scalar b -> b
forall t. Mode t => Scalar t -> t
auto a
Scalar b
wrappedM b -> b -> b
forall a. Num a => a -> a -> a
- (b
_E b -> b -> b
forall a. Num a => a -> a -> a
- Scalar b -> b
forall t. Mode t => Scalar t -> t
auto a
Scalar b
e b -> b -> b
forall a. Num a => a -> a -> a
* b -> b
forall a. Floating a => a -> a
sin b
_E)
        wrappedE :: Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
wrappedE = (a -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> Maybe a -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Fractional a => a -> Angle a
rad (Maybe a -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> ([a] -> Maybe a)
-> [a]
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Element [a] -> Element [a]) -> [a] -> Maybe (Element [a])
forall a.
(Converge a, Ord (Element a)) =>
(Element a -> Element a) -> a -> Maybe (Element a)
convergeErr (Id a -> a
forall a. Id a -> a
runId (Id a -> a) -> (a -> Id a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Id a -> Id a
forall a. Num a => a -> a
abs (Id a -> Id a) -> (a -> Id a) -> a -> Id a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Id a -> Id a
forall b. (Mode b, Floating b, Scalar b ~ a) => b -> b
err (Id a -> Id a) -> (a -> Id a) -> a -> Id a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.  a -> Id a
forall a. a -> Id a
Id) ([a] -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> [a] -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall a b. (a -> b) -> a -> b
$
                   (forall s. AD s (Tower a) -> AD s (Tower a)) -> a -> [a]
forall a.
Fractional a =>
(forall s. AD s (Tower a) -> AD s (Tower a)) -> a -> [a]
findZeroNoEq forall s. AD s (Tower a) -> AD s (Tower a)
forall b. (Mode b, Floating b, Scalar b ~ a) => b -> b
err a
initialGuess
        _E :: Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
_E = (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| ((Integer -> a)
-> Qu '[] 'DefaultLCSU Integer -> Qu '[] 'DefaultLCSU a
forall a b (u :: [Factor *]) (l :: LCSU *).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit Integer -> a
forall a. Num a => Integer -> a
fromInteger Qu '[] 'DefaultLCSU Integer
n Qu '[] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu (Normalize ('[] @+ '[ 'F PlaneAngle One])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
turn)) (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
 -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
wrappedE

-- | 'eccentricAnomalyAtMeanAnomaly' specialized to 'Float'.
--
-- This function is used to calculate the initial guess for
-- 'eccentricAnomalyAtMeanAnomaly'.
eccentricAnomalyAtMeanAnomalyFloat :: Orbit Float -> Angle Float -> Maybe (Angle Float)
eccentricAnomalyAtMeanAnomalyFloat :: Orbit Float -> Angle Float -> Maybe (Angle Float)
eccentricAnomalyAtMeanAnomalyFloat o :: Orbit Float
o _M :: Angle Float
_M = case Orbit Float -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit Float
o of
                                            Elliptic -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float)
forall a. a -> Maybe a
Just Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
Angle Float
_E
                                            _ -> Maybe (Angle Float)
forall a. Maybe a
Nothing
  where wrappedM :: Float
wrappedM = (Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
Angle Float
_M Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
forall a (u :: [Factor *]) (l :: LCSU *).
Real a =>
Qu u l a -> Qu u l a -> Qu u l a
`mod'` Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
forall a. Floating a => PlaneAngle a
turn) Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float -> Radian -> Float
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si|rad|]
        e :: Float
e = Orbit Float -> Unitless Float
forall a. Orbit a -> Unitless a
eccentricity Orbit Float
o Qu '[] 'DefaultLCSU Float -> Number -> Float
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si||]
        sinM :: Float
sinM = Float -> Float
forall a. Floating a => a -> a
sin Float
wrappedM
        cosM :: Float
cosM = Float -> Float
forall a. Floating a => a -> a
cos Float
wrappedM
        -- Use a better initial guess
        -- http://alpheratz.net/dynamics/twobody/KeplerIterations_summary.pdf
        initialGuess :: Float
initialGuess = Float
wrappedM Float -> Float -> Float
forall a. Num a => a -> a -> a
+
                       Float
e Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
sinM Float -> Float -> Float
forall a. Num a => a -> a -> a
+
                       Float
e Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
e Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
sinM Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
cosM Float -> Float -> Float
forall a. Num a => a -> a -> a
+
                       0.5 Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
e Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
e Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
e Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
sinM Float -> Float -> Float
forall a. Num a => a -> a -> a
* (3 Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
cosM Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
cosM Float -> Float -> Float
forall a. Num a => a -> a -> a
- 1)
        _E :: Angle Float
        _E :: Angle Float
_E = Float -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
forall a. Fractional a => a -> Angle a
rad (Float -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float)
-> ([Float] -> Float)
-> [Float]
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Float] -> Float
forall a. [a] -> a
last ([Float] -> Float) -> ([Float] -> [Float]) -> [Float] -> Float
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [Float] -> [Float]
forall a. Int -> [a] -> [a]
take 5 ([Float] -> Angle Float) -> [Float] -> Angle Float
forall a b. (a -> b) -> a -> b
$
             (forall s. AD s (Tower Float) -> AD s (Tower Float))
-> Float -> [Float]
forall a.
(Fractional a, Eq a) =>
(forall s. AD s (Tower a) -> AD s (Tower a)) -> a -> [a]
findZero (\_E :: AD s (Tower Float)
_E -> Scalar (AD s (Tower Float)) -> AD s (Tower Float)
forall t. Mode t => Scalar t -> t
auto Float
Scalar (AD s (Tower Float))
wrappedM AD s (Tower Float) -> AD s (Tower Float) -> AD s (Tower Float)
forall a. Num a => a -> a -> a
- (AD s (Tower Float)
_E AD s (Tower Float) -> AD s (Tower Float) -> AD s (Tower Float)
forall a. Num a => a -> a -> a
- Scalar (AD s (Tower Float)) -> AD s (Tower Float)
forall t. Mode t => Scalar t -> t
auto Float
Scalar (AD s (Tower Float))
e AD s (Tower Float) -> AD s (Tower Float) -> AD s (Tower Float)
forall a. Num a => a -> a -> a
* AD s (Tower Float) -> AD s (Tower Float)
forall a. Floating a => a -> a
sin AD s (Tower Float)
_E))
                      Float
initialGuess

-- | Calculate the eccentric anomaly, E, of an orbiting body when it has true
-- anomaly, ν.
--
-- The number of orbits represented by the anomalies is preserved;
-- i.e. ν `div` 2π = E `div` 2π
--
-- Returns Nothing if given a parabolic or hyperbolic orbit.
eccentricAnomalyAtTrueAnomaly :: (Floating a, Real a)
                              => Orbit a -> Angle a -> Maybe (Angle a)
eccentricAnomalyAtTrueAnomaly :: Orbit a -> Angle a -> Maybe (Angle a)
eccentricAnomalyAtTrueAnomaly o :: Orbit a
o ν :: Angle a
ν = case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
                                       Elliptic -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall a. a -> Maybe a
Just Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
_E
                                       _ -> Maybe (Angle a)
forall a. Maybe a
Nothing
  where (n :: Qu '[] 'DefaultLCSU Integer
n, wrappedν :: Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
wrappedν) = Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
ν Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> (Qu '[] 'DefaultLCSU Integer,
    Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall a b (u :: [Factor *]) (l :: LCSU *).
(Real a, Integral b) =>
Qu u l a -> Qu u l a -> (Qu '[] l b, Qu u l a)
`divMod'` Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
turn
        cosν :: a
cosν = a -> a
forall a. Floating a => a -> a
cos (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
ν Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Radian -> a
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si|rad|])
        -- sinν = sin (wrappedν # [si|rad|])
        e :: a
e = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o Qu '[] 'DefaultLCSU a -> Number -> a
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si||]
        wrappedE :: Angle a
wrappedE = a -> Angle a
forall a. Fractional a => a -> Angle a
rad (a -> Angle a) -> a -> Angle a
forall a b. (a -> b) -> a -> b
$ a -> a
forall a. Floating a => a -> a
acos ((a
e a -> a -> a
forall a. Num a => a -> a -> a
+ a
cosν) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
e a -> a -> a
forall a. Num a => a -> a -> a
* a
cosν))
        -- wrappedE = rad $ atan2 (sqrt (1 - e*e) * sinν) (e + cosν)
        _E :: Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
_E = if Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
wrappedν Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Bool
forall a. Ord a => a -> a -> Bool
< Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
halfTurn
               then ((Integer -> a)
-> Qu '[] 'DefaultLCSU Integer -> Qu '[] 'DefaultLCSU a
forall a b (u :: [Factor *]) (l :: LCSU *).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit Integer -> a
forall a. Num a => Integer -> a
fromInteger Qu '[] 'DefaultLCSU Integer
n Qu '[] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu (Normalize ('[] @+ '[ 'F PlaneAngle One])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
turn) Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
wrappedE
               else ((Integer -> a)
-> Qu '[] 'DefaultLCSU Integer -> Qu '[] 'DefaultLCSU a
forall a b (u :: [Factor *]) (l :: LCSU *).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit Integer -> a
forall a. Num a => Integer -> a
fromInteger (Qu '[] 'DefaultLCSU Integer
nQu '[] 'DefaultLCSU Integer
-> Qu '[] 'DefaultLCSU Integer -> Qu '[] 'DefaultLCSU Integer
forall a. Num a => a -> a -> a
+1) Qu '[] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu (Normalize ('[] @+ '[ 'F PlaneAngle One])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
turn) Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|-| Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
wrappedE

---------
-- To hyperbolic
---------

hyperbolicAnomalyAtTime
  :: forall a
   . (Converge [a], RealFloat a)
  => Orbit a
  -> Time a
  -> Maybe (AngleH a)
hyperbolicAnomalyAtTime :: Orbit a -> Time a -> Maybe (AngleH a)
hyperbolicAnomalyAtTime o :: Orbit a
o =
  Orbit a -> Angle a -> Maybe (AngleH a)
forall a.
(Converge [a], RealFloat a) =>
Orbit a -> Angle a -> Maybe (AngleH a)
hyperbolicAnomalyAtMeanAnomaly Orbit a
o (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
 -> Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a))
-> (Qu '[ 'F Time One] 'DefaultLCSU a
    -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> Qu '[ 'F Time One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Orbit a -> Time a -> Angle a
forall a. (Floating a, Ord a) => Orbit a -> Time a -> Angle a
meanAnomalyAtTime Orbit a
o

hyperbolicAnomalyAtMeanAnomaly
  :: forall a
   . (Converge [a], RealFloat a)
  => Orbit a
  -> Angle a
  -> Maybe (AngleH a)
hyperbolicAnomalyAtMeanAnomaly :: Orbit a -> Angle a -> Maybe (AngleH a)
hyperbolicAnomalyAtMeanAnomaly o :: Orbit a
o _M :: Angle a
_M = case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Hyperbolic -> Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a)
Maybe (AngleH a)
_H
  _          -> Maybe (AngleH a)
forall a. Maybe a
Nothing
 where
  e :: a
e                       = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o Qu '[] 'DefaultLCSU a -> Number -> a
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si||]
  _M' :: a
_M'                     = Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
_M Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Radian -> a
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si|rad|]
  _MDouble :: Double
_MDouble                = a -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac a
_M'
  Just initialGuessDouble :: Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
initialGuessDouble = Orbit Double -> Angle Double -> Maybe (AngleH Double)
hyperbolicAnomalyAtMeanAnomalyDouble
    ((a -> Double) -> Orbit a -> Orbit Double
forall a b. (a -> b) -> Orbit a -> Orbit b
unsafeMapOrbit a -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac Orbit a
o)
    (Double -> Angle Double
forall a. Fractional a => a -> Angle a
rad Double
_MDouble)
  initialGuess :: a
initialGuess = Double -> a
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Double -> a)
-> (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
    -> Double)
-> Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
-> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
-> RadianHyperbolic -> Double
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# RadianHyperbolic
RadianHyperbolic) (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double -> a)
-> Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double -> a
forall a b. (a -> b) -> a -> b
$ Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
initialGuessDouble
  err :: (Mode b, Floating b, Scalar b ~ a) => b -> b
  err :: b -> b
err _H :: b
_H = Scalar b -> b
forall t. Mode t => Scalar t -> t
auto a
Scalar b
_M' b -> b -> b
forall a. Num a => a -> a -> a
- (Scalar b -> b
forall t. Mode t => Scalar t -> t
auto a
Scalar b
e b -> b -> b
forall a. Num a => a -> a -> a
* b -> b
forall a. Floating a => a -> a
sinh b
_H b -> b -> b
forall a. Num a => a -> a -> a
- b
_H)
  _H :: Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a)
_H = (a -> Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a)
-> Maybe a
-> Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
forall a. Fractional a => a -> AngleH a
rdh (Maybe a
 -> Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a))
-> ([a] -> Maybe a)
-> [a]
-> Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Element [a] -> Element [a]) -> [a] -> Maybe (Element [a])
forall a.
(Converge a, Ord (Element a)) =>
(Element a -> Element a) -> a -> Maybe (Element a)
convergeErr (Id a -> a
forall a. Id a -> a
runId (Id a -> a) -> (a -> Id a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Id a -> Id a
forall a. Num a => a -> a
abs (Id a -> Id a) -> (a -> Id a) -> a -> Id a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Id a -> Id a
forall b. (Mode b, Floating b, Scalar b ~ a) => b -> b
err (Id a -> Id a) -> (a -> Id a) -> a -> Id a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> Id a
forall a. a -> Id a
Id) ([a] -> Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a))
-> [a] -> Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a)
forall a b. (a -> b) -> a -> b
$ (forall s. AD s (Tower a) -> AD s (Tower a)) -> a -> [a]
forall a.
Fractional a =>
(forall s. AD s (Tower a) -> AD s (Tower a)) -> a -> [a]
findZeroNoEq
    forall s. AD s (Tower a) -> AD s (Tower a)
forall b. (Mode b, Floating b, Scalar b ~ a) => b -> b
err
    a
initialGuess

-- | Calculate the hyperbolic anomaly, H, at a given mean anomaly. Unline
-- 'eccentricAnomalyAtMeanAnomalyFloat' this uses double precision floats to
-- help avoid overflowing.
hyperbolicAnomalyAtMeanAnomalyDouble
  :: Orbit Double -> Angle Double -> Maybe (AngleH Double)
hyperbolicAnomalyAtMeanAnomalyDouble :: Orbit Double -> Angle Double -> Maybe (AngleH Double)
hyperbolicAnomalyAtMeanAnomalyDouble o :: Orbit Double
o _M :: Angle Double
_M = case Orbit Double -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit Double
o of
  Hyperbolic -> case Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
_H of
    -- If you hit this, a better initial guess would probably help
    Qu x :: Double
x | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x -> String
-> Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double)
forall a. HasCallStack => String -> a
error "NaN while trying to find hyperbolic anomaly"
    _              -> Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
-> Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double)
forall a. a -> Maybe a
Just Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
_H
  _ -> Maybe (AngleH Double)
forall a. Maybe a
Nothing
 where
  -- Perhaps use something like https://www.researchgate.net/publication/226007277_A_Method_Solving_Kepler%27s_Equation_for_Hyperbolic_Case
  e :: Double
e            = Orbit Double -> Unitless Double
forall a. Orbit a -> Unitless a
eccentricity Orbit Double
o Qu '[] 'DefaultLCSU Double -> Number -> Double
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si||]
  _M' :: Double
_M'          = Qu '[ 'F PlaneAngle One] 'DefaultLCSU Double
Angle Double
_M Qu '[ 'F PlaneAngle One] 'DefaultLCSU Double -> Radian -> Double
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si|rad|]
  -- TODO: A better guess here
  initialGuess :: Double
initialGuess = Double
_M'
  _H :: Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
_H           = Double -> Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
forall a. Fractional a => a -> AngleH a
rdh (Double -> Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double)
-> ([Double] -> Double)
-> [Double]
-> Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Double
forall a. [a] -> a
last ([Double] -> Double)
-> ([Double] -> [Double]) -> [Double] -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take 200 ([Double]
 -> Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double)
-> [Double]
-> Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
forall a b. (a -> b) -> a -> b
$ (forall s. AD s ForwardDouble -> AD s ForwardDouble)
-> Double -> [Double]
Newton.findZero
    (\_H :: AD s ForwardDouble
_H -> Scalar (AD s ForwardDouble) -> AD s ForwardDouble
forall t. Mode t => Scalar t -> t
auto Double
Scalar (AD s ForwardDouble)
_M' AD s ForwardDouble -> AD s ForwardDouble -> AD s ForwardDouble
forall a. Num a => a -> a -> a
- (Scalar (AD s ForwardDouble) -> AD s ForwardDouble
forall t. Mode t => Scalar t -> t
auto Double
Scalar (AD s ForwardDouble)
e AD s ForwardDouble -> AD s ForwardDouble -> AD s ForwardDouble
forall a. Num a => a -> a -> a
* AD s ForwardDouble -> AD s ForwardDouble
forall a. Floating a => a -> a
sinh AD s ForwardDouble
_H AD s ForwardDouble -> AD s ForwardDouble -> AD s ForwardDouble
forall a. Num a => a -> a -> a
- AD s ForwardDouble
_H))
    Double
initialGuess

-- | Returns the hyperbolic anomaly, H, for an orbit at true anomaly ν.
--
-- Returns 'Nothing' when given an 'Elliptic' or 'Parabolic' orbit, or a true
-- anomaly out of the range of the hyperbolic orbit.
hyperbolicAnomalyAtTrueAnomaly
  :: (Floating a, Ord a) => Orbit a -> Angle a -> Maybe (AngleH a)
hyperbolicAnomalyAtTrueAnomaly :: Orbit a -> Angle a -> Maybe (AngleH a)
hyperbolicAnomalyAtTrueAnomaly o :: Orbit a
o ν :: Angle a
ν = case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  _ | Just d :: Angle a
d <- Orbit a -> Maybe (Angle a)
forall a. (Floating a, Ord a) => Orbit a -> Maybe (Angle a)
hyperbolicDepartureAngle Orbit a
o, Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a (l :: LCSU *) (u :: [Factor *]).
Num a =>
Qu u l a -> Qu u l a
qAbs Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
ν Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Bool
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Ord n) =>
Qu d1 l n -> Qu d2 l n -> Bool
|>| Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
d -> Maybe (AngleH a)
forall a. Maybe a
Nothing
  Hyperbolic -> Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a)
forall a. a -> Maybe a
Just Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
_H
  _          -> Maybe (AngleH a)
forall a. Maybe a
Nothing
 where
  e :: Unitless a
e     = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
  coshH :: Qu '[] 'DefaultLCSU a
coshH = (Angle a -> Unitless a
forall a. Floating a => Angle a -> Unitless a
qCos Angle a
ν Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Num a => a -> a -> a
+ Qu '[] 'DefaultLCSU a
Unitless a
e) Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Fractional a => a -> a -> a
/ (1 Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Num a => a -> a -> a
+ Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Num a => a -> a -> a
* Angle a -> Unitless a
forall a. Floating a => Angle a -> Unitless a
qCos Angle a
ν)
  sign :: a
sign  = a -> a
forall a. Num a => a -> a
signum (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
ν Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Radian -> a
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si|rad|])
  _H :: Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
_H    = a
sign a
-> Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
forall n (b :: [Factor *]) (l :: LCSU *).
Num n =>
n -> Qu b l n -> Qu b l n
*| Unitless a -> AngleH a
forall a. Floating a => Unitless a -> AngleH a
qArcCosh Qu '[] 'DefaultLCSU a
Unitless a
coshH

---------
-- To true
---------

-- | Calculate the true anomaly, ν, of a body at time since periapse, t.
trueAnomalyAtTime
  :: forall a . (Converge [a], RealFloat a) => Orbit a -> Time a -> Angle a
trueAnomalyAtTime :: Orbit a -> Time a -> Angle a
trueAnomalyAtTime o :: Orbit a
o t :: Time a
t = case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Elliptic   -> Orbit a -> Angle a -> Angle a
forall a.
(Converge [a], RealFloat a) =>
Orbit a -> Angle a -> Angle a
trueAnomalyAtMeanAnomaly Orbit a
o Angle a
_M
  Hyperbolic -> Orbit a -> Angle a -> Angle a
forall a.
(Converge [a], RealFloat a) =>
Orbit a -> Angle a -> Angle a
trueAnomalyAtMeanAnomaly Orbit a
o Angle a
_M
  Parabolic ->
    let _A :: Qu
  (Normalize ('[ 'F Time ('P 'Zero)] @+ '[ 'F Time One]))
  'DefaultLCSU
  a
_A = (3 Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[] @- '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| 2) Qu '[] 'DefaultLCSU a
-> Qu '[ 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Qu (Normalize ('[] @+ '[ 'F Time ('P 'Zero)])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[ 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu ('[ 'F Time ('P ('P 'Zero))] @/ 'S One) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu '[ 'F Length ('S ('S One))] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
         @- '[ 'F Length ('S ('S One))]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| (2 Qu '[] 'DefaultLCSU a
-> Qu '[ 'F Length ('S ('S One))] 'DefaultLCSU a
-> Qu
     (Normalize ('[] @+ '[ 'F Length ('S ('S One))])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu
     (Normalize
        (Normalize ('[ 'F Length One] @+ '[ 'F Length One])
         @+ '[ 'F Length One]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *).
Num n =>
Qu a l n -> Qu (Normalize (Normalize (a @+ a) @+ a)) l n
qCube (Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
l Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F Length One] @- '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| 2))) Qu '[ 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Qu '[ 'F Time One] 'DefaultLCSU a
-> Qu
     (Normalize ('[ 'F Time ('P 'Zero)] @+ '[ 'F Time One]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[ 'F Time One] 'DefaultLCSU a
Time a
t
        _B :: Qu ('[] @/ 'S ('S One)) 'DefaultLCSU a
_B = Qu '[] 'DefaultLCSU a -> Qu ('[] @/ 'S ('S One)) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *).
Floating n =>
Qu a l n -> Qu (a @/ 'S ('S One)) l n
qCubeRoot (Qu '[] 'DefaultLCSU a
Qu
  (Normalize ('[ 'F Time ('P 'Zero)] @+ '[ 'F Time One]))
  'DefaultLCSU
  a
_A Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| Qu '[] 'DefaultLCSU a -> Qu ('[] @/ 'S One) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Qu '[] 'DefaultLCSU a -> Qu (Normalize ('[] @+ '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *).
Num n =>
Qu a l n -> Qu (Normalize (a @+ a)) l n
qSq Qu '[] 'DefaultLCSU a
Qu
  (Normalize ('[ 'F Time ('P 'Zero)] @+ '[ 'F Time One]))
  'DefaultLCSU
  a
_A Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| 1))
    in  2 Qu '[] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu (Normalize ('[] @+ '[ 'F PlaneAngle One])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Unitless a -> Angle a
forall a. Floating a => Unitless a -> Angle a
qArcTan (Qu '[] 'DefaultLCSU a
Qu ('[] @/ 'S ('S One)) 'DefaultLCSU a
_B Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Num a => a -> a -> a
- Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Fractional a => a -> a
recip Qu '[] 'DefaultLCSU a
Qu ('[] @/ 'S ('S One)) 'DefaultLCSU a
_B)
 where
  μ :: Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ  = Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
forall a.
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
primaryGravitationalParameter Orbit a
o
  l :: Distance a
l  = Orbit a -> Distance a
forall a. Num a => Orbit a -> Distance a
semiLatusRectum Orbit a
o
  _M :: Angle a
_M = Orbit a -> Time a -> Angle a
forall a. (Floating a, Ord a) => Orbit a -> Time a -> Angle a
meanAnomalyAtTime Orbit a
o Time a
t

-- | Calculate the true anomaly, ν, of an orbiting body when it has the given
-- mean anomaly, _M.
trueAnomalyAtMeanAnomaly
  :: (Converge [a], RealFloat a) => Orbit a -> Angle a -> Angle a
trueAnomalyAtMeanAnomaly :: Orbit a -> Angle a -> Angle a
trueAnomalyAtMeanAnomaly o :: Orbit a
o _M :: Angle a
_M = case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Elliptic -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. HasCallStack => Maybe a -> a
fromJust
    (Orbit a -> Angle a -> Maybe (Angle a)
forall a. RealFloat a => Orbit a -> Angle a -> Maybe (Angle a)
trueAnomalyAtEccentricAnomaly Orbit a
o (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
 -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
    -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall (m :: * -> *) b c a.
Monad m =>
(b -> m c) -> (a -> m b) -> a -> m c
<=< Orbit a -> Angle a -> Maybe (Angle a)
forall a.
(Converge [a], Floating a, Real a) =>
Orbit a -> Angle a -> Maybe (Angle a)
eccentricAnomalyAtMeanAnomaly Orbit a
o (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
 -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall a b. (a -> b) -> a -> b
$ Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
_M)
  Hyperbolic -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. HasCallStack => Maybe a -> a
fromJust
    (Orbit a -> AngleH a -> Maybe (Angle a)
forall a.
(Ord a, Floating a) =>
Orbit a -> AngleH a -> Maybe (Angle a)
trueAnomalyAtHyperbolicAnomaly Orbit a
o (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
 -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
    -> Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a))
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall (m :: * -> *) b c a.
Monad m =>
(b -> m c) -> (a -> m b) -> a -> m c
<=< Orbit a -> Angle a -> Maybe (AngleH a)
forall a.
(Converge [a], RealFloat a) =>
Orbit a -> Angle a -> Maybe (AngleH a)
hyperbolicAnomalyAtMeanAnomaly Orbit a
o (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
 -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a))
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall a b. (a -> b) -> a -> b
$ Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
_M)
  _ -> String -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. HasCallStack => String -> a
error "trueAnomalyAtMeanAnomaly is not defined for Parabolic orbits"

-- | Calculate the true anomaly, ν, of an orbiting body when it has the given
-- eccentric anomaly, _E.
--
-- The number of orbits represented by the anomalies is preserved;
-- i.e. ν `div` 2π = E `div` 2π
trueAnomalyAtEccentricAnomaly :: RealFloat a
                              => Orbit a -- ^ An elliptic orbit
                              -> Angle a -- ^ The eccentric anomaly _E
                              -> Maybe (Angle a) -- ^ The true anomaly, ν
trueAnomalyAtEccentricAnomaly :: Orbit a -> Angle a -> Maybe (Angle a)
trueAnomalyAtEccentricAnomaly o :: Orbit a
o _E :: Angle a
_E = case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
                                       Elliptic -> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall a. a -> Maybe a
Just Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
ν
                                       _        -> Maybe (Angle a)
forall a. Maybe a
Nothing
  where (n :: Qu '[] 'DefaultLCSU a
n, wrappedE :: a
wrappedE) = (Qu '[] 'DefaultLCSU Integer -> Qu '[] 'DefaultLCSU a)
-> (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> a)
-> (Qu '[] 'DefaultLCSU Integer,
    Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> (Qu '[] 'DefaultLCSU a, a)
forall (p :: * -> * -> *) a b c d.
Bifunctor p =>
(a -> b) -> (c -> d) -> p a c -> p b d
bimap ((Integer -> a)
-> Qu '[] 'DefaultLCSU Integer -> Qu '[] 'DefaultLCSU a
forall a b (u :: [Factor *]) (l :: LCSU *).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit Integer -> a
forall a. Num a => Integer -> a
fromInteger) (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a -> Radian -> a
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si|rad|]) ((Qu '[] 'DefaultLCSU Integer,
  Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
 -> (Qu '[] 'DefaultLCSU a, a))
-> (Qu '[] 'DefaultLCSU Integer,
    Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
-> (Qu '[] 'DefaultLCSU a, a)
forall a b. (a -> b) -> a -> b
$
                        Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
_E Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> (Qu '[] 'DefaultLCSU Integer,
    Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
forall a b (u :: [Factor *]) (l :: LCSU *).
(Real a, Integral b) =>
Qu u l a -> Qu u l a -> (Qu '[] l b, Qu u l a)
`divMod'` Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
turn
        e :: a
e = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o Qu '[] 'DefaultLCSU a -> Number -> a
forall (dim :: [Factor *]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si||]
        wrappedν :: Angle a
wrappedν = a -> Angle a
forall a. Fractional a => a -> Angle a
rad (a -> Angle a) -> a -> Angle a
forall a b. (a -> b) -> a -> b
$ 2 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a -> a
forall a. RealFloat a => a -> a -> a
atan2 (a -> a
forall a. Floating a => a -> a
sqrt (1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
e) a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
sin (a
wrappedE a -> a -> a
forall a. Fractional a => a -> a -> a
/ 2))
                                   (a -> a
forall a. Floating a => a -> a
sqrt (1 a -> a -> a
forall a. Num a => a -> a -> a
- a
e) a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
cos (a
wrappedE a -> a -> a
forall a. Fractional a => a -> a -> a
/ 2))
        ν :: Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
ν = Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall a. Floating a => PlaneAngle a
turn Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[] 'DefaultLCSU a
n Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
Angle a
wrappedν

trueAnomalyAtHyperbolicAnomaly
  :: (Ord a, Floating a) => Orbit a -> AngleH a -> Maybe (Angle a)
trueAnomalyAtHyperbolicAnomaly :: Orbit a -> AngleH a -> Maybe (Angle a)
trueAnomalyAtHyperbolicAnomaly o :: Orbit a
o _H :: AngleH a
_H = case Orbit a -> Classification
forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Hyperbolic -> Qu ('F PlaneAngle One : Normalize '[]) 'DefaultLCSU a
-> Maybe (Qu ('F PlaneAngle One : Normalize '[]) 'DefaultLCSU a)
forall a. a -> Maybe a
Just Qu ('F PlaneAngle One : Normalize '[]) 'DefaultLCSU a
Qu (Normalize ('[] @+ '[ 'F PlaneAngle One])) 'DefaultLCSU a
ν
  _          -> Maybe (Angle a)
forall a. Maybe a
Nothing
 where
  e :: Unitless a
e         = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
  tanνOver2 :: Qu '[] 'DefaultLCSU a
tanνOver2 = Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Floating a => a -> a
sqrt ((Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Num a => a -> a -> a
+ 1) Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Fractional a => a -> a -> a
/ (Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Num a => a -> a -> a
- 1)) Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall a. Num a => a -> a -> a
* AngleH a -> Unitless a
forall a. Floating a => AngleH a -> Unitless a
qTanh (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
AngleH a
_H Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu
     (Normalize ('[ 'F PlaneAngleHyperbolic One] @- '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| 2)
  ν :: Qu (Normalize ('[] @+ '[ 'F PlaneAngle One])) 'DefaultLCSU a
ν         = 2 Qu '[] 'DefaultLCSU a
-> Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
-> Qu (Normalize ('[] @+ '[ 'F PlaneAngle One])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Unitless a -> Angle a
forall a. Floating a => Unitless a -> Angle a
qArcTan Qu '[] 'DefaultLCSU a
Unitless a
tanνOver2

----------------------------------------------------------------
-- Other orbital properties
----------------------------------------------------------------

-- | The distance, r, from the primary body to the orbiting body at a particular
-- true anomaly.
radiusAtTrueAnomaly :: (Ord a, Floating a) => Orbit a -> Angle a -> Distance a
radiusAtTrueAnomaly :: Orbit a -> Angle a -> Distance a
radiusAtTrueAnomaly o :: Orbit a
o trueAnomaly :: Angle a
trueAnomaly = case Orbit a -> Maybe (Distance a)
forall a. (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)
semiMajorAxis Orbit a
o of
  Just _  -> Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
l Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F Length One] @- '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| (1 Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| Qu '[] 'DefaultLCSU a
Unitless a
e Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[] @+ '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Angle a -> Unitless a
forall a. Floating a => Angle a -> Unitless a
qCos Angle a
ν)
  Nothing -> (Qu '[ 'F Length ('S One), 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S One), 'F Time ('P 'Zero)]
         @+ '[ 'F Length ('S One), 'F Time ('P 'Zero)]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *).
Num n =>
Qu a l n -> Qu (Normalize (a @+ a)) l n
qSq Qu '[ 'F Length ('S One), 'F Time ('P 'Zero)] 'DefaultLCSU a
Quantity ((Meter :^ Succ (Succ 'Zero)) :* (Second :^ Pred 'Zero)) a
h Qu
  '[ 'F Length ('S ('S ('S One))), 'F Time ('P ('P 'Zero))]
  'DefaultLCSU
  a
-> Qu
     '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S ('S ('S One))), 'F Time ('P ('P 'Zero))]
         @- '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ) Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[ 'F Length One] @+ '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| (1 Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu (Normalize ('[] @- '[])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| (1 Qu '[] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a -> Qu '[] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| Angle a -> Unitless a
forall a. Floating a => Angle a -> Unitless a
qCos Angle a
ν))
 where
  h :: Quantity ((Meter :^ Succ (Succ 'Zero)) :* (Second :^ Pred 'Zero)) a
h = Orbit a
-> Quantity
     ((Meter :^ Succ (Succ 'Zero)) :* (Second :^ Pred 'Zero)) a
forall a.
Floating a =>
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ 'Zero)) :* (Second :^ Pred 'Zero)) a
specificAngularMomentum Orbit a
o
  e :: Unitless a
e = Orbit a -> Unitless a
forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
  ν :: Angle a
ν = Angle a
trueAnomaly
  μ :: Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ = Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
forall a.
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
primaryGravitationalParameter Orbit a
o
  l :: Distance a
l = Orbit a -> Distance a
forall a. Num a => Orbit a -> Distance a
semiLatusRectum Orbit a
o

-- | What is the speed, v, of a body at a particular true anomaly
speedAtTrueAnomaly :: (Ord a, Floating a) => Orbit a -> Angle a -> Speed a
speedAtTrueAnomaly :: Orbit a -> Angle a -> Speed a
speedAtTrueAnomaly o :: Orbit a
o trueAnomaly :: Angle a
trueAnomaly = case Orbit a -> Maybe (Distance a)
forall a. (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)
semiMajorAxis Orbit a
o of
  Nothing -> Qu '[ 'F Length ('S One), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu
     ('[ 'F Length ('S One), 'F Time ('P ('P 'Zero))] @/ 'S One)
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] @+ '[]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| 2 Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
         @- '[ 'F Length One]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
r)
  Just a :: Distance a
a  -> Qu '[ 'F Length ('S One), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu
     ('[ 'F Length ('S One), 'F Time ('P ('P 'Zero))] @/ 'S One)
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu '[ 'F Length ('P 'Zero)] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
         @+ '[ 'F Length ('P 'Zero)]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| (2 Qu '[] 'DefaultLCSU a
-> Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu (Normalize ('[] @- '[ 'F Length One])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
r Qu '[ 'F Length ('P 'Zero)] 'DefaultLCSU a
-> Qu '[ 'F Length ('P 'Zero)] 'DefaultLCSU a
-> Qu '[ 'F Length ('P 'Zero)] 'DefaultLCSU a
forall (d1 :: [Factor *]) (d2 :: [Factor *]) n (l :: LCSU *).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|-| 1 Qu '[] 'DefaultLCSU a
-> Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu (Normalize ('[] @- '[ 'F Length One])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
a))
 where
  ν :: Angle a
ν = Angle a
trueAnomaly
  μ :: Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ = Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
forall a.
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
primaryGravitationalParameter Orbit a
o
  r :: Distance a
r = Orbit a -> Angle a -> Distance a
forall a. (Ord a, Floating a) => Orbit a -> Angle a -> Distance a
radiusAtTrueAnomaly Orbit a
o Angle a
ν

-- | Specific angular momentum, h, is the angular momentum per unit mass
specificAngularMomentum :: Floating a => Orbit a -> Quantity [si|m^2 s^-1|] a
specificAngularMomentum :: Orbit a
-> Quantity
     ((Meter :^ Succ (Succ 'Zero)) :* (Second :^ Pred 'Zero)) a
specificAngularMomentum o :: Orbit a
o = Qu
  '[ 'F Length ('S ('S ('S One))), 'F Time ('P ('P 'Zero))]
  'DefaultLCSU
  a
-> Qu
     ('[ 'F Length ('S ('S ('S One))), 'F Time ('P ('P 'Zero))]
      @/ 'S One)
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
         @+ '[ 'F Length One]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
l)
  where
    μ :: Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ = Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
forall a.
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
primaryGravitationalParameter Orbit a
o
    l :: Distance a
l = Orbit a -> Distance a
forall a. Num a => Orbit a -> Distance a
semiLatusRectum Orbit a
o

-- | Specific orbital energy, ε, is the orbital energy per unit mass
specificOrbitalEnergy
  :: (Ord a, Floating a) => Orbit a -> Quantity [si|J / kg|] a
specificOrbitalEnergy :: Orbit a -> Quantity (Joule :/ (Kilo :@ Gram)) a
specificOrbitalEnergy o :: Orbit a
o = case Orbit a -> Maybe (Distance a)
forall a. (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)
semiMajorAxis Orbit a
o of
  Just a :: Distance a
a  -> Qu '[ 'F Length ('S One), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu
     '[ 'F Length ('S One), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
forall n (d :: [Factor *]) (l :: LCSU *).
Num n =>
Qu d l n -> Qu d l n
qNegate (Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
         @- '[ 'F Length One]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| (2 Qu '[] 'DefaultLCSU a
-> Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu (Normalize ('[] @+ '[ 'F Length One])) 'DefaultLCSU a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
a))
  Nothing -> Quantity (Joule :/ (Kilo :@ Gram)) a
forall n (dimspec :: [Factor *]) (l :: LCSU *).
Num n =>
Qu dimspec l n
zero
  where μ :: Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ = Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
forall a.
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
primaryGravitationalParameter Orbit a
o

-- | Specific potential energy, εp, is the potential energy per unit mass at a
-- particular true anomaly
specificPotentialEnergyAtTrueAnomaly
  :: (Ord a, Floating a) => Orbit a -> Angle a -> Quantity [si|J / kg|] a
specificPotentialEnergyAtTrueAnomaly :: Orbit a -> Angle a -> Quantity (Joule :/ (Kilo :@ Gram)) a
specificPotentialEnergyAtTrueAnomaly o :: Orbit a
o ν :: Angle a
ν = Qu '[ 'F Length ('S One), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu
     '[ 'F Length ('S One), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
forall n (d :: [Factor *]) (l :: LCSU *).
Num n =>
Qu d l n -> Qu d l n
qNegate (Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
         @- '[ 'F Length One]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
r)
 where
  r :: Distance a
r = Orbit a -> Angle a -> Distance a
forall a. (Ord a, Floating a) => Orbit a -> Angle a -> Distance a
radiusAtTrueAnomaly Orbit a
o Angle a
ν
  μ :: Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ = Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
forall a.
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
primaryGravitationalParameter Orbit a
o

-- | Specific kinetic energy, εk, is the kinetic energy per unit mass at a
-- particular true anomaly
specificKineticEnergyAtTrueAnomaly
  :: (Ord a, Floating a) => Orbit a -> Angle a -> Quantity [si|J / kg|] a
specificKineticEnergyAtTrueAnomaly :: Orbit a -> Angle a -> Quantity (Joule :/ (Kilo :@ Gram)) a
specificKineticEnergyAtTrueAnomaly o :: Orbit a
o ν :: Angle a
ν = Qu '[ 'F Length One, 'F Time ('P 'Zero)] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length One, 'F Time ('P 'Zero)]
         @+ '[ 'F Length One, 'F Time ('P 'Zero)]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *).
Num n =>
Qu a l n -> Qu (Normalize (a @+ a)) l n
qSq (Orbit a -> Angle a -> Speed a
forall a. (Ord a, Floating a) => Orbit a -> Angle a -> Speed a
speedAtTrueAnomaly Orbit a
o Angle a
ν) Qu '[ 'F Length ('S One), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu '[] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S One), 'F Time ('P ('P 'Zero))] @- '[]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| 2

----------------------------------------------------------------
-- Utils
----------------------------------------------------------------

-- | The escape velocity for a primary with specified gravitational parameter
-- at a particular distance.
escapeVelocityAtDistance
  :: (Floating a) => Quantity [si| m^3 s^-2 |] a -> Distance a -> Speed a
escapeVelocityAtDistance :: Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
-> Distance a -> Speed a
escapeVelocityAtDistance μ :: Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ r :: Distance a
r = Qu '[ 'F Length ('S One), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu
     ('[ 'F Length ('S One), 'F Time ('P ('P 'Zero))] @/ 'S One)
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (2 Qu '[] 'DefaultLCSU a
-> Qu
     '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[] @+ '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ Qu
  '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU a
-> Qu '[ 'F Length One] 'DefaultLCSU a
-> Qu
     (Normalize
        ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
         @- '[ 'F Length One]))
     'DefaultLCSU
     a
forall n (a :: [Factor *]) (l :: LCSU *) (b :: [Factor *]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Qu '[ 'F Length One] 'DefaultLCSU a
Distance a
r)