{-# language CPP #-}
{-# 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
                                                )
#if MIN_VERSION_ad(4,5,0)
import           Numeric.AD.Rank1.Halley        ( findZero
                                                , findZeroNoEq
                                                )
#else
import           Numeric.AD.Halley              ( findZero
                                                , findZeroNoEq
                                                )
#endif

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.
                       forall a. 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.
                     , forall a. Orbit a -> Distance a
periapsis                     :: !(Distance a)
                       -- | The 'inclinationSpecifier' describes the angle
                       -- between the obtital plane and the reference plane.
                     , forall a. 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.
                     , forall a. 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.
                     , forall a.
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
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
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.
                                         forall a. InclinationSpecifier a -> Angle a
longitudeOfAscendingNode :: !(Angle a)
                                         -- | The orbit's inclination, i.
                                         --
                                         -- The angle between the reference
                                         -- plane and the orbital plane
                                       , forall a. InclinationSpecifier a -> Angle a
inclination              :: !(Angle a)
                                       }
                              -- | The orbit lies in the reference plane
                            | NonInclined
  deriving (Int -> InclinationSpecifier a -> ShowS
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
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.
                                        forall a. 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
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
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
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]
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
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 :: forall a b (u :: [Factor (*)]) (l :: LCSU (*)).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit a -> b
f = forall (d :: [Factor (*)]) (l :: LCSU (*)) n.
UnsafeQu d l n -> Qu d l n
qu forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> b
f forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (d :: [Factor (*)]) (l :: LCSU (*)) n.
Qu d l n -> UnsafeQu d l n
UnsafeQu

unsafeMapOrbit :: (a -> b) -> Orbit a -> Orbit b
unsafeMapOrbit :: forall a b. (a -> b) -> Orbit a -> Orbit b
unsafeMapOrbit a -> b
f (Orbit Unitless a
e Distance a
q InclinationSpecifier a
i PeriapsisSpecifier a
p Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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 (forall a b (u :: [Factor (*)]) (l :: LCSU (*)).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit a -> b
f Unitless a
e)
                                           (forall a b (u :: [Factor (*)]) (l :: LCSU (*)).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit a -> b
f Distance a
q)
                                           (forall a b.
(a -> b) -> InclinationSpecifier a -> InclinationSpecifier b
unsafeMapInclinationSpecifier a -> b
f InclinationSpecifier a
i)
                                           (forall a b.
(a -> b) -> PeriapsisSpecifier a -> PeriapsisSpecifier b
unsafeMapPeriapsisSpecifier a -> b
f PeriapsisSpecifier a
p)
                                           (forall a b (u :: [Factor (*)]) (l :: LCSU (*)).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit a -> b
f Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ)

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

unsafeMapPeriapsisSpecifier :: (a -> b)
                            -> PeriapsisSpecifier a -> PeriapsisSpecifier b
unsafeMapPeriapsisSpecifier :: forall a b.
(a -> b) -> PeriapsisSpecifier a -> PeriapsisSpecifier b
unsafeMapPeriapsisSpecifier a -> b
f PeriapsisSpecifier a
p = case PeriapsisSpecifier a
p of
  PeriapsisSpecifier a
Circular    -> forall a. PeriapsisSpecifier a
Circular
  Eccentric Angle a
a -> forall a. Angle a -> PeriapsisSpecifier a
Eccentric (forall a b (u :: [Factor (*)]) (l :: LCSU (*)).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit a -> b
f 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 :: forall a. (Ord a, Num a) => Orbit a -> Bool
isValid Orbit a
o = Unitless a
e forall a. Ord a => a -> a -> Bool
>= Qu '[] 'DefaultLCSU a
0 Bool -> Bool -> Bool
&&
            ((Unitless a
e forall a. Eq a => a -> a -> Bool
== Qu '[] 'DefaultLCSU a
0) Bool -> Bool -> Bool
`iff` (forall a. Orbit a -> PeriapsisSpecifier a
periapsisSpecifier Orbit a
o forall a. Eq a => a -> a -> Bool
== forall a. PeriapsisSpecifier a
Circular)) Bool -> Bool -> Bool
&&
            Distance a
q forall a. Ord a => a -> a -> Bool
> forall n (dimspec :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu dimspec l n
zero Bool -> Bool -> Bool
&&
            Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  a
μ forall a. Ord a => a -> a -> Bool
> forall n (dimspec :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu dimspec l n
zero
  where
    iff :: Bool -> Bool -> Bool
iff = forall a. Eq a => a -> a -> Bool
(==) :: Bool -> Bool -> Bool
    e :: Unitless a
e = forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
    q :: Distance a
q = forall a. Orbit a -> Distance a
periapsis Orbit a
o
    μ :: 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 :: forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o | Unitless a
e forall a. Ord a => a -> a -> Bool
< Qu '[] 'DefaultLCSU a
1     = Classification
Elliptic
           | Unitless a
e forall a. Eq a => a -> a -> Bool
== Qu '[] 'DefaultLCSU a
1    = Classification
Parabolic
           | Unitless a
e forall a. Ord a => a -> a -> Bool
> Qu '[] 'DefaultLCSU a
1     = Classification
Hyperbolic
           | Bool
otherwise = forall a. HasCallStack => String -> a
error String
"classify: NaN eccentricity"
  where e :: Unitless a
e = 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 :: forall a. (Floating a, Real a) => Orbit a -> Orbit a
normalizeOrbit (Orbit Unitless a
e Distance a
q InclinationSpecifier a
inc PeriapsisSpecifier a
per Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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
  (InclinationSpecifier a
inc', Bool
flipped) = case InclinationSpecifier a
inc of
    InclinationSpecifier a
NonInclined              -> (forall a. InclinationSpecifier a
NonInclined, Bool
False)
    Inclined Angle a
_ Angle a
i | Angle a
i forall a. Eq a => a -> a -> Bool
== forall n (dimspec :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu dimspec l n
zero -> (forall a. InclinationSpecifier a
NonInclined, Bool
False)
    Inclined Angle a
_Ω Angle a
i ->
      let iR :: Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
iR  = Angle a
i forall a (u :: [Factor (*)]) (l :: LCSU (*)).
Real a =>
Qu u l a -> Qu u l a -> Qu u l a
`mod'` forall a. Floating a => PlaneAngle a
turn
          i' :: Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
i'  = if Bool
flipped then forall a. Floating a => PlaneAngle a
turn 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 Angle 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
|+| forall a. Floating a => PlaneAngle a
halfTurn else Angle a
_Ω) forall a (u :: [Factor (*)]) (l :: LCSU (*)).
Real a =>
Qu u l a -> Qu u l a -> Qu u l a
`mod'` forall a. Floating a => PlaneAngle a
turn
      in  (forall a. Angle a -> Angle a -> InclinationSpecifier a
Inclined Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
_Ω' Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
i', Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
iR forall a. Ord a => a -> a -> Bool
>= forall a. Floating a => PlaneAngle a
halfTurn)

  per' :: PeriapsisSpecifier a
per' = case PeriapsisSpecifier a
per of
    PeriapsisSpecifier a
Circular | Bool
flipped   -> forall a. Angle a -> PeriapsisSpecifier a
Eccentric forall a. Floating a => PlaneAngle a
halfTurn
             | Bool
otherwise -> forall a. PeriapsisSpecifier a
Circular
    Eccentric Angle a
ω
      | Angle a
ω forall a. Eq a => a -> a -> Bool
== forall n (dimspec :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu dimspec l n
zero, Unitless a
e forall a. Eq a => a -> a -> Bool
== Qu '[] 'DefaultLCSU a
0, Bool -> Bool
not Bool
flipped
      -> forall a. PeriapsisSpecifier a
Circular
      | Bool
otherwise
      -> forall a. Angle a -> PeriapsisSpecifier a
Eccentric forall a b. (a -> b) -> a -> b
$ (if Bool
flipped then Angle 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
|+| forall a. Floating a => PlaneAngle a
halfTurn else Angle a
ω) forall a (u :: [Factor (*)]) (l :: LCSU (*)).
Real a =>
Qu u l a -> Qu u l a -> Qu u l a
`mod'` 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 :: forall a. (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)
semiMajorAxis Orbit a
o =
  case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
    Classification
Parabolic -> forall a. Maybe a
Nothing
    Classification
_         -> forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ Distance a
q 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 '[] 'DefaultLCSU a
1 forall (d1 :: [Factor (*)]) (d2 :: [Factor (*)]) n (l :: LCSU (*)).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|-| Unitless a
e)
  where
    q :: Distance a
q = forall a. Orbit a -> Distance a
periapsis Orbit a
o
    e :: Unitless a
e = 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 :: forall a. (Floating a, Ord a) => Orbit a -> Distance a
semiMinorAxis Orbit a
o =
  case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
    Classification
Elliptic   -> Distance a
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
|*| forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Qu '[] 'DefaultLCSU a
1 forall (d1 :: [Factor (*)]) (d2 :: [Factor (*)]) n (l :: LCSU (*)).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|-| Unitless a
e forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2::Int))
    Classification
Parabolic  -> forall n (dimspec :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu dimspec l n
zero
    Classification
Hyperbolic -> Distance a
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
|*| forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Unitless a
e forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2::Int) 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
1)
  where
    e :: Unitless a
e = forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
    Just Distance a
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 :: forall a. Num a => Orbit a -> Distance a
semiLatusRectum Orbit a
orbit = Unitless a
e forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Distance a
q forall (d1 :: [Factor (*)]) (d2 :: [Factor (*)]) n (l :: LCSU (*)).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| Distance a
q
  where q :: Distance a
q = forall a. Orbit a -> Distance a
periapsis Orbit a
orbit
        e :: Unitless a
e = 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 :: forall a. (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)
apoapsis Orbit a
o =
  case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
    Classification
Elliptic -> forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ Distance a
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
1 forall (d1 :: [Factor (*)]) (d2 :: [Factor (*)]) n (l :: LCSU (*)).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| Unitless a
e)
    Classification
_        -> forall a. Maybe a
Nothing
  where
    Just Distance a
a = forall a. (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)
semiMajorAxis Orbit a
o
    e :: Unitless a
e = 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 :: forall a.
(Floating a, Ord a) =>
Orbit a -> Quantity (Radian :/ Second) a
meanMotion Orbit a
o =
  case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
    Classification
Elliptic   -> forall (b :: [Factor (*)]) a.
Qu b 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ b)) 'DefaultLCSU a
addRad forall a b. (a -> b) -> a -> b
$ forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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
|/| forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu a l n -> Qu (Normalize (Normalize (a @+ a) @+ a)) l n
qCube Distance a
a)
    Classification
Hyperbolic -> forall (b :: [Factor (*)]) a.
Qu b 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ b)) 'DefaultLCSU a
addRad forall a b. (a -> b) -> a -> b
$ forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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
|/| forall n (d :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu d l n -> Qu d l n
qNegate (forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu a l n -> Qu (Normalize (Normalize (a @+ a) @+ a)) l n
qCube Distance a
a))
    Classification
Parabolic  -> forall (b :: [Factor (*)]) a.
Qu b 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ b)) 'DefaultLCSU a
addRad forall a b. (a -> b) -> a -> b
$ Qu '[] 'DefaultLCSU a
2 forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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
|/| forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu a l n -> Qu (Normalize (Normalize (a @+ a) @+ a)) l n
qCube Distance a
l)
  where
    Just Distance a
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
μ = forall a.
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
primaryGravitationalParameter Orbit a
o
    l :: Distance a
l = 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 :: forall a. (Floating a, Ord a) => Orbit a -> Maybe (Time a)
period Orbit a
o =
  case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
    Classification
Elliptic -> forall a. a -> Maybe a
Just Qu
  (Normalize
     ('[ 'F PlaneAngle One]
      @- '[ 'F PlaneAngle One, 'F Time ('P 'Zero)]))
  'DefaultLCSU
  a
p
    Classification
_ -> forall a. Maybe a
Nothing
  where
    n :: Quantity (Radian :/ Second) a
n = 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 = forall a. Floating a => PlaneAngle a
turn forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| 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 :: forall a.
(Ord a, Floating a) =>
Orbit a -> Quantity ((Meter :^ Succ (Succ 'Zero)) :/ Second) a
arealVelocity Orbit a
o = forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Distance a
l forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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 '[] 'DefaultLCSU a
2
  where l :: Distance a
l = forall a. Num a => Orbit a -> Distance a
semiLatusRectum Orbit a
o
        μ :: 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 :: forall a. (Floating a, Ord a) => Orbit a -> Maybe (Angle a)
hyperbolicDepartureAngle Orbit a
o =
  case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
    Classification
Hyperbolic ->
      let e :: Unitless a
e = forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
          θ :: Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
θ = forall (b :: [Factor (*)]) a.
Qu b 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ b)) 'DefaultLCSU a
addRad forall a b. (a -> b) -> a -> b
$ forall a. Floating a => a -> a
acos (-Qu '[] 'DefaultLCSU a
1 forall a. Fractional a => a -> a -> a
/ Unitless a
e)
      in forall a. a -> Maybe a
Just Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
θ
    Classification
Parabolic -> forall a. a -> Maybe a
Just (forall a. Floating a => PlaneAngle a
turn 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 '[] 'DefaultLCSU a
2)
    Classification
_ -> 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 :: forall a. (Floating a, Ord a) => Orbit a -> Maybe (Angle a)
hyperbolicApproachAngle = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall n (d :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu d l n -> Qu d l n
qNegate forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 :: forall a. (Floating a, Ord a) => Orbit a -> Angle a -> Time a
timeAtMeanAnomaly Orbit a
o Angle a
_M = Angle a
_M forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Quantity (Radian :/ Second) a
n
  where n :: Quantity (Radian :/ Second) a
n = 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 :: forall a.
(Floating a, Ord a) =>
Orbit a -> Angle a -> Maybe (Time a)
timeAtEccentricAnomaly Orbit a
o = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. (Floating a, Ord a) => Orbit a -> Angle a -> Time a
timeAtMeanAnomaly Orbit a
o) forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 :: forall a.
(Floating a, Ord a) =>
Orbit a -> AngleH a -> Maybe (Time a)
timeAtHyperbolicAnomaly Orbit a
o =
  forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. (Floating a, Ord a) => Orbit a -> Angle a -> Time a
timeAtMeanAnomaly Orbit a
o) forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 :: forall a.
(Real a, Floating a) =>
Orbit a -> Angle a -> Maybe (Time a)
timeAtTrueAnomaly Orbit a
o Angle a
ν = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Classification
_ | Just Angle a
d <- forall a. (Floating a, Ord a) => Orbit a -> Maybe (Angle a)
hyperbolicDepartureAngle Orbit a
o, forall a (l :: LCSU (*)) (u :: [Factor (*)]).
Num a =>
Qu u l a -> Qu u l a
qAbs Angle a
ν forall (d1 :: [Factor (*)]) (d2 :: [Factor (*)]) n (l :: LCSU (*)).
(d1 @~ d2, Ord n) =>
Qu d1 l n -> Qu d2 l n -> Bool
|>| Angle a
d -> forall a. Maybe a
Nothing
  Classification
Parabolic ->
    let _D :: Unitless a
_D = forall a. Floating a => Angle a -> Unitless a
qTan (Angle 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 '[] 'DefaultLCSU a
2)
        t :: Qu
  (Normalize
     (Normalize
        ('[]
         @@+ Reorder
               (Normalize
                  (Normalize
                     (Normalize ('[ 'F Length One] @@+ '[ 'F Length One])
                      @@+ Reorder
                            '[ 'F Length One]
                            (Normalize ('[ 'F Length One] @@+ '[ 'F Length One])))
                   @- '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))])
                @/ 'S One)
               '[])
      @+ '[]))
  'DefaultLCSU
  a
t  = Qu '[] 'DefaultLCSU a
0.5 forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu a l n -> Qu (Normalize (Normalize (a @+ a) @+ a)) l n
qCube Distance a
l forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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
_D forall (d1 :: [Factor (*)]) (d2 :: [Factor (*)]) n (l :: LCSU (*)).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| (forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu a l n -> Qu (Normalize (Normalize (a @+ a) @+ a)) l n
qCube Unitless a
_D 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 '[] 'DefaultLCSU a
3))
    in  forall a. a -> Maybe a
Just Qu
  (Normalize
     (Normalize
        ('[]
         @@+ Reorder
               (Normalize
                  (Normalize
                     (Normalize ('[ 'F Length One] @@+ '[ 'F Length One])
                      @@+ Reorder
                            '[ 'F Length One]
                            (Normalize ('[ 'F Length One] @@+ '[ 'F Length One])))
                   @- '[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))])
                @/ 'S One)
               '[])
      @+ '[]))
  'DefaultLCSU
  a
t
  Classification
_ -> forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. (Floating a, Ord a) => Orbit a -> Angle a -> Time a
timeAtMeanAnomaly Orbit a
o) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a.
(Real a, Floating a) =>
Orbit a -> Angle a -> Maybe (Angle a)
meanAnomalyAtTrueAnomaly Orbit a
o forall a b. (a -> b) -> a -> b
$ Angle a
ν
 where
  μ :: 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 = 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 :: forall a. (Floating a, Ord a) => Orbit a -> Time a -> Angle a
meanAnomalyAtTime Orbit a
o Time a
t = Time a
t forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Quantity (Radian :/ Second) a
n
  where n :: Quantity (Radian :/ Second) a
n = 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 :: forall a.
(Floating a, Ord a) =>
Orbit a -> Angle a -> Maybe (Angle a)
meanAnomalyAtEccentricAnomaly Orbit a
o Angle a
_E = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
                                       Classification
Elliptic -> forall a. a -> Maybe a
Just Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
_M
                                       Classification
_ -> forall a. Maybe a
Nothing
  where e :: Unitless a
e = forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
        untypedE :: Qu
  (Normalize ('[ 'F PlaneAngle One] @- '[ 'F PlaneAngle One]))
  'DefaultLCSU
  a
untypedE = forall (u :: [Factor (*)]) a.
Qu u 'DefaultLCSU a
-> Qu (Normalize (u @- '[ 'F PlaneAngle One])) 'DefaultLCSU a
delRad Angle a
_E
        _M :: Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
_M = forall (b :: [Factor (*)]) a.
Qu b 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ b)) 'DefaultLCSU a
addRad (Qu
  (Normalize ('[ 'F PlaneAngle One] @- '[ 'F PlaneAngle One]))
  'DefaultLCSU
  a
untypedE forall (d1 :: [Factor (*)]) (d2 :: [Factor (*)]) n (l :: LCSU (*)).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|-| Unitless a
e forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| forall a. Floating a => a -> a
sin 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 :: forall a.
(Floating a, Ord a) =>
Orbit a -> AngleH a -> Maybe (Angle a)
meanAnomalyAtHyperbolicAnomaly Orbit a
o AngleH a
_H = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Classification
Hyperbolic -> forall a. a -> Maybe a
Just Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
_M
  Classification
_          -> forall a. Maybe a
Nothing
 where
  e :: Unitless a
e  = forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
  _M :: Qu (Normalize ('[ 'F PlaneAngle One] @+ '[])) 'DefaultLCSU a
_M = forall (b :: [Factor (*)]) a.
Qu b 'DefaultLCSU a
-> Qu (Normalize ('[ 'F PlaneAngle One] @+ b)) 'DefaultLCSU a
addRad forall a b. (a -> b) -> a -> b
$ Unitless a
e forall a. Num a => a -> a -> a
* forall a. Floating a => AngleH a -> Unitless a
qSinh AngleH a
_H forall a. Num a => a -> a -> a
- forall n (l :: LCSU (*)). n -> Qu '[] l n
quantity (AngleH a
_H 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 :: forall a.
(Real a, Floating a) =>
Orbit a -> Angle a -> Maybe (Angle a)
meanAnomalyAtTrueAnomaly Orbit a
o Angle a
ν = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Classification
Parabolic -> forall a. Maybe a
Nothing
  Classification
Elliptic -> forall a.
(Floating a, Ord a) =>
Orbit a -> Angle a -> Maybe (Angle a)
meanAnomalyAtEccentricAnomaly Orbit a
o forall (m :: * -> *) b c a.
Monad m =>
(b -> m c) -> (a -> m b) -> a -> m c
<=<
              forall a.
(Floating a, Real a) =>
Orbit a -> Angle a -> Maybe (Angle a)
eccentricAnomalyAtTrueAnomaly Orbit a
o forall a b. (a -> b) -> a -> b
$ Angle a
ν
  Classification
Hyperbolic -> forall a.
(Floating a, Ord a) =>
Orbit a -> AngleH a -> Maybe (Angle a)
meanAnomalyAtHyperbolicAnomaly Orbit a
o forall (m :: * -> *) b c a.
Monad m =>
(b -> m c) -> (a -> m b) -> a -> m c
<=<
                forall a.
(Floating a, Ord a) =>
Orbit a -> Angle a -> Maybe (AngleH a)
hyperbolicAnomalyAtTrueAnomaly Orbit a
o forall a b. (a -> b) -> a -> b
$ 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 :: forall a.
(Converge [a], Floating a, Real a) =>
Orbit a -> Time a -> Maybe (Angle a)
eccentricAnomalyAtTime Orbit a
o Time a
t = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Classification
Elliptic -> forall a.
(Converge [a], Floating a, Real a) =>
Orbit a -> Angle a -> Maybe (Angle a)
eccentricAnomalyAtMeanAnomaly Orbit a
o forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (Floating a, Ord a) => Orbit a -> Time a -> Angle a
meanAnomalyAtTime Orbit a
o forall a b. (a -> b) -> a -> b
$ Time a
t
  Classification
_ -> 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 :: forall a.
(Converge [a], Floating a, Real a) =>
Orbit a -> Angle a -> Maybe (Angle a)
eccentricAnomalyAtMeanAnomaly Orbit a
o Angle a
_M = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
                                       Classification
Elliptic -> Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
_E
                                       Classification
_ -> forall a. Maybe a
Nothing
  where (Qu '[] 'DefaultLCSU Integer
n, a
wrappedM) = forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second (forall (dim :: [Factor (*)]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si|rad|]) (Angle a
_M 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'` forall a. Floating a => PlaneAngle a
turn)
        e :: a
e = forall a. Orbit a -> Unitless a
eccentricity Orbit a
o 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 = forall a. Fractional a => a -> Angle a
rad forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (Real a, Fractional b) => a -> b
realToFrac forall a b. (a -> b) -> a -> b
$ a
wrappedM
        oFloat :: Orbit Float
oFloat = forall a b. (a -> b) -> Orbit a -> Orbit b
unsafeMapOrbit forall a b. (Real a, Fractional b) => a -> b
realToFrac Orbit a
o
        initialGuessFloat :: Angle Float
        Just Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
Angle Float
initialGuessFloat = Orbit Float -> Angle Float -> Maybe (Angle Float)
eccentricAnomalyAtMeanAnomalyFloat Orbit Float
oFloat Qu '[ 'F PlaneAngle One] 'DefaultLCSU Float
_MFloat
        initialGuess :: a
initialGuess = forall a b. (Real a, Fractional b) => a -> b
realToFrac forall b c a. (b -> c) -> (a -> b) -> a -> c
. (forall (dim :: [Factor (*)]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si|rad|]) forall a b. (a -> b) -> a -> b
$ Angle Float
initialGuessFloat
        err :: (Mode b, Floating b, Scalar b ~ a) => b -> b
        err :: forall b. (Mode b, Floating b, Scalar b ~ a) => b -> b
err b
_E = forall t. Mode t => Scalar t -> t
auto a
wrappedM forall a. Num a => a -> a -> a
- (b
_E forall a. Num a => a -> a -> a
- forall t. Mode t => Scalar t -> t
auto a
e forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin b
_E)
        wrappedE :: Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
wrappedE = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a. Fractional a => a -> Angle a
rad forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a.
(Converge a, Ord (Element a)) =>
(Element a -> Element a) -> a -> Maybe (Element a)
convergeErr (forall a. Id a -> a
runId forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => a -> a
abs forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall b. (Mode b, Floating b, Scalar b ~ a) => b -> b
err forall b c a. (b -> c) -> (a -> b) -> a -> c
.  forall a. a -> Id a
Id) forall a b. (a -> b) -> a -> b
$
                   forall a. Fractional a => (Tower a -> Tower a) -> a -> [a]
findZeroNoEq forall b. (Mode b, Floating b, Scalar b ~ a) => b -> b
err a
initialGuess
        _E :: Maybe (Qu '[ 'F PlaneAngle One] 'DefaultLCSU a)
_E = (forall (d1 :: [Factor (*)]) (d2 :: [Factor (*)]) n (l :: LCSU (*)).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| (forall a b (u :: [Factor (*)]) (l :: LCSU (*)).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit forall a. Num a => Integer -> a
fromInteger Qu '[] 'DefaultLCSU Integer
n forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| forall a. Floating a => PlaneAngle a
turn)) 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 Orbit Float
o Angle Float
_M = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit Float
o of
                                            Classification
Elliptic -> forall a. a -> Maybe a
Just Angle Float
_E
                                            Classification
_ -> forall a. Maybe a
Nothing
  where wrappedM :: Float
wrappedM = (Angle Float
_M forall a (u :: [Factor (*)]) (l :: LCSU (*)).
Real a =>
Qu u l a -> Qu u l a -> Qu u l a
`mod'` forall a. Floating a => PlaneAngle a
turn) forall (dim :: [Factor (*)]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si|rad|]
        e :: Float
e = forall a. Orbit a -> Unitless a
eccentricity Orbit Float
o forall (dim :: [Factor (*)]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si||]
        sinM :: Float
sinM = forall a. Floating a => a -> a
sin Float
wrappedM
        cosM :: Float
cosM = 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 forall a. Num a => a -> a -> a
+
                       Float
e forall a. Num a => a -> a -> a
* Float
sinM forall a. Num a => a -> a -> a
+
                       Float
e forall a. Num a => a -> a -> a
* Float
e forall a. Num a => a -> a -> a
* Float
sinM forall a. Num a => a -> a -> a
* Float
cosM forall a. Num a => a -> a -> a
+
                       Float
0.5 forall a. Num a => a -> a -> a
* Float
e forall a. Num a => a -> a -> a
* Float
e forall a. Num a => a -> a -> a
* Float
e forall a. Num a => a -> a -> a
* Float
sinM forall a. Num a => a -> a -> a
* (Float
3 forall a. Num a => a -> a -> a
* Float
cosM forall a. Num a => a -> a -> a
* Float
cosM forall a. Num a => a -> a -> a
- Float
1)
        _E :: Angle Float
        _E :: Angle Float
_E = forall a. Fractional a => a -> Angle a
rad forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. [a] -> a
last forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Int -> [a] -> [a]
take Int
5 forall a b. (a -> b) -> a -> b
$
             forall a. (Fractional a, Eq a) => (Tower a -> Tower a) -> a -> [a]
findZero (\Tower Float
_E -> forall t. Mode t => Scalar t -> t
auto Float
wrappedM forall a. Num a => a -> a -> a
- (Tower Float
_E forall a. Num a => a -> a -> a
- forall t. Mode t => Scalar t -> t
auto Float
e forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin 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 :: forall a.
(Floating a, Real a) =>
Orbit a -> Angle a -> Maybe (Angle a)
eccentricAnomalyAtTrueAnomaly Orbit a
o Angle a
ν = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
                                       Classification
Elliptic -> forall a. a -> Maybe a
Just Qu (Normalize '[ 'F PlaneAngle One]) 'DefaultLCSU a
_E
                                       Classification
_ -> forall a. Maybe a
Nothing
  where (Qu '[] 'DefaultLCSU Integer
n, Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
wrappedν) = Angle 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'` forall a. Floating a => PlaneAngle a
turn
        cosν :: a
cosν = forall a. Floating a => a -> a
cos (Angle 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 = forall a. Orbit a -> Unitless a
eccentricity Orbit a
o forall (dim :: [Factor (*)]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si||]
        wrappedE :: Angle a
wrappedE = forall a. Fractional a => a -> Angle a
rad forall a b. (a -> b) -> a -> b
$ forall a. Floating a => a -> a
acos ((a
e forall a. Num a => a -> a -> a
+ a
cosν) forall a. Fractional a => a -> a -> a
/ (a
1 forall a. Num a => a -> a -> a
+ a
e forall a. Num a => a -> a -> a
* a
cosν))
        -- wrappedE = rad $ atan2 (sqrt (1 - e*e) * sinν) (e + cosν)
        _E :: Qu (Normalize '[ 'F PlaneAngle One]) 'DefaultLCSU a
_E = if Qu '[ 'F PlaneAngle One] 'DefaultLCSU a
wrappedν forall a. Ord a => a -> a -> Bool
< forall a. Floating a => PlaneAngle a
halfTurn
               then (forall a b (u :: [Factor (*)]) (l :: LCSU (*)).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit forall a. Num a => Integer -> a
fromInteger Qu '[] 'DefaultLCSU Integer
n forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| forall a. Floating a => PlaneAngle a
turn) 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
wrappedE
               else (forall a b (u :: [Factor (*)]) (l :: LCSU (*)).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit forall a. Num a => Integer -> a
fromInteger (Qu '[] 'DefaultLCSU Integer
nforall a. Num a => a -> a -> a
+Qu '[] 'DefaultLCSU Integer
1) forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| forall a. Floating a => PlaneAngle a
turn) 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
wrappedE

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

hyperbolicAnomalyAtTime
  :: forall a
   . (Converge [a], RealFloat a)
  => Orbit a
  -> Time a
  -> Maybe (AngleH a)
hyperbolicAnomalyAtTime :: forall a.
(Converge [a], RealFloat a) =>
Orbit a -> Time a -> Maybe (AngleH a)
hyperbolicAnomalyAtTime Orbit a
o =
  forall a.
(Converge [a], RealFloat a) =>
Orbit a -> Angle a -> Maybe (AngleH a)
hyperbolicAnomalyAtMeanAnomaly Orbit a
o forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 :: forall a.
(Converge [a], RealFloat a) =>
Orbit a -> Angle a -> Maybe (AngleH a)
hyperbolicAnomalyAtMeanAnomaly Orbit a
o Angle a
_M = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Classification
Hyperbolic -> Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a)
_H
  Classification
_          -> forall a. Maybe a
Nothing
 where
  e :: a
e                       = forall a. Orbit a -> Unitless a
eccentricity Orbit a
o forall (dim :: [Factor (*)]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si||]
  _M' :: a
_M'                     = Angle a
_M forall (dim :: [Factor (*)]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si|rad|]
  _MDouble :: Double
_MDouble                = forall a b. (Real a, Fractional b) => a -> b
realToFrac a
_M'
  Just AngleH Double
initialGuessDouble = Orbit Double -> Angle Double -> Maybe (AngleH Double)
hyperbolicAnomalyAtMeanAnomalyDouble
    (forall a b. (a -> b) -> Orbit a -> Orbit b
unsafeMapOrbit forall a b. (Real a, Fractional b) => a -> b
realToFrac Orbit a
o)
    (forall a. Fractional a => a -> Angle a
rad Double
_MDouble)
  initialGuess :: a
initialGuess = forall a b. (Real a, Fractional b) => a -> b
realToFrac forall b c a. (b -> c) -> (a -> b) -> a -> c
. (forall (dim :: [Factor (*)]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# RadianHyperbolic
RadianHyperbolic) forall a b. (a -> b) -> a -> b
$ AngleH Double
initialGuessDouble
  err :: (Mode b, Floating b, Scalar b ~ a) => b -> b
  err :: forall b. (Mode b, Floating b, Scalar b ~ a) => b -> b
err b
_H = forall t. Mode t => Scalar t -> t
auto a
_M' forall a. Num a => a -> a -> a
- (forall t. Mode t => Scalar t -> t
auto a
e forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sinh b
_H forall a. Num a => a -> a -> a
- b
_H)
  _H :: Maybe (Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a)
_H = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a. Fractional a => a -> AngleH a
rdh forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a.
(Converge a, Ord (Element a)) =>
(Element a -> Element a) -> a -> Maybe (Element a)
convergeErr (forall a. Id a -> a
runId forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => a -> a
abs forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall b. (Mode b, Floating b, Scalar b ~ a) => b -> b
err forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. a -> Id a
Id) forall a b. (a -> b) -> a -> b
$ forall a. Fractional a => (Tower a -> Tower a) -> a -> [a]
findZeroNoEq
    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 Orbit Double
o Angle Double
_M = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit Double
o of
  Classification
Hyperbolic -> case Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
_H of
    -- If you hit this, a better initial guess would probably help
    Qu Double
x | forall a. RealFloat a => a -> Bool
isNaN Double
x -> forall a. HasCallStack => String -> a
error String
"NaN while trying to find hyperbolic anomaly"
    Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
_              -> forall a. a -> Maybe a
Just Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU Double
_H
  Classification
_ -> 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            = forall a. Orbit a -> Unitless a
eccentricity Orbit Double
o forall (dim :: [Factor (*)]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si||]
  _M' :: Double
_M'          = Angle Double
_M 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           = forall a. Fractional a => a -> AngleH a
rdh forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. [a] -> a
last forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Int -> [a] -> [a]
take Int
200 forall a b. (a -> b) -> a -> b
$ (forall s. AD s ForwardDouble -> AD s ForwardDouble)
-> Double -> [Double]
Newton.findZero
    (\AD s ForwardDouble
_H -> forall t. Mode t => Scalar t -> t
auto Double
_M' forall a. Num a => a -> a -> a
- (forall t. Mode t => Scalar t -> t
auto Double
e forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sinh AD s ForwardDouble
_H 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 :: forall a.
(Floating a, Ord a) =>
Orbit a -> Angle a -> Maybe (AngleH a)
hyperbolicAnomalyAtTrueAnomaly Orbit a
o Angle a
ν = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Classification
_ | Just Angle a
d <- forall a. (Floating a, Ord a) => Orbit a -> Maybe (Angle a)
hyperbolicDepartureAngle Orbit a
o, forall a (l :: LCSU (*)) (u :: [Factor (*)]).
Num a =>
Qu u l a -> Qu u l a
qAbs Angle a
ν forall (d1 :: [Factor (*)]) (d2 :: [Factor (*)]) n (l :: LCSU (*)).
(d1 @~ d2, Ord n) =>
Qu d1 l n -> Qu d2 l n -> Bool
|>| Angle a
d -> forall a. Maybe a
Nothing
  Classification
Hyperbolic -> forall a. a -> Maybe a
Just Qu '[ 'F PlaneAngleHyperbolic One] 'DefaultLCSU a
_H
  Classification
_          -> forall a. Maybe a
Nothing
 where
  e :: Unitless a
e     = forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
  coshH :: Qu '[] 'DefaultLCSU a
coshH = (forall a. Floating a => Angle a -> Unitless a
qCos Angle a
ν forall a. Num a => a -> a -> a
+ Unitless a
e) forall a. Fractional a => a -> a -> a
/ (Qu '[] 'DefaultLCSU a
1 forall a. Num a => a -> a -> a
+ Unitless a
e forall a. Num a => a -> a -> a
* forall a. Floating a => Angle a -> Unitless a
qCos Angle a
ν)
  sign :: a
sign  = forall a. Num a => a -> a
signum (Angle 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 forall n (b :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
n -> Qu b l n -> Qu b l n
*| forall a. Floating a => Unitless a -> AngleH a
qArcCosh Qu '[] 'DefaultLCSU 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 :: forall a.
(Converge [a], RealFloat a) =>
Orbit a -> Time a -> Angle a
trueAnomalyAtTime Orbit a
o Time a
t = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Classification
Elliptic   -> forall a.
(Converge [a], RealFloat a) =>
Orbit a -> Angle a -> Angle a
trueAnomalyAtMeanAnomaly Orbit a
o Angle a
_M
  Classification
Hyperbolic -> forall a.
(Converge [a], RealFloat a) =>
Orbit a -> Angle a -> Angle a
trueAnomalyAtMeanAnomaly Orbit a
o Angle a
_M
  Classification
Parabolic ->
    let _A :: Qu
  (Normalize
     (Normalize
        (Normalize ('[] @- '[])
         @@+ Reorder
               (Normalize
                  ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
                   @- Normalize
                        ('[]
                         @@+ Reorder
                               (Normalize
                                  (Normalize
                                     (Normalize ('[ 'F Length One] @- '[])
                                      @@+ Normalize ('[ 'F Length One] @- '[]))
                                   @@+ Reorder
                                         (Normalize ('[ 'F Length One] @- '[]))
                                         (Normalize
                                            (Normalize ('[ 'F Length One] @- '[])
                                             @@+ Normalize ('[ 'F Length One] @- '[])))))
                               '[]))
                @/ 'S One)
               (Normalize ('[] @- '[])))
      @+ '[ 'F Time One]))
  'DefaultLCSU
  a
_A = (Qu '[] 'DefaultLCSU a
3 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 '[] 'DefaultLCSU a
2) forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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 '[] 'DefaultLCSU a
2 forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu a l n -> Qu (Normalize (Normalize (a @+ a) @+ a)) l n
qCube (Distance a
l 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 '[] 'DefaultLCSU a
2))) forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Time a
t
        _B :: Qu
  (Normalize
     (Normalize
        (Normalize ('[] @- '[])
         @@+ Reorder
               (Normalize
                  ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
                   @- Normalize
                        ('[]
                         @@+ Reorder
                               (Normalize
                                  (Normalize
                                     (Normalize ('[ 'F Length One] @- '[])
                                      @@+ Normalize ('[ 'F Length One] @- '[]))
                                   @@+ Reorder
                                         (Normalize ('[ 'F Length One] @- '[]))
                                         (Normalize
                                            (Normalize ('[ 'F Length One] @- '[])
                                             @@+ Normalize ('[ 'F Length One] @- '[])))))
                               '[]))
                @/ 'S One)
               (Normalize ('[] @- '[])))
      @@+ '[ 'F Time One])
   @/ 'S ('S One))
  'DefaultLCSU
  a
_B = forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Floating n =>
Qu a l n -> Qu (a @/ 'S ('S One)) l n
qCubeRoot (Qu
  (Normalize
     (Normalize
        (Normalize ('[] @- '[])
         @@+ Reorder
               (Normalize
                  ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
                   @- Normalize
                        ('[]
                         @@+ Reorder
                               (Normalize
                                  (Normalize
                                     (Normalize ('[ 'F Length One] @- '[])
                                      @@+ Normalize ('[ 'F Length One] @- '[]))
                                   @@+ Reorder
                                         (Normalize ('[ 'F Length One] @- '[]))
                                         (Normalize
                                            (Normalize ('[ 'F Length One] @- '[])
                                             @@+ Normalize ('[ 'F Length One] @- '[])))))
                               '[]))
                @/ 'S One)
               (Normalize ('[] @- '[])))
      @+ '[ 'F Time One]))
  'DefaultLCSU
  a
_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
|+| forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu a l n -> Qu (Normalize (a @+ a)) l n
qSq Qu
  (Normalize
     (Normalize
        (Normalize ('[] @- '[])
         @@+ Reorder
               (Normalize
                  ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
                   @- Normalize
                        ('[]
                         @@+ Reorder
                               (Normalize
                                  (Normalize
                                     (Normalize ('[ 'F Length One] @- '[])
                                      @@+ Normalize ('[ 'F Length One] @- '[]))
                                   @@+ Reorder
                                         (Normalize ('[ 'F Length One] @- '[]))
                                         (Normalize
                                            (Normalize ('[ 'F Length One] @- '[])
                                             @@+ Normalize ('[ 'F Length One] @- '[])))))
                               '[]))
                @/ 'S One)
               (Normalize ('[] @- '[])))
      @+ '[ 'F Time One]))
  'DefaultLCSU
  a
_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
1))
    in  Qu '[] 'DefaultLCSU a
2 forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| forall a. Floating a => Unitless a -> Angle a
qArcTan (Qu
  (Normalize
     (Normalize
        (Normalize ('[] @- '[])
         @@+ Reorder
               (Normalize
                  ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
                   @- Normalize
                        ('[]
                         @@+ Reorder
                               (Normalize
                                  (Normalize
                                     (Normalize ('[ 'F Length One] @- '[])
                                      @@+ Normalize ('[ 'F Length One] @- '[]))
                                   @@+ Reorder
                                         (Normalize ('[ 'F Length One] @- '[]))
                                         (Normalize
                                            (Normalize ('[ 'F Length One] @- '[])
                                             @@+ Normalize ('[ 'F Length One] @- '[])))))
                               '[]))
                @/ 'S One)
               (Normalize ('[] @- '[])))
      @@+ '[ 'F Time One])
   @/ 'S ('S One))
  'DefaultLCSU
  a
_B forall a. Num a => a -> a -> a
- forall a. Fractional a => a -> a
recip Qu
  (Normalize
     (Normalize
        (Normalize ('[] @- '[])
         @@+ Reorder
               (Normalize
                  ('[ 'F Length ('S ('S One)), 'F Time ('P ('P 'Zero))]
                   @- Normalize
                        ('[]
                         @@+ Reorder
                               (Normalize
                                  (Normalize
                                     (Normalize ('[ 'F Length One] @- '[])
                                      @@+ Normalize ('[ 'F Length One] @- '[]))
                                   @@+ Reorder
                                         (Normalize ('[ 'F Length One] @- '[]))
                                         (Normalize
                                            (Normalize ('[ 'F Length One] @- '[])
                                             @@+ Normalize ('[ 'F Length One] @- '[])))))
                               '[]))
                @/ 'S One)
               (Normalize ('[] @- '[])))
      @@+ '[ 'F Time One])
   @/ 'S ('S One))
  'DefaultLCSU
  a
_B)
 where
  μ :: 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  = forall a. Num a => Orbit a -> Distance a
semiLatusRectum Orbit a
o
  _M :: Angle a
_M = 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 :: forall a.
(Converge [a], RealFloat a) =>
Orbit a -> Angle a -> Angle a
trueAnomalyAtMeanAnomaly Orbit a
o Angle a
_M = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Classification
Elliptic -> forall a. HasCallStack => Maybe a -> a
fromJust
    (forall a. RealFloat a => Orbit a -> Angle a -> Maybe (Angle a)
trueAnomalyAtEccentricAnomaly Orbit a
o forall (m :: * -> *) b c a.
Monad m =>
(b -> m c) -> (a -> m b) -> a -> m c
<=< forall a.
(Converge [a], Floating a, Real a) =>
Orbit a -> Angle a -> Maybe (Angle a)
eccentricAnomalyAtMeanAnomaly Orbit a
o forall a b. (a -> b) -> a -> b
$ Angle a
_M)
  Classification
Hyperbolic -> forall a. HasCallStack => Maybe a -> a
fromJust
    (forall a.
(Ord a, Floating a) =>
Orbit a -> AngleH a -> Maybe (Angle a)
trueAnomalyAtHyperbolicAnomaly Orbit a
o forall (m :: * -> *) b c a.
Monad m =>
(b -> m c) -> (a -> m b) -> a -> m c
<=< forall a.
(Converge [a], RealFloat a) =>
Orbit a -> Angle a -> Maybe (AngleH a)
hyperbolicAnomalyAtMeanAnomaly Orbit a
o forall a b. (a -> b) -> a -> b
$ Angle a
_M)
  Classification
_ -> forall a. HasCallStack => String -> a
error String
"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 :: forall a. RealFloat a => Orbit a -> Angle a -> Maybe (Angle a)
trueAnomalyAtEccentricAnomaly Orbit a
o Angle a
_E = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
                                       Classification
Elliptic -> forall a. a -> Maybe a
Just Qu (Normalize '[ 'F PlaneAngle One]) 'DefaultLCSU a
ν
                                       Classification
_        -> forall a. Maybe a
Nothing
  where (Qu '[] 'DefaultLCSU a
n, a
wrappedE) = forall (p :: * -> * -> *) a b c d.
Bifunctor p =>
(a -> b) -> (c -> d) -> p a c -> p b d
bimap (forall a b (u :: [Factor (*)]) (l :: LCSU (*)).
(a -> b) -> Qu u l a -> Qu u l b
unsafeMapUnit forall a. Num a => Integer -> a
fromInteger) (forall (dim :: [Factor (*)]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si|rad|]) forall a b. (a -> b) -> a -> b
$
                        Angle a
_E 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'` forall a. Floating a => PlaneAngle a
turn
        e :: a
e = forall a. Orbit a -> Unitless a
eccentricity Orbit a
o forall (dim :: [Factor (*)]) unit n.
(ValidDLU dim 'DefaultLCSU unit, Fractional n) =>
Qu dim 'DefaultLCSU n -> unit -> n
# [si||]
        wrappedν :: Angle a
wrappedν = forall a. Fractional a => a -> Angle a
rad forall a b. (a -> b) -> a -> b
$ a
2 forall a. Num a => a -> a -> a
* forall a. RealFloat a => a -> a -> a
atan2 (forall a. Floating a => a -> a
sqrt (a
1 forall a. Num a => a -> a -> a
+ a
e) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin (a
wrappedE forall a. Fractional a => a -> a -> a
/ a
2))
                                   (forall a. Floating a => a -> a
sqrt (a
1 forall a. Num a => a -> a -> a
- a
e) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos (a
wrappedE forall a. Fractional a => a -> a -> a
/ a
2))
        ν :: Qu (Normalize '[ 'F PlaneAngle One]) 'DefaultLCSU a
ν = forall a. Floating a => PlaneAngle a
turn 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 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
wrappedν

trueAnomalyAtHyperbolicAnomaly
  :: (Ord a, Floating a) => Orbit a -> AngleH a -> Maybe (Angle a)
trueAnomalyAtHyperbolicAnomaly :: forall a.
(Ord a, Floating a) =>
Orbit a -> AngleH a -> Maybe (Angle a)
trueAnomalyAtHyperbolicAnomaly Orbit a
o AngleH a
_H = case forall a. (Num a, Ord a) => Orbit a -> Classification
classify Orbit a
o of
  Classification
Hyperbolic -> forall a. a -> Maybe a
Just Qu (Normalize ('[] @+ '[ 'F PlaneAngle One])) 'DefaultLCSU a
ν
  Classification
_          -> forall a. Maybe a
Nothing
 where
  e :: Unitless a
e         = forall a. Orbit a -> Unitless a
eccentricity Orbit a
o
  tanνOver2 :: Qu '[] 'DefaultLCSU a
tanνOver2 = forall a. Floating a => a -> a
sqrt ((Unitless a
e forall a. Num a => a -> a -> a
+ Qu '[] 'DefaultLCSU a
1) forall a. Fractional a => a -> a -> a
/ (Unitless a
e forall a. Num a => a -> a -> a
- Qu '[] 'DefaultLCSU a
1)) forall a. Num a => a -> a -> a
* forall a. Floating a => AngleH a -> Unitless a
qTanh (AngleH a
_H 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 '[] 'DefaultLCSU a
2)
  ν :: Qu (Normalize ('[] @+ '[ 'F PlaneAngle One])) 'DefaultLCSU a
ν         = Qu '[] 'DefaultLCSU a
2 forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| forall a. Floating a => Unitless a -> Angle a
qArcTan Qu '[] 'DefaultLCSU 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 :: forall a. (Ord a, Floating a) => Orbit a -> Angle a -> Distance a
radiusAtTrueAnomaly Orbit a
o Angle a
trueAnomaly = case forall a. (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)
semiMajorAxis Orbit a
o of
  Just Distance a
_  -> Distance a
l 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 '[] 'DefaultLCSU a
1 forall (d1 :: [Factor (*)]) (d2 :: [Factor (*)]) n (l :: LCSU (*)).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| Unitless a
e forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| forall a. Floating a => Angle a -> Unitless a
qCos Angle a
ν)
  Maybe (Distance a)
Nothing -> (forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu a l n -> Qu (Normalize (a @+ a)) l n
qSq Quantity ((Meter :^ Succ (Succ 'Zero)) :* (Second :^ Pred 'Zero)) a
h forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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
1 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 '[] 'DefaultLCSU a
1 forall (d1 :: [Factor (*)]) (d2 :: [Factor (*)]) n (l :: LCSU (*)).
(d1 @~ d2, Num n) =>
Qu d1 l n -> Qu d2 l n -> Qu d1 l n
|+| forall a. Floating a => Angle a -> Unitless a
qCos Angle a
ν))
 where
  h :: Quantity ((Meter :^ Succ (Succ 'Zero)) :* (Second :^ Pred 'Zero)) a
h = forall a.
Floating a =>
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ 'Zero)) :* (Second :^ Pred 'Zero)) a
specificAngularMomentum Orbit a
o
  e :: Unitless a
e = 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
μ = forall a.
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ (Succ 'Zero)))
      :* (Second :^ Pred (Pred 'Zero)))
     a
primaryGravitationalParameter Orbit a
o
  l :: Distance a
l = 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 :: forall a. (Ord a, Floating a) => Orbit a -> Angle a -> Speed a
speedAtTrueAnomaly Orbit a
o Angle a
trueAnomaly = case forall a. (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)
semiMajorAxis Orbit a
o of
  Maybe (Distance a)
Nothing -> forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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
2 forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Distance a
r)
  Just Distance a
a  -> forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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
2 forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Distance a
r 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
1 forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Fractional n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @- b)) l n
|/| Distance a
a))
 where
  ν :: Angle a
ν = Angle a
trueAnomaly
  μ :: 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 = 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 :: forall a.
Floating a =>
Orbit a
-> Quantity
     ((Meter :^ Succ (Succ 'Zero)) :* (Second :^ Pred 'Zero)) a
specificAngularMomentum Orbit a
o = forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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
|*| Distance a
l)
  where
    μ :: 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 = 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 :: forall a.
(Ord a, Floating a) =>
Orbit a -> Quantity (Joule :/ (Kilo :@ Gram)) a
specificOrbitalEnergy Orbit a
o = case forall a. (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)
semiMajorAxis Orbit a
o of
  Just Distance a
a  -> forall n (d :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu d l n -> Qu d l n
qNegate (Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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 '[] 'DefaultLCSU a
2 forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Distance a
a))
  Maybe (Distance a)
Nothing -> 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
μ = 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 :: forall a.
(Ord a, Floating a) =>
Orbit a -> Angle a -> Quantity (Joule :/ (Kilo :@ Gram)) a
specificPotentialEnergyAtTrueAnomaly Orbit a
o Angle a
ν = forall n (d :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu d l n -> Qu d l n
qNegate (Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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
|/| Distance a
r)
 where
  r :: Distance a
r = 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
μ = 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 :: forall a.
(Ord a, Floating a) =>
Orbit a -> Angle a -> Quantity (Joule :/ (Kilo :@ Gram)) a
specificKineticEnergyAtTrueAnomaly Orbit a
o Angle a
ν = forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Num n =>
Qu a l n -> Qu (Normalize (a @+ a)) l n
qSq (forall a. (Ord a, Floating a) => Orbit a -> Angle a -> Speed a
speedAtTrueAnomaly Orbit a
o Angle 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 '[] 'DefaultLCSU a
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 :: forall a.
Floating a =>
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
μ Distance a
r = forall n (a :: [Factor (*)]) (l :: LCSU (*)).
Floating n =>
Qu a l n -> Qu (a @/ 'S One) l n
qSqrt (Qu '[] 'DefaultLCSU a
2 forall n (a :: [Factor (*)]) (l :: LCSU (*)) (b :: [Factor (*)]).
Num n =>
Qu a l n -> Qu b l n -> Qu (Normalize (a @+ b)) l n
|*| Quantity
  ((Meter :^ Succ (Succ (Succ 'Zero)))
   :* (Second :^ Pred (Pred 'Zero)))
  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
|/| Distance a
r)