orbits 0.3 → 0.4

raw patch · 17 files changed

+1811/−147 lines, 17 filesdep +lensdep +linearPVP ok

version bump matches the API change (PVP)

Dependencies added: lens, linear

API changes (from Hackage documentation)

Files

• changelog.md +6/−0
• default.nix +5/−4
• orbits.cabal +21/−15
• readme.md +9/−2
• src/Data/Metrology/Extra.hs +112/−5
• src/Physics/Orbit.hs +305/−51
• src/Physics/Orbit/Metrology.hs +27/−0
• src/Physics/Orbit/Sol.hs +86/−0
• src/Physics/Orbit/StateVectors.hs +559/−0
• test/Data/CReal/QuickCheck.hs +12/−3
• test/Data/Metrology/Extra.hs +140/−0
• test/Data/Metrology/QuickCheck.hs +2/−1
• test/Linear/QuickCheck.hs +18/−0
• test/Physics/Orbit/QuickCheck.hs +83/−6
• test/Test.hs +166/−59
• test/Test/QuickCheck/Extra.hs +59/−1
• test/Test/StateVectors.hs +201/−0

@@ -2,5 +2,11 @@  ## WIP +## 0.4 - 2020-09-29+  - Hyperbolic anomaly support+  - Several new utility functions+  - Data for Sol+  - Conversions to and from state vectors+ ## 0.3   - Switch to `units` from `uom-plugin`

@@ -1,4 +1,6 @@-{ pkgs ? import <nixpkgs> { }, compiler ? null, hoogle ? true }:+{ nixpkgsSrc ? builtins.fetchTarball+  "https://github.com/NixOS/nixpkgs/archive/1179840f9a88b8a548f4b11d1a03aa25a790c379.tar.gz"+, pkgs ? import nixpkgsSrc { }, compiler ? null, hoogle ? true }:  let   src = pkgs.nix-gitignore.gitignoreSource [ ] ./.;@@ -15,9 +17,8 @@       overrides = self: super:         {           exact-real = markUnbroken (dontCheck (doJailbreak super.exact-real));-          units-defs = self.callCabal2nix "" (builtins.fetchTarball-            "https://hackage.haskell.org/package/units-defs-2.2/units-defs-2.2.tar.gz")-            { };+          # reanimate = self.callCabal2nix "" ../../src/reanimate {};+          # reanimate-svg = self.callCabal2nix "" ../../src/reanimate-svg {};         } // pkgs.lib.optionalAttrs hoogle {           ghc = super.ghc // { withPackages = super.ghc.withHoogle; };           ghcWithPackages = self.ghc.withPackages;

@@ -1,17 +1,17 @@ cabal-version: 1.24 --- This file has been generated from package.yaml by hpack version 0.33.0.+-- This file has been generated from package.yaml by hpack version 0.33.1. -- -- see: https://github.com/sol/hpack ----- hash: 19ff836740d3a6d31f4b16ef06879a924c919fabab962046e35730122d362dff+-- hash: b0d20db4cd07be0ae62542dc01fb9fddb137451e34eb3a68439405d4619e86d9  name:           orbits-version:        0.3+version:        0.4 synopsis:       Types and functions for Kepler orbits. category:       Physics-homepage:       https://github.com/expipiplus1/orbit#readme-bug-reports:    https://github.com/expipiplus1/orbit/issues+homepage:       https://github.com/expipiplus1/orbits#readme+bug-reports:    https://github.com/expipiplus1/orbits/issues author:         Joe Hermaszewski maintainer:     Joe Hermaszewski <keep.it.real@monoid.al> copyright:      2020 Joe Hermaszewski@@ -26,7 +26,7 @@  source-repository head   type: git-  location: https://github.com/expipiplus1/orbit+  location: https://github.com/expipiplus1/orbits  custom-setup   setup-depends:@@ -38,16 +38,21 @@   exposed-modules:       Data.Constants.Mechanics.Extra       Physics.Orbit+      Physics.Orbit.Metrology+      Physics.Orbit.Sol+      Physics.Orbit.StateVectors   other-modules:       Data.Metrology.Extra   hs-source-dirs:       src-  default-extensions: DataKinds GeneralizedNewtypeDeriving QuasiQuotes ScopedTypeVariables TypeOperators+  default-extensions: DataKinds FlexibleContexts GeneralizedNewtypeDeriving LambdaCase QuasiQuotes ScopedTypeVariables TemplateHaskell TypeApplications TypeFamilies TypeOperators ViewPatterns   ghc-options: -Wall -O2   build-depends:       ad >=4.3.2     , base >=4.8 && <5     , exact-real >=0.12+    , lens+    , linear     , units     , units-defs >=2.2   default-language: Haskell2010@@ -59,8 +64,8 @@          hs-source-dirs:       test/doctest-  default-extensions: DataKinds GeneralizedNewtypeDeriving QuasiQuotes ScopedTypeVariables TypeOperators-  ghc-options: -Wall+  default-extensions: DataKinds FlexibleContexts GeneralizedNewtypeDeriving LambdaCase QuasiQuotes ScopedTypeVariables TemplateHaskell TypeApplications TypeFamilies TypeOperators ViewPatterns+  ghc-options: -Wall -O2   build-depends:       base     , doctest@@ -71,25 +76,26 @@   main-is: Test.hs   other-modules:       Data.CReal.QuickCheck+      Data.Metrology.Extra       Data.Metrology.QuickCheck+      Linear.QuickCheck       Physics.Orbit.QuickCheck       Test.QuickCheck.Extra+      Test.StateVectors       WrappedAngle-      Data.Constants.Mechanics.Extra-      Data.Metrology.Extra-      Physics.Orbit       Paths_orbits   hs-source-dirs:       test-      src-  default-extensions: DataKinds GeneralizedNewtypeDeriving QuasiQuotes ScopedTypeVariables TypeOperators-  ghc-options: -Wall -threaded+  default-extensions: DataKinds FlexibleContexts GeneralizedNewtypeDeriving LambdaCase QuasiQuotes ScopedTypeVariables TemplateHaskell TypeApplications TypeFamilies TypeOperators ViewPatterns+  ghc-options: -Wall -O2 -threaded   build-depends:       QuickCheck     , ad     , base     , checkers     , exact-real+    , lens+    , linear     , orbits     , random     , tagged

@@ -1,5 +1,4 @@-orbit-=====+# orbits  *For my uncle Zbys who watched the planets and stars.* @@ -27,7 +26,15 @@ | t      | Time since periapse              |                                | | M      | Mean anomaly                     |                                | | E      | Eccentric anomaly                | Only for elliptic orbits       |+| H      | Hyperbolic anomaly               | Only for hyperbolic orbits     | | ν      | True anomaly                     |                                |+| h      | Specific angular momentum        |                                |+| ε      | Specific orbital energy          |                                |+| εp     | Specific potential energy        |                                |+| εk     | Specific kinetic energy          |                                |+| v      | Orbital speed or velocity        |                                |+| r      | The radius to the orbiting body  |                                |+| φ      | Flight path angle                |                                |   Note that in the Haskell source uppercase symbols such as Ω and M are written

@@ -1,10 +1,10 @@-module Data.Metrology.Extra-  ( mod'-  , div'-  , divMod'-  ) where+{-# language QuasiQuotes #-} +module Data.Metrology.Extra where++import           Control.Applicative import           Data.Coerce                    ( coerce )+import           Data.Constants.Mechanics.Extra ( ) import qualified Data.Fixed                    as F                                                 ( div'                                                 , divMod'@@ -12,6 +12,11 @@                                                 ) import           Data.Metrology import           Data.Metrology.Unsafe          ( Qu(..) )+import           Data.Units.SI.Parser+import           Linear.Metric+import           Linear.V3+import           Linear.Vector+import           Physics.Orbit.Metrology  mod' :: forall a u l . Real a => Qu u l a -> Qu u l a -> Qu u l a mod' = coerce (F.mod' :: a -> a -> a)@@ -31,3 +36,105 @@   -> Qu u l a   -> (Qu '[] l b, Qu u l a) divMod' = coerce (F.divMod' :: a -> a -> (b, a))++rad :: Fractional a => a -> Angle a+rad = (% [si|rad|])++rdh :: Fractional a => a -> AngleH a+rdh = (% RadianHyperbolic)++qCos :: Floating a => Angle a -> Unitless a+qCos θ = quantity $ cos (θ # [si|rad|])++qSin :: Floating a => Angle a -> Unitless a+qSin θ = quantity $ sin (θ # [si|rad|])++qTan :: Floating a => Angle a -> Unitless a+qTan θ = quantity $ tan (θ # [si|rad|])++qArcTan :: Floating a => Unitless a -> Angle a+qArcTan = rad . atan . (# [si||])++qArcTan2 :: RealFloat a => Unitless a -> Unitless a -> Angle a+qArcTan2 x y = rad (atan2 (x # [si||]) (y # [si||]))++qArcCos :: Floating a => Unitless a -> Angle a+qArcCos = rad . acos . (# [si||])++qRecip+  :: forall u l a . Fractional a => Qu u l a -> Qu (Normalize ('[] @- u)) l a+qRecip = coerce (recip @a)++qTanh :: Floating a => AngleH a -> Unitless a+qTanh = quantity . tanh . (# RadianHyperbolic)++qSinh :: Floating a => AngleH a -> Unitless a+qSinh = quantity . sinh . (# RadianHyperbolic)++qCosh :: Floating a => AngleH a -> Unitless a+qCosh = quantity . cosh . (# RadianHyperbolic)++qArcCosh :: Floating a => Unitless a -> AngleH a+qArcCosh = rdh . acosh . (# [si||])++qAbs :: forall a l u . Num a => Qu u l a -> Qu u l a+qAbs = coerce (abs @a)++qCross+  :: Num n+  => V3 (Qu a l n)+  -> V3 (Qu b l n)+  -> V3 (Qu (Normalize (a @@+ Reorder b a)) l n)+qCross (V3 a b c) (V3 d e f) =+  V3 (b |*| f |-| c |*| e) (c |*| d |-| a |*| f) (a |*| e |-| b |*| d)++qNorm :: forall u l a . Floating a => V3 (Qu u l a) -> Qu u l a+qNorm = coerce (norm @V3 @a)++-- qNormalize+--   :: forall u l a . (Floating a, Epsilon a) => V3 (Qu u l a) -> V3 (Qu '[] l a)+-- qNormalize = coerce (normalize @a @V3)+qNormalize+  :: Floating n+  => V3 (Qu b l n)+  -> V3+       ( Qu+           ( Normalize+               (Normalize ('[] @- b) @@+ Reorder b (Normalize ('[] @- b)))+           )+           l+           n+       )+qNormalize x = (qRecip (qNorm x) |*|) <$> x++qDot+  :: forall u v l a. Num a+  => V3 (Qu u l a)+  -> V3 (Qu v l a)+  -> Qu (Normalize (u @@+ Reorder v u)) l a+qDot = coerce (dot @V3 @a)++qQuadrance+  :: forall u l a+   . Num a+  => V3 (Qu u l a)+  -> Qu (Normalize (u @@+ Reorder u u)) l a+qQuadrance = coerce (quadrance @V3 @a)++(|^/|) :: (Functor f, Fractional n) =>+            f (Qu b l n)+            -> Qu u l n+            -> f (Qu+                    (Normalize+                       (Normalize ('[] @- u) @@+ Reorder b (Normalize ('[] @- u))))+                    l+                    n)+x |^/| y = (qRecip y |*|) <$> x++(|^-^|)+  :: forall f u l a+   . (Additive f, Applicative f, Num a)+  => f (Qu u l a)+  -> f (Qu u l a)+  -> f (Qu u l a)+(|^-^|) = liftA2 (|-|)

@@ -1,9 +1,4 @@-{-# LANGUAGE FlexibleContexts      #-}-{-# LANGUAGE FlexibleInstances     #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE QuasiQuotes           #-}-{-# LANGUAGE TypeFamilies          #-}-{-# OPTIONS_GHC -fno-warn-orphans #-}+{-# language QuasiQuotes #-}  -- | Types and functions for dealing with Kepler orbits. module Physics.Orbit@@ -17,6 +12,7 @@     -- ** Utilities   , isValid   , classify+  , normalizeOrbit     -- ** Orbital elements   , apoapsis   , meanMotion@@ -33,11 +29,13 @@     -- *** To time since periapse   , timeAtMeanAnomaly   , timeAtEccentricAnomaly+  , timeAtHyperbolicAnomaly   , timeAtTrueAnomaly      -- *** To mean anomaly   , meanAnomalyAtTime   , meanAnomalyAtEccentricAnomaly+  , meanAnomalyAtHyperbolicAnomaly   , meanAnomalyAtTrueAnomaly      -- *** To eccentric anomaly@@ -46,17 +44,39 @@   , 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'@@ -71,10 +91,13 @@                                                 , 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          ( UnsafeQu(..) )+import           Data.Metrology.Unsafe          ( Qu(..)+                                                , UnsafeQu(..)+                                                ) import           Data.Units.SI.Parser import           Numeric.AD                     ( Mode                                                 , Scalar@@ -84,25 +107,13 @@                                                 , findZeroNoEq                                                 ) import           Numeric.AD.Internal.Identity   ( Id(..) )+import qualified Numeric.AD.Newton.Double      as Newton+import           Physics.Orbit.Metrology  -------------------------------------------------------------------------------- -- Types -------------------------------------------------------------------------------- -type Quantity u = MkQu_ULN u 'DefaultLCSU--- | A measure in seconds.-type Time     = Quantity [si|s|]--- | A measure in meters.-type Distance = Quantity [si| m |]--- | A measure in meters per second.-type Speed    = Quantity [si| m s^-1 |]--- | A measure in kilograms.-type Mass     = Quantity [si| kg |]--- | A measure in radians.-type Angle    = Quantity [si| rad |]--- | A unitless measure.-type Unitless = Quantity [si||]- -- | Data type defining an orbit parameterized by the type used to -- represent values data Orbit a = Orbit { -- | The orbit's eccentricity, e.@@ -179,13 +190,14 @@                                         -- angle of the periapsis relative to                                         -- the reference direction in the                                         -- orbital plane.-                                        argumentOfPeriapsis :: !(Angle a) }+                                        argumentOfPeriapsis :: !(Angle a)+                                      }                             -- | The orbit has an eccentricity of 0 so the                             -- 'argumentOfPeriapsis' is indeterminate.                           | Circular   deriving (Show, Eq) --- | What for the orbit's geometry takes. This is dependant only on the+-- | What form the orbit's geometry takes. This is dependant only on the -- 'eccentricity', e >= 0, of the orbit. data Classification = -- | 0 <= e < 1                       --@@ -224,9 +236,9 @@ -- Functions -------------------------------------------------------------------------------- --- | Return true is the orbit is valid and false if it is invalid. The behavior--- of all the other functions in this module is undefined when given an invalid--- orbit.+-- | 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 o = e >= 0 &&             ((e == 0) `iff` (periapsisSpecifier o == Circular)) &&@@ -238,16 +250,43 @@     q = periapsis o     μ = primaryGravitationalParameter o --- | 'classify' is a funciton which returns the orbit's class.+-- | What shape is the orbit classify :: (Num a, Ord a) => Orbit a -> Classification-classify o-  | e < 1 = Elliptic-  | e == 1 = Parabolic-  | e > 1 = Hyperbolic-  | otherwise = error "classify"-  where-    e = eccentricity o+classify o | e < 1     = Elliptic+           | e == 1    = Parabolic+           | e > 1     = Hyperbolic+           | otherwise = error "classify: NaN eccentricity"+  where e = eccentricity o +-- | Return an equivalent orbit such that+--+-- - i ∈ [0..π)+-- - Ω ∈ [0..2π)+-- - ω ∈ [0..2π)+-- - inclinationSpecifier == NonInclined if i = 0+-- - periapsisSpecifier == Circular if e == 0 and ω == 0+normalizeOrbit :: (Floating a, Real a) => Orbit a -> Orbit a+normalizeOrbit (Orbit e q inc per μ) = Orbit e q inc' per' μ+ where+  -- Were we actually given a descending node and have to flip things+  (inc', flipped) = case inc of+    NonInclined              -> (NonInclined, False)+    Inclined _ i | i == zero -> (NonInclined, False)+    Inclined _Ω i ->+      let iR  = i `mod'` turn+          i'  = if flipped then turn |-| iR else iR+          _Ω' = (if flipped then _Ω |+| halfTurn else _Ω) `mod'` turn+      in  (Inclined _Ω' i', iR >= halfTurn)++  per' = case per of+    Circular | flipped   -> Eccentric halfTurn+             | otherwise -> Circular+    Eccentric ω+      | ω == zero, e == 0, not flipped+      -> Circular+      | otherwise+      -> Eccentric $ (if flipped then ω |+| halfTurn else ω) `mod'` 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.@@ -262,7 +301,7 @@  -- | 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.+-- case of a parabolic orbit 'semiMinorAxis' returns @0m@. semiMinorAxis :: (Floating a, Ord a) => Orbit a -> Distance a semiMinorAxis o =   case classify o of@@ -358,6 +397,14 @@ hyperbolicApproachAngle :: (Floating a, Ord a) => Orbit a -> Maybe (Angle a) hyperbolicApproachAngle = fmap qNegate . 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.@@ -377,11 +424,37 @@ timeAtEccentricAnomaly :: (Floating a, Ord a) => Orbit a -> Angle a -> Maybe (Time a) timeAtEccentricAnomaly o = fmap (timeAtMeanAnomaly o) . meanAnomalyAtEccentricAnomaly 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 o =+  fmap (timeAtMeanAnomaly o) . meanAnomalyAtHyperbolicAnomaly o+ -- | Calculate the time since periapse given the true anomaly, ν, of an -- orbiting body.-timeAtTrueAnomaly :: (Real a, Floating a) => Orbit a -> Angle a -> Maybe (Time a)-timeAtTrueAnomaly o = fmap (timeAtMeanAnomaly o) . meanAnomalyAtTrueAnomaly o+--+-- 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 o ν = case classify o of+  _ | Just d <- hyperbolicDepartureAngle o, qAbs ν |>| d -> Nothing+  Parabolic ->+    let _D = qTan (ν |/| 2)+        t  = 0.5 |*| qSqrt (qCube l |/| μ) |*| (_D |+| (qCube _D |/| 3))+    in  Just t+  _ -> fmap (timeAtMeanAnomaly o) . meanAnomalyAtTrueAnomaly o $ ν+ where+  μ = primaryGravitationalParameter o+  l = semiLatusRectum 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.@@ -409,20 +482,40 @@         untypedE = delRad _E         _M = addRad (untypedE |-| e |*| sin 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 o _H = case classify o of+  Hyperbolic -> Just _M+  _          -> Nothing+ where+  e  = eccentricity o+  _M = addRad $ e * qSinh _H - quantity (_H # 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π ----- Currently only implemented for elliptic orbits.+-- 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 o = case classify o of+meanAnomalyAtTrueAnomaly o ν = case classify o of+  Parabolic -> Nothing   Elliptic -> meanAnomalyAtEccentricAnomaly o <=<-              eccentricAnomalyAtTrueAnomaly o-  _ -> error "TODO: meanAnomalyAtTrueAnomaly"+              eccentricAnomalyAtTrueAnomaly o $ ν+  Hyperbolic -> meanAnomalyAtHyperbolicAnomaly o <=<+                hyperbolicAnomalyAtTrueAnomaly o $ ν +---------+-- To eccentric+---------+ -- | Calculate the eccentric anomaly, E, of an elliptic orbit at time t. -- -- 'eccentricAnomalyAtTime' returns Nothing when given a parabolic or@@ -508,17 +601,108 @@                then (unsafeMapUnit fromInteger n |*| turn) |+| wrappedE                else (unsafeMapUnit fromInteger (n+1) |*| turn) |-| wrappedE +---------+-- To hyperbolic+---------++hyperbolicAnomalyAtTime+  :: forall a+   . (Converge [a], RealFloat a)+  => Orbit a+  -> Time a+  -> Maybe (AngleH a)+hyperbolicAnomalyAtTime o =+  hyperbolicAnomalyAtMeanAnomaly o . meanAnomalyAtTime o++hyperbolicAnomalyAtMeanAnomaly+  :: forall a+   . (Converge [a], RealFloat a)+  => Orbit a+  -> Angle a+  -> Maybe (AngleH a)+hyperbolicAnomalyAtMeanAnomaly o _M = case classify o of+  Hyperbolic -> _H+  _          -> Nothing+ where+  e                       = eccentricity o # [si||]+  _M'                     = _M # [si|rad|]+  _MDouble                = realToFrac _M'+  Just initialGuessDouble = hyperbolicAnomalyAtMeanAnomalyDouble+    (unsafeMapOrbit realToFrac o)+    (rad _MDouble)+  initialGuess = realToFrac . (# RadianHyperbolic) $ initialGuessDouble+  err :: (Mode b, Floating b, Scalar b ~ a) => b -> b+  err _H = auto _M' - (auto e * sinh _H - _H)+  _H = fmap rdh . convergeErr (runId . abs . err . Id) $ findZeroNoEq+    err+    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 o _M = case classify o of+  Hyperbolic -> case _H of+    -- If you hit this, a better initial guess would probably help+    Qu x | isNaN x -> error "NaN while trying to find hyperbolic anomaly"+    _              -> Just _H+  _ -> Nothing+ where+  -- Perhaps use something like https://www.researchgate.net/publication/226007277_A_Method_Solving_Kepler%27s_Equation_for_Hyperbolic_Case+  e            = eccentricity o # [si||]+  _M'          = _M # [si|rad|]+  -- TODO: A better guess here+  initialGuess = _M'+  _H           = rdh . last . take 200 $ Newton.findZero+    (\_H -> auto _M' - (auto e * sinh _H - _H))+    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 o ν = case classify o of+  _ | Just d <- hyperbolicDepartureAngle o, qAbs ν |>| d -> Nothing+  Hyperbolic -> Just _H+  _          -> Nothing+ where+  e     = eccentricity o+  coshH = (qCos ν + e) / (1 + e * qCos ν)+  sign  = signum (ν # [si|rad|])+  _H    = sign *| qArcCosh coshH++---------+-- To true+---------+ -- | Calculate the true anomaly, ν, of a body at time since periapse, t.-trueAnomalyAtTime :: (Converge [a], RealFloat a)-                  => Orbit a -> Time a -> Maybe (Angle a)-trueAnomalyAtTime o = trueAnomalyAtMeanAnomaly o . meanAnomalyAtTime o+trueAnomalyAtTime+  :: forall a . (Converge [a], RealFloat a) => Orbit a -> Time a -> Angle a+trueAnomalyAtTime o t = case classify o of+  Elliptic   -> trueAnomalyAtMeanAnomaly o _M+  Hyperbolic -> trueAnomalyAtMeanAnomaly o _M+  Parabolic ->+    let _A = (3 |/| 2) |*| qSqrt (μ |/| (2 |*| qCube (l |/| 2))) |*| t+        _B = qCubeRoot (_A |+| qSqrt (qSq _A |+| 1))+    in  2 |*| qArcTan (_B - recip _B)+ where+  μ  = primaryGravitationalParameter o+  l  = semiLatusRectum o+  _M = meanAnomalyAtTime o 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 -> Maybe (Angle a)-trueAnomalyAtMeanAnomaly o = trueAnomalyAtEccentricAnomaly o <=<-                             eccentricAnomalyAtMeanAnomaly o+trueAnomalyAtMeanAnomaly+  :: (Converge [a], RealFloat a) => Orbit a -> Angle a -> Angle a+trueAnomalyAtMeanAnomaly o _M = case classify o of+  Elliptic -> fromJust+    (trueAnomalyAtEccentricAnomaly o <=< eccentricAnomalyAtMeanAnomaly o $ _M)+  Hyperbolic -> fromJust+    (trueAnomalyAtHyperbolicAnomaly o <=< hyperbolicAnomalyAtMeanAnomaly o $ _M)+  _ -> error "trueAnomalyAtMeanAnomaly is not defined for Parabolic orbits"  -- | Calculate the true anomaly, ν, of an orbiting body when it has the given -- eccentric anomaly, _E.@@ -536,12 +720,82 @@                         _E `divMod'` turn         e = eccentricity o # [si||]         wrappedν = rad $ 2 * atan2 (sqrt (1 + e) * sin (wrappedE / 2))-                                        (sqrt (1 - e) * cos (wrappedE / 2))+                                   (sqrt (1 - e) * cos (wrappedE / 2))         ν = turn |*| n |+| wrappedν +trueAnomalyAtHyperbolicAnomaly+  :: (Ord a, Floating a) => Orbit a -> AngleH a -> Maybe (Angle a)+trueAnomalyAtHyperbolicAnomaly o _H = case classify o of+  Hyperbolic -> Just ν+  _          -> Nothing+ where+  e         = eccentricity o+  tanνOver2 = sqrt ((e + 1) / (e - 1)) * qTanh (_H |/| 2)+  ν         = 2 |*| qArcTan 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 o trueAnomaly = case semiMajorAxis o of+  Just _  -> l |/| (1 |+| e |*| qCos ν)+  Nothing -> (qSq h |/| μ) |*| (1 |/| (1 |+| qCos ν))+ where+  h = specificAngularMomentum o+  e = eccentricity o+  ν = trueAnomaly+  μ = primaryGravitationalParameter o+  l = semiLatusRectum 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 o trueAnomaly = case semiMajorAxis o of+  Nothing -> qSqrt (μ |*| 2 |/| r)+  Just a  -> qSqrt (μ |*| (2 |/| r |-| 1 |/| a))+ where+  ν = trueAnomaly+  μ = primaryGravitationalParameter o+  r = radiusAtTrueAnomaly o ν++-- | Specific angular momentum, h, is the angular momentum per unit mass+specificAngularMomentum :: Floating a => Orbit a -> Quantity [si|m^2 s^-1|] a+specificAngularMomentum o = qSqrt (μ |*| l)+  where+    μ = primaryGravitationalParameter o+    l = semiLatusRectum o++-- | Specific orbital energy, ε, is the orbital energy per unit mass+specificOrbitalEnergy+  :: (Ord a, Floating a) => Orbit a -> Quantity [si|J / kg|] a+specificOrbitalEnergy o = case semiMajorAxis o of+  Just a  -> qNegate (μ |/| (2 |*| a))+  Nothing -> zero+  where μ = primaryGravitationalParameter 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 o ν = qNegate (μ |/| r)+ where+  r = radiusAtTrueAnomaly o ν+  μ = primaryGravitationalParameter 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 o ν = qSq (speedAtTrueAnomaly o ν) |/| 2++---------------------------------------------------------------- -- Utils ---------------------------------------------------------------- -rad :: Fractional a => a -> Angle a-rad = (% [si|rad|])+-- | 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 μ r = qSqrt (2 |*| μ |/| r)

@@ -0,0 +1,27 @@+{-# language QuasiQuotes #-}++module Physics.Orbit.Metrology where++import           Data.Metrology+import           Data.Metrology.TH+import           Data.Units.SI.Parser++declareDimension "PlaneAngleHyperbolic"+declareCanonicalUnit "RadianHyperbolic" [t| PlaneAngleHyperbolic |] (Just "rdh")+type instance DefaultUnitOfDim PlaneAngleHyperbolic = RadianHyperbolic++type Quantity u = MkQu_ULN u 'DefaultLCSU+-- | A measure in seconds.+type Time     = Quantity [si|s|]+-- | A measure in meters.+type Distance = Quantity [si| m |]+-- | A measure in meters per second.+type Speed    = Quantity [si| m s^-1 |]+-- | A measure in kilograms.+type Mass     = Quantity [si| kg |]+-- | A measure in radians.+type Angle    = Quantity [si| rad |]+-- | A measure in radians (hyperbolic)+type AngleH   = Quantity RadianHyperbolic+-- | A unitless measure.+type Unitless = Quantity [si||]

@@ -0,0 +1,86 @@+module Physics.Orbit.Sol+  where++import           Data.Constants.Mechanics+import           Data.Metrology+import           Data.Units.Astronomical+import           Data.Units.SI.Parser+import           Physics.Orbit++solMass :: Fractional a => Mass a+solMass = 1988500e24 % [si|kg|]++solGraviationalParameter :: Fractional a => Quantity [si| m^3 s^-2 |] a+solGraviationalParameter = solMass |*| gravity_G++venusOrbit :: Fractional a => Orbit a+venusOrbit = Orbit+  { eccentricity                  = 0.006772+  , periapsis                     = 0.718440 % AstronomicalUnit+  , inclinationSpecifier          = Inclined+                                      { longitudeOfAscendingNode = 76.680 % [si|deg|]+                                      , inclination = 2.19 % [si|deg|]+                                      }+  , periapsisSpecifier = Eccentric { argumentOfPeriapsis = 54.884 % [si|deg|] }+  , primaryGravitationalParameter = solGraviationalParameter+  }++earthOrbit :: Fractional a => Orbit a+earthOrbit = Orbit+  { eccentricity                  = 0.01671123+  , periapsis                     = 0.9832899 % AstronomicalUnit+  , inclinationSpecifier          = Inclined+                                      { longitudeOfAscendingNode = 348.73936 % [si|deg|]+                                      , inclination = 1.578690 % [si|deg|]+                                      }+  , periapsisSpecifier = Eccentric { argumentOfPeriapsis = 114.20783 % [si|deg|]+                                   }+  , primaryGravitationalParameter = solGraviationalParameter+  }++marsOrbit :: Fractional a => Orbit a+marsOrbit = Orbit+  { eccentricity                  = 0.0934+  , periapsis                     = 1.382 % AstronomicalUnit+  , inclinationSpecifier          = Inclined+                                      { longitudeOfAscendingNode = 49.558 % [si|deg|]+                                      , inclination = 1.67 % [si|deg|]+                                      }+  , periapsisSpecifier = Eccentric { argumentOfPeriapsis = 286.502 % [si|deg|] }+  , primaryGravitationalParameter = solGraviationalParameter+  }++halleyOrbit :: Fractional a => Orbit a+halleyOrbit = Orbit+  { eccentricity                  = 0.96714+  , periapsis                     = 0.586 % AstronomicalUnit+  , inclinationSpecifier          = Inclined+                                      { longitudeOfAscendingNode = 58.42 % [si|deg|]+                                      , inclination = 162.26 % [si|deg|]+                                      }+  , periapsisSpecifier = Eccentric { argumentOfPeriapsis = 111.33 % [si|deg|] }+  , primaryGravitationalParameter = solGraviationalParameter+  }++++-- | The fastest comet in the west. Nice for testing as it's on a hyperbolic+-- trajectory. See https://en.wikipedia.org/wiki/C/1980_E1+--+-- Orbital data from:+-- http://ssd.jpl.nasa.gov/horizons.cgi?CGISESSID=6c2730c1201457522760d3f26b7d1f00#results+c1980E1Orbit :: Fractional a => Orbit a+c1980E1Orbit = Orbit+  { eccentricity                  = 1.057731876173255+  , periapsis                     = 3.363937831611605 % AstronomicalUnit+  , inclinationSpecifier          = Inclined+    { longitudeOfAscendingNode = 114.5581951921299 % [si|deg|]+    , inclination              = 1.661712630614323 % [si|deg|]+    }+  , periapsisSpecifier            = Eccentric+    { argumentOfPeriapsis = 135.0826233919265 % [si|deg|]+    }+  , primaryGravitationalParameter = solGraviationalParameter+  }++

@@ -0,0 +1,559 @@+{-# language QuasiQuotes #-}++module Physics.Orbit.StateVectors+  ( -- *** Types+    StateVectors(..)+  , Position+  , Velocity+    -- *** Conversion to state vectors+  , stateVectorsAtTrueAnomaly+  , positionAtTrueAnomaly+  , positionInPlaneAtTrueAnomaly+  , velocityAtTrueAnomaly+  , velocityInPlaneAtTrueAnomaly+    -- *** Conversion from state vectors+  , elementsFromStateVectors+  , eccentricityVector+  , trueAnomalyAtPosition+    -- *** Rotations to and from orbital plane+  , orbitalPlaneQuaternion+  , rotateToPlane+  , rotateFromPlane+    -- *** other utilities+  , flightPathAngleAtTrueAnomaly+  , specificAngularMomentumVector+  ) where++import           Control.Lens.Operators         ( (^.) )+import           Data.Coerce+import           Data.Constants.Mechanics.Extra+import           Data.Metrology+import           Data.Metrology.Extra+import           Data.Metrology.Unsafe          ( Qu(..) )+import           Data.Units.SI.Parser+import           Linear.Conjugate+import           Linear.Quaternion+import           Linear.V3+import           Physics.Orbit++type Position a = V3 (Distance a)+type Velocity a = V3 (Speed a)++data StateVectors a = StateVectors+  { position :: Position a+  , velocity :: Velocity a+  }+  deriving (Show, Eq)++----------------------------------------------------------------+-- Conversiont to state vectors+----------------------------------------------------------------++-- | Get the position in space of a body after rotating it according to the+-- inclination and periapsis specifier.+positionAtTrueAnomaly+  :: (Conjugate a, RealFloat a) => Orbit a -> Angle a -> Position a+positionAtTrueAnomaly o = rotateFromPlane o . positionInPlaneAtTrueAnomaly o++-- | Get the position of a body relative to the orbital plane+positionInPlaneAtTrueAnomaly+  :: (Ord a, Floating a) => Orbit a -> Angle a -> Position a+positionInPlaneAtTrueAnomaly o ν = r+ where+  radius = radiusAtTrueAnomaly o ν+  r      = V3 (qCos ν |*| radius) (qSin ν |*| radius) zero++-- | Get the velocity in space of a body after rotating it according to the+-- inclination and periapsis specifier.+velocityAtTrueAnomaly+  :: (Conjugate a, RealFloat a) => Orbit a -> Angle a -> Velocity a+velocityAtTrueAnomaly o = rotateFromPlane o . velocityInPlaneAtTrueAnomaly o++-- | The in-plane velocity of a body+velocityInPlaneAtTrueAnomaly+  :: (Ord a, Floating a) => Orbit a -> Angle a -> Velocity a+velocityInPlaneAtTrueAnomaly o ν = v+ where+  μ    = primaryGravitationalParameter o+  e    = eccentricity o+  r    = radiusAtTrueAnomaly o ν+  h    = specificAngularMomentum o+  cosν = qCos ν+  sinν = qSin ν+  vr   = μ |*| e |*| sinν |/| h+  vtA  = h |/| r+  v    = V3 (vr |*| cosν |-| vtA |*| sinν) (vr |*| sinν |+| vtA |*| cosν) zero++stateVectorsAtTrueAnomaly+  :: (Conjugate a, RealFloat a) => Orbit a -> Angle a -> StateVectors a+stateVectorsAtTrueAnomaly o ν = StateVectors r v+ where+  r = positionAtTrueAnomaly o ν+  v = velocityAtTrueAnomaly o ν++----------------------------------------------------------------+-- Conversion from state vectors+----------------------------------------------------------------++-- Thanks to https://downloads.rene-schwarz.com/download/M002-Cartesian_State_Vectors_to_Keplerian_Orbit_Elements.pdf+elementsFromStateVectors+  :: (Ord a, Floating a, Conjugate a, RealFloat a, Show a)+  => Quantity [si| m^3 s^-2 |] a+  -> StateVectors a+  -> (Orbit a, Angle a)+elementsFromStateVectors μ sv@(StateVectors r v) = (o, ν)+ where+  o     = Orbit e q inclinationSpecifier' periapsisSpecifier' μ++  h     = specificAngularMomentumVector sv+  n     = V3 (qNegate (h ^. _y)) (h ^. _x) zero++  e'    = eccentricityVector μ sv+  e     = qNorm e'+  eNorm = (recip e *) <$> e'++  aInv  = (2 |/| qNorm r) |-| (qQuadrance v |/| μ)+  a     = qRecip aInv+  q     = if aInv == zero -- parabolic trajectory+    then qQuadrance h |/| (2 |*| μ)+    else a |*| (1 - e)++  ν = if e == zero+    then -- fall back to the slower version if this is a circular orbit+         trueAnomalyAtPosition o r+    else+      let cosν = eNorm `qDot` qNormalize r+      in  if r `qDot` v >= zero then qArcCos cosν else turn |-| qArcCos cosν++  inclinationSpecifier' =+    let i    = qArcCos ((h ^. _z) |/| qNorm h)+        cosΩ = n ^. _x |/| qNorm n+        _Ω   = if n ^. _y >= zero then qArcCos cosΩ else turn |-| qArcCos cosΩ+    in  if h ^. _x == zero && h ^. _y == zero+          then NonInclined+          else Inclined _Ω i++  -- If the orbit is not inclined, ω is relative to the reference direction+  -- [1,0,0]+  periapsisSpecifier' =+    let cosω = case inclinationSpecifier' of+          Inclined _ _ -> qNormalize n `qDot` eNorm+          NonInclined  -> eNorm ^. _x+        -- ω = if (e' ^. _z) >= zero then qArcCos cosω else turn |-| qArcCos cosω+        ω = case inclinationSpecifier' of+          Inclined _ _ ->+            if (e' ^. _z) >= zero then qArcCos cosω else turn |-| qArcCos cosω+          NonInclined ->+            let sinω = eNorm ^. _y in qArcTan2 sinω cosω `mod'` turn+    in  if e == zero then Circular else Eccentric ω++-- | Calculate the true anomaly, ν, of a body at position, r, given its orbital+-- elements.+trueAnomalyAtPosition+  :: (Conjugate a, RealFloat a) => Orbit a -> Position a -> Angle a+trueAnomalyAtPosition o r = ν+ where+  V3 (Qu x) (Qu y) _ = rotateToPlane o r+  ν                  = atan2 y x % [si|rad|]++-- | Calculate the momentum vector, h, given state vectors+specificAngularMomentumVector+  :: Num a => StateVectors a -> V3 (Quantity [si|m^2 / s|] a)+specificAngularMomentumVector (StateVectors r v) = r `qCross` v++-- | Calculate the eccentricity vector, e, given state vectors+eccentricityVector+  :: Floating a+  => Quantity [si| m^3 s^-2 |] a+  -> StateVectors a+  -> V3 (Unitless a)+eccentricityVector μ sv@(StateVectors r v) = e+ where+  e = (v `qCross` h) |^/| μ |^-^| qNormalize r+  h = specificAngularMomentumVector sv++----------------------------------------------------------------+-- Rotations to and from the orbital plane+----------------------------------------------------------------++-- | Rotate a position relative to the orbital plane according to the+-- inclination specifier and periapsis specifier.+--+-- The orbital plane is perpendicular to the z axis+rotateFromPlane+  :: (Conjugate a, RealFloat a)+  => Orbit a+  -> V3 (Qu u l a)+  -> V3 (Qu u l a)+rotateFromPlane = qRotate . orbitalPlaneQuaternion++-- | Rotate a position such that is is relative to the orbital plane according+-- to the inclination specifier and periapsis specifier.+--+-- The orbital plane is perpendicular to the z axis+rotateToPlane+  :: (Conjugate a, RealFloat a) => Orbit a -> V3 (Qu u l a) -> V3 (Qu u l a)+rotateToPlane = qRotate . conjugate . orbitalPlaneQuaternion++-- | A quaternion representing the rotation of the orbital plane+orbitalPlaneQuaternion :: RealFloat a => Orbit a -> Quaternion a+orbitalPlaneQuaternion o = lon * per+ where+  per = case periapsisSpecifier o of+    Eccentric ω -> rotateZ ω+    Circular    -> 1+  lon = case inclinationSpecifier o of+    Inclined _Ω i -> rotateZ _Ω * rotateX i+    NonInclined   -> 1++----------------------------------------------------------------+-- Orbit Utils+----------------------------------------------------------------++-- | Get the flight path angle, φ, of a body a a specific true anomaly. This is+-- the angle of the body's motion relative to a vector perpendicular to the+-- radius.+flightPathAngleAtTrueAnomaly+  :: (Real a, Floating a) => Orbit a -> Angle a -> Angle a+flightPathAngleAtTrueAnomaly o ν = sign (qArcCos cosφ)+ where+  cosφ = h |/| (r |*| v)+  sign = if (ν `mod'` turn) < halfTurn then id else qNegate+  r    = radiusAtTrueAnomaly o ν+  v    = speedAtTrueAnomaly o ν+  h    = specificAngularMomentum o++----------------------------------------------------------------+-- Utils+----------------------------------------------------------------++qRotate+  :: forall a q+   . (Coercible (q a) a, Conjugate a, RealFloat a)+  => Quaternion a+  -> V3 (q a)+  -> V3 (q a)+qRotate = coerce (rotate @a)++rotateX :: Floating a => Angle a -> Quaternion a+rotateX θ = Quaternion (cos half) (V3 (sin half) 0 0)+  where half = (θ # [si|rad|]) / 2++_rotateY :: Floating a => Angle a -> Quaternion a+_rotateY θ = Quaternion (cos half) (V3 0 (sin half) 0)+  where half = (θ # [si|rad|]) / 2++rotateZ :: Floating a => Angle a -> Quaternion a+rotateZ θ = Quaternion (cos half) (V3 0 0 (sin half))+  where half = (θ # [si|rad|]) / 2++{-+++orbitalPlaneQuaternion :: Orbit -> Quaternion Double+orbitalPlaneQuaternion Elliptic{..} = l * p+  where p = case periapsisSpecifier of+              Eccentric ω -> rotateZ ω+              Circular -> noRotation+        l = case longitudeSpecifier of+              Inclined{..} -> rotateZ longitudeOfAscendingNode * rotateX inclination+              NonInclined -> noRotation++rotateToWorld :: Orbit -> V3 Double -> V3 Double+rotateToWorld orbit = rotate (orbitalPlaneQuaternion orbit)++rotateToPlane :: Orbit -> V3 Double -> V3 Double+rotateToPlane orbit = rotate (conjugate (orbitalPlaneQuaternion orbit))++positionAtTrueAnomaly :: Orbit -> Angle -> V3 Double+positionAtTrueAnomaly orbit trueAnomaly = rotateToWorld orbit r+  where ν = trueAnomaly+        d = radiusAtTrueAnomaly orbit ν+        r = V3 (cos ν) (sin ν) 0 ^* d++velocityAtTrueAnomaly :: Orbit -> Angle -> V3 Double+velocityAtTrueAnomaly orbit trueAnomaly = rotateToWorld orbit v+  where ν = trueAnomaly+        μ = primaryGravitationalParameter orbit+        e = eccentricity orbit+        h = sqrt (μ * a * (1 - e^2))+        a = semiMajorAxis orbit+        r = radiusAtTrueAnomaly orbit trueAnomaly+        vr = μ * e * sin ν / h+        vtA = h / r+        v = V3 (vr * cos ν - vtA * sin ν) (vr * sin ν + vtA * cos ν) 0++trueAnomalyAtPosition :: Orbit -> V3 Double -> Angle+trueAnomalyAtPosition orbit r = ν+  where V3 x y _ = rotateToPlane orbit r+        ν = atan2 y x++-- also equal to sqrt(μ/a^3)+averageAngularVelocity :: Orbit -> Angle+averageAngularVelocity orbit = 2 * pi / p+  where p = period orbit++distance :: Orbit -> Angle -> Distance+distance orbit@Elliptic{..} trueAnomaly = semilatusRectum orbit / (1 + eccentricity * cos trueAnomaly)+++{-+eccentricAnomaly :: Orbit -> Angle -> Angle+eccentricAnomaly Elliptic{..} trueAnomaly = acos ((eccentricity + cosTrue)/(1 + eccentricity * cosTrue))+  where cosTrue = cos trueAnomaly++meanAnomaly :: Orbit -> Angle -> Angle+meanAnomaly orbit@Elliptic{..} trueAnomaly = e - eccentricity * sin e+  where e = eccentricAnomaly orbit trueAnomaly+  -}++eccentricityVector :: Orbit -> Angle -> V3 Double+eccentricityVector orbit trueAnomaly = eccentricityVectorFromState μ sv+  where sv = stateVectorsFromOrbit orbit trueAnomaly+        μ = primaryGravitationalParameter orbit++eccentricityVectorFromState :: Double -> StateVectors -> V3 Double+eccentricityVectorFromState primaryGravitationalParameter StateVectors{..} = (v `cross` h) ^/ μ - normalize r+  where μ = primaryGravitationalParameter+        r = position+        v = velocity+        h = r `cross` v++trueAnomalyFromState :: Orbit -> StateVectors -> Angle+trueAnomalyFromState orbit stateVectors = if r `dot` v >= 0 then ν else 2 * pi - ν+  where e = eccentricityVectorFromState μ stateVectors+        r = position stateVectors+        v = velocity stateVectors+        ν = acos $ (e `dot` r) / (norm e * norm r)+        μ = primaryGravitationalParameter orbit++orbitalSpeed :: Orbit -> Angle -> Double+orbitalSpeed orbit trueAnomaly = v+  where d = Orbit.distance orbit trueAnomaly+        --ν = trueAnomaly+        μ = primaryGravitationalParameter orbit+        a = semiMajorAxis orbit+        v = if | isElliptic orbit ||+                 isHyperbolic orbit -> sqrt (μ * (2 / d - 1 / a))+               | isParabolic orbit -> sqrt (μ * 2 / d)++velocityAngleFromPrograde :: Orbit -> Angle -> Angle+velocityAngleFromPrograde orbit trueAnomaly = φ+  where ν = trueAnomaly+        e = eccentricity orbit+        φ = if | isElliptic orbit ||+                 isHyperbolic orbit -> atan2 (e * sin ν) (1 + e * cos ν)+               | isParabolic orbit -> ν / 2++-- | Calculate the state vectors relative to the orbital plane+--+-- The Z dimension is perpendicular to the orbital plane and hence is+-- always zero+orbitalPlaneStateVectors :: Orbit -> Angle -> StateVectors+orbitalPlaneStateVectors orbit trueAnomaly = StateVectors r v+  where d = Orbit.distance orbit trueAnomaly+        ν = trueAnomaly+        r = V3 (d * cos ν) (d * sin ν) 0+        e = eccentricity orbit+        --a = semiMajorAxis orbit+        u = V3 (1 + e * cos ν) (e * sin ν) 0+        v = u ^* orbitalSpeed orbit trueAnomaly+        --n = averageAngularVelocity orbit+        --v = V3 (- sin ν) (e + cos ν) 0 ^* (n * a / sqrt (1 - e^2))++rotateX :: Angle -> Quaternion Double+rotateX = axisAngle $ V3 1 0 0++rotateY :: Angle -> Quaternion Double+rotateY = axisAngle $ V3 0 1 0++rotateZ :: Angle -> Quaternion Double+rotateZ = axisAngle $ V3 0 0 1++noRotation :: Num a => Quaternion a+noRotation = Quaternion 1 (V3 0 0 0)++stateVectorsFromOrbit :: Orbit -> Angle -> StateVectors+stateVectorsFromOrbit orbit trueAnomaly = StateVectors r v+  where ν = trueAnomaly+        r = positionAtTrueAnomaly orbit ν+        v = velocityAtTrueAnomaly orbit ν+        {-+        o = orbitalPlaneStateVectors orbit trueAnomaly+        r' = position o+        v' = velocity o+        p = case periapsisSpecifier of+              Eccentric ω -> rotateZ ω+              Circular -> noRotation+        l = case longitudeSpecifier of+              Inclined{..} -> rotateZ longitudeOfAscendingNode * rotateX inclination+              NonInclined -> noRotation+        r = (l * p) `rotate` r'+        v = (l * p) `rotate` v'+        -}++orbitFromStateVectors :: StateVectors -> Double -> (Orbit, Angle)+orbitFromStateVectors sv@StateVectors{..} primaryGravitationalParameter = (orbit, ν)+  where r = position+        v = velocity+        μ = primaryGravitationalParameter+        -- `h` is the specific relative angular momentum+        h@(V3 _ _ hz) = r `cross` v+        -- `an` is the vector pointing towards the ascending node+        -- Todo, handle inclinations of 90 degrees here+        an@(V3 anx any _) = let an' = V3 0 0 1 `cross` h+                            in if nearZero (norm an') then V3 1 0 0 else an'+        -- `ev` is the eccentricity vector+        ev@(V3 evx evy evz) = eccentricityVectorFromState μ sv+        -- ε is the specific orbital energY+        ε = quadrance v / 2 - μ / norm r+        -- `a` is the semimajor axis+        --a = μ * norm r / (2 * μ - norm r * quadrance v)+        a = let a' = μ / (2 * ε)+            in if | isElliptic orbit ||+                    isHyperbolic orbit -> -a'+                  | isParabolic orbit -> error "parabolic orbits don't have a well defined semi-major axis"+        -- `e` is the eccentricity, Sometimes these numbers come out a tiny+        -- bit negative so clamp with 0+        --e = sqrt (max 0 $ 1 - quadrance h / (μ * a))+        e = norm ev+        -- `i` is the inclination+        i = acos $ hz / norm h+        -- `lan` is the longitude of the ascending node, sometimes known as Ω+        lan = let lan' = acos (anx / norm an)+              in if any >= 0 then lan' else 2 * pi - lan'+        -- `ω` is the argument of periapsis+        ω = let ω' = acos ((an `dot` ev)/(norm an * norm ev))+            in if evz < 0 then 2 * pi - ω' else ω'+        --ω = let ω' = atan2 evy evx+            --in if (r `cross` v < 0) then 2 * pi - ω' else ω'+        -- `ν` is the true anomaly+        ν = let ν' = acos ((ev `dot` r)/(norm ev * norm r))+            in if | isElliptic orbit -> if r `dot` v < 0 then 2 * pi - ν' else ν'+                  | isHyperbolic orbit -> if r `dot` v < 0 then -ν' else ν'+        orbit = Elliptic{ eccentricity = e+               , semiMajorAxis = a+               , longitudeSpecifier = if nearZero i then NonInclined+                                                    else Inclined { inclination = i+                                                                  , longitudeOfAscendingNode = lan+                                                                  }+               , periapsisSpecifier = if nearZero e then Circular+                                                    else Eccentric{argumentOfPeriapsis = ω}+               , primaryGravitationalParameter = μ}++lambert :: V3 Double -> V3 Double -> Double -> Double -> [(V3 Double, V3 Double)]+lambert r1 r2 primaryGravitationalParameter transferTime = [(v1, v2)]+  where μ = primaryGravitationalParameter+        h = r1 `cross` r2+        cosθ = (r1 `dot` r2) / (norm r1 * norm r2)+        θ = let θ' = acos cosθ+            in if | (h^._z) >= 0 -> θ' -- Todo, fixme+                  | otherwise -> 2 * π - θ'+        d = if | 0 <= θ && θ <= π -> 1+               | π < θ && θ <= 2 * π -> -1+        τ = d * sqrt (norm r1 * norm r2 * (1 + cosθ)) / (norm r1 + norm r2)+        s = sqrt $ ((norm r1 + norm r2)^3) / μ+        n = 0+        wse k = let v = k - sqrt 2+                in sqrt 2/3 - v/5 + 2/35*sqrt 2*v^2 - 2/63*v^3 + 2/231*sqrt 2*v^4 -+                   2/429*v^5 + 8/6435*sqrt 2*v^6 - 8/12155*v^7 + 8/46189*sqrt 2*v^8 -+                   8/88179*v^9 + 16/676039*sqrt 2*v^10 - 16/1300075*v^11 ++                   16/5014575*sqrt 2*v^12 - 16/9694845*v^13 ++                   128/300540195*sqrt 2*v^14 - 128/583401555*v^15 ++                   128/2268783825*sqrt 2*v^16+        tof n k = (tofk, tof'k, tof''k)+                  where tofk = s * sqrt (1 - k * τ) * (τ + (1 - k * τ) * w) -- 26+                        tof'k = -tofk / (2 * c) + s * τ * sqrt (c * τ) * (w' * c - w)+                        tof''k = -tofk / (4 * c^2) + s * τ * sqrt (c * τ) * (w / c + c * w'' - 3 * w')+                        c = (1 - k * τ) / τ+                        ε = 2e-2+                        w = if | k < sqrt 2 - ε ->+                                  ((1 - signum k) * π + signum k * acos (1 - m) + 2 * π * n) /+                                  sqrt (m^3) -+                                  k/m+                               | k > sqrt 2 + ε -> - acosh (1 - m) / sqrt (-m^3) - k / m+                               | otherwise -> ws -- 27+                        w' = if | k < sqrt 2 - ε -> (-2 + 3 * w * k) / m+                                | k > sqrt 2 + ε -> (-2 + 3 * w * k) / (-m)+                                | otherwise -> ws'+                        w'' = if | k < sqrt 2 - ε -> (5 * w' * k + 3 * w) / m+                                 | k > sqrt 2 + ε -> (5 * w' * k + 3 * w) / (-m)+                                 | otherwise -> ws''+                        (ws:ws':ws'':_) = diffs wse k+                        m = 2 - k^2+        --Right k = traceShowId $ newton (\k -> let (a, b, _) = tof (traceShowId k) in (a - transferTime, b)) (-sqrt 2) (sqrt 2) 1e-6+        initialGuess = 0+        --isValid = (&&) <$> (not . isNaN) <*> (-sqrt 2<)+        --ks = filter (isValid . snd) . zip initialGuesses $ halley (\k -> let (a,b,c) = tof n k in (a - transferTime, b, c)) <$> initialGuesses+        --Right k = newton (\k -> let (a,b,_) = tof n k in (a - transferTime, b)) (-sqrt 2) (sqrt 2) 1e-6+        --k = -1.414284878632464+        --k = snd . head $ ks+        k = halley (\k -> let (a,b,c) = tof n k in (a - transferTime, b, c)) initialGuess+        --kMinTime n = (\(Right y) -> y) $ newton (\k -> let (_, y',y'') = tof n k in (y', y'')) (-1) 1 1e-6+        --kbs = kMinTime <$> [1..]+        --tbs = (^._1) . uncurry tof <$> zip [1..] kbs+        f = 1 - (1 - k * τ) * (norm r1 + norm r2) / norm r1 -- 1 - (1 - k * τ) / norm r1+        --g' = 1 - (1 - k * τ) / norm r2+        g' = 1 - (1 - k * τ) * (norm r1 + norm r2) / norm r2+        --g = s * τ * sqrt ((1 - k * τ) * μ) -- τ * (norm r1 + norm r2) * sqrt (1 - k * τ)+        g = s * τ * sqrt (1 - k * τ)+        v1 = (r2 - f *^ r1) ^/ g -- ^* sqrt μ+        v2 = (g' *^ r2 - r1) ^/ g+        {-debugInfo = "kbs: " ++ show (take 5 kbs) +++                    "\ntbs: " ++ show (take 5 tbs) +++                    -- "\nks: " ++ show ks +++                    "\nd: " ++ show d +++                    "\nτ (tau): " ++ show τ +++                    "\nθ (theta): " ++ show θ +++                    "\nk: " ++ show k +++                    "\nn: " ++ show n +++                    "\ntof: " ++ show (tof n k) +++                    "\nf: " ++ show f +++                    "\ng: " ++ show g +++                    "\nr1: " ++ show r1 +++                    "\nr2: " ++ show r2 +++                    "\nv1: " ++ show v1 +++                    "\nv2: " ++ show v2 +++                    "\norbits: "-}+        (orbit1, ν1) = traceShowId $ Debug.Trace.trace debugInfo $ orbitFromStateVectors (StateVectors r1 v1) μ+        (orbit2, ν2) = traceShowId $ orbitFromStateVectors (StateVectors r2 v2) μ+        ma1 = meanAnomalyAtTrueAnomaly orbit1 ν1+        ma2 = meanAnomalyAtTrueAnomaly orbit2 ν2++isValid :: V3 Double -> Bool+isValid = noneOf each isNaN++ballisticTransfer :: (Double -> StateVectors) -> (Double -> StateVectors) -> Double -> Double -> Double -> Double -> (Burn, Burn)+ballisticTransfer fo1 fo2 primaryGravitationalParameter departureMin departureMax maxTransferTime = (b1, b2)+  where (b1, b2, _) = minimumBy (compare `on` (^._3)) ts+        ts = do let numDepartureSamples = 100+                    numArrivalSamples = 100+                    departureInterval = departureMax - departureMin+                d <- [0..numDepartureSamples-1]+                a <- [0..numArrivalSamples-1]+                let departureTime = departureMin + departureInterval * d / (numDepartureSamples - 1)+                    transferTime = maxTransferTime * a / (numArrivalSamples - 1)+                    arrivalTime = departureTime + transferTime+                    StateVectors{position = r1, velocity = v1} = fo1 departureTime+                    StateVectors{position = r2, velocity = v2} = fo2 arrivalTime+                (v1', v2') <- lambert r1 r2 primaryGravitationalParameter transferTime+                guard $ noneOf each isNaN v1'+                guard $ noneOf each isNaN v2'+                let b1 = v1' - v1+                    b2 = v2' - v2+                    δv1 = norm b1+                    δv2 = norm b2+                    δv = δv1 + δv2+                pure (Burn departureTime b1, Burn arrivalTime b2, δv)++toManoeuvreReferenceFrame :: Orbit -> Angle -> V3 Double -> V3 Double+toManoeuvreReferenceFrame orbit trueAnomaly = (m !*)+  where ν = trueAnomaly+        r = normalize $ positionAtTrueAnomaly orbit ν+        v = normalize $ velocityAtTrueAnomaly orbit ν+        prograde = normalize $ v+        normal = normalize $ prograde `cross` (-r)+        radial = normalize $ prograde `cross` normal+        m = V3 prograde normal radial++-}

@@ -6,10 +6,19 @@  import           Data.CReal import           GHC.TypeLits-import           Test.QuickCheck.Arbitrary (Arbitrary (..),-                                            arbitrarySizedFractional,-                                            shrinkRealFrac)+import           Linear.Conjugate+import           Linear.Epsilon+import           Test.QuickCheck.Arbitrary      ( Arbitrary(..)+                                                , arbitrarySizedFractional+                                                , shrinkRealFrac+                                                )  instance KnownNat n => Arbitrary (CReal n) where   arbitrary = arbitrarySizedFractional   shrink = shrinkRealFrac++instance TrivialConjugate (CReal n) where+instance Conjugate (CReal n) where++instance Epsilon (CReal n) where+  nearZero = const False

@@ -0,0 +1,140 @@+{-# language QuasiQuotes #-}++module Data.Metrology.Extra where++import           Control.Applicative+import           Data.Coerce                    ( coerce )+import           Data.Constants.Mechanics.Extra ( )+import qualified Data.Fixed                    as F+                                                ( div'+                                                , divMod'+                                                , mod'+                                                )+import           Data.Metrology+import           Data.Metrology.Unsafe          ( Qu(..) )+import           Data.Units.SI.Parser+import           Linear.Metric+import           Linear.V3+import           Linear.Vector+import           Physics.Orbit.Metrology++mod' :: forall a u l . Real a => Qu u l a -> Qu u l a -> Qu u l a+mod' = coerce (F.mod' :: a -> a -> a)++div'+  :: forall a b u v l+   . (Real a, Integral b)+  => Qu u l a+  -> Qu v l a+  -> Qu (Normalize (u @- v)) l b+div' = coerce (F.div' :: a -> a -> b)++divMod'+  :: forall a b u l+   . (Real a, Integral b)+  => Qu u l a+  -> Qu u l a+  -> (Qu '[] l b, Qu u l a)+divMod' = coerce (F.divMod' :: a -> a -> (b, a))++rad :: Fractional a => a -> Angle a+rad = (% [si|rad|])++rdh :: Fractional a => a -> AngleH a+rdh = (% RadianHyperbolic)++qCos :: Floating a => Angle a -> Unitless a+qCos θ = quantity $ cos (θ # [si|rad|])++qSin :: Floating a => Angle a -> Unitless a+qSin θ = quantity $ sin (θ # [si|rad|])++qTan :: Floating a => Angle a -> Unitless a+qTan θ = quantity $ tan (θ # [si|rad|])++qArcTan :: Floating a => Unitless a -> Angle a+qArcTan = rad . atan . (# [si||])++qArcTan2 :: RealFloat a => Unitless a -> Unitless a -> Angle a+qArcTan2 x y = rad (atan2 (x # [si||]) (y # [si||]))++qArcCos :: Floating a => Unitless a -> Angle a+qArcCos = rad . acos . (# [si||])++qRecip+  :: forall u l a . Fractional a => Qu u l a -> Qu (Normalize ('[] @- u)) l a+qRecip = coerce (recip @a)++qTanh :: Floating a => AngleH a -> Unitless a+qTanh = quantity . tanh . (# RadianHyperbolic)++qSinh :: Floating a => AngleH a -> Unitless a+qSinh = quantity . sinh . (# RadianHyperbolic)++qCosh :: Floating a => AngleH a -> Unitless a+qCosh = quantity . cosh . (# RadianHyperbolic)++qArcCosh :: Floating a => Unitless a -> AngleH a+qArcCosh = rdh . acosh . (# [si||])++qAbs :: forall a l u . Num a => Qu u l a -> Qu u l a+qAbs = coerce (abs @a)++qCross+  :: Num n+  => V3 (Qu a l n)+  -> V3 (Qu b l n)+  -> V3 (Qu (Normalize (a @@+ Reorder b a)) l n)+qCross (V3 a b c) (V3 d e f) =+  V3 (b |*| f |-| c |*| e) (c |*| d |-| a |*| f) (a |*| e |-| b |*| d)++qNorm :: forall u l a . Floating a => V3 (Qu u l a) -> Qu u l a+qNorm = coerce (norm @V3 @a)++-- qNormalize+--   :: forall u l a . (Floating a, Epsilon a) => V3 (Qu u l a) -> V3 (Qu '[] l a)+-- qNormalize = coerce (normalize @a @V3)+qNormalize+  :: Floating n+  => V3 (Qu b l n)+  -> V3+       ( Qu+           ( Normalize+               (Normalize ('[] @- b) @@+ Reorder b (Normalize ('[] @- b)))+           )+           l+           n+       )+qNormalize x = (qRecip (qNorm x) |*|) <$> x++qDot+  :: forall u v l a. Num a+  => V3 (Qu u l a)+  -> V3 (Qu v l a)+  -> Qu (Normalize (u @@+ Reorder v u)) l a+qDot = coerce (dot @V3 @a)++qQuadrance+  :: forall u l a+   . Num a+  => V3 (Qu u l a)+  -> Qu (Normalize (u @@+ Reorder u u)) l a+qQuadrance = coerce (quadrance @V3 @a)++(|^/|) :: (Functor f, Fractional n) =>+            f (Qu b l n)+            -> Qu u l n+            -> f (Qu+                    (Normalize+                       (Normalize ('[] @- u) @@+ Reorder b (Normalize ('[] @- u))))+                    l+                    n)+x |^/| y = (qRecip y |*|) <$> x++(|^-^|)+  :: forall f u l a+   . (Additive f, Applicative f, Num a)+  => f (Qu u l a)+  -> f (Qu u l a)+  -> f (Qu u l a)+(|^-^|) = liftA2 (|-|)

@@ -17,6 +17,7 @@                                                 )  newtype PositiveQuantity a = PositiveQuantity { getPositiveQuantity :: a }+  deriving(Show)  deriving instance Arbitrary a => Arbitrary (Qu u l a) @@ -24,7 +25,7 @@  instance (Num a, Ord a, Arbitrary a) => Arbitrary (PositiveQuantity (Qu u l a)) where   arbitrary = PositiveQuantity . Qu . getPositive <$> arbitrary-  shrink (PositiveQuantity x) = PositiveQuantity <$> shrink x+  shrink (PositiveQuantity (Qu x)) = [ PositiveQuantity (Qu 1) | x /= 1 ]  instance (Eq a) => EqProp (Qu u l a) where   (=-=) = eq

@@ -0,0 +1,18 @@+{-# OPTIONS_GHC -fno-warn-orphans #-}++module Linear.QuickCheck () where++import           Linear.V3+import           Test.QuickCheck+import           Test.QuickCheck.Checkers       ( EqProp(..)+                                                )++instance (EqProp a) => EqProp (V3 a) where+  (V3 x1 x2 x3) =-= (V3 y1 y2 y3) = (x1, x2, x3) =-= (y1, y2, y3)++instance Arbitrary a => Arbitrary (V3 a) where+  arbitrary = uncurry3 V3 <$> arbitrary+  shrink (V3 x y z) = fmap (uncurry3 V3) . shrink $ (x, y, z)++uncurry3 :: (a -> b -> c -> d) -> (a, b, c) -> d+uncurry3 f (a, b, c) = f a b c

@@ -1,6 +1,7 @@ {-# LANGUAGE DataKinds       #-} {-# LANGUAGE QuasiQuotes     #-} {-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE PatternSynonyms #-} {-# OPTIONS_GHC -fno-warn-orphans #-}  module Physics.Orbit.QuickCheck@@ -8,25 +9,38 @@   , EllipticOrbit(..)   , ParabolicOrbit(..)   , HyperbolicOrbit(..)+  , CanonicalOrbit(..)+  , pattern CircularOrbitF+  , pattern EllipticOrbitF+  , pattern ParabolicOrbitF+  , pattern HyperbolicOrbitF   , unitOrbit+  , overAllClasses   ) where +import           Data.Constants.Mechanics.Extra import           Data.Metrology-import           Data.Metrology.Unsafe+import           Data.Metrology.Extra           ( mod' ) import           Data.Metrology.QuickCheck+import           Data.Metrology.Unsafe import           Data.Units.SI.Parser+import           Linear.V3 import           Physics.Orbit                  ( Distance                                                 , InclinationSpecifier(..)                                                 , Orbit(..)                                                 , PeriapsisSpecifier(..)                                                 , Unitless                                                 )+import           Physics.Orbit.StateVectors import           System.Random                  ( Random ) import           Test.QuickCheck                ( Arbitrary(..)+                                                , Testable                                                 , choose                                                 , oneof                                                 , suchThat                                                 )+import           Test.Tasty                     ( TestTree )+import           Test.Tasty.QuickCheck          ( testProperty )  {-# ANN module ("HLint: ignore Reduce duplication" :: String) #-} @@ -42,6 +56,24 @@ newtype HyperbolicOrbit a = HyperbolicOrbit {getHyperbolicOrbit :: Orbit a}   deriving(Show, Eq) +-- | An orbit where all angles are in [0..2π) or [0..π)+--+-- Also not a weird orbit like circular or non inclined+newtype CanonicalOrbit a = CanonicalOrbit {getCanonicalOrbit :: Orbit a}+  deriving(Show, Eq)++pattern CircularOrbitF :: Orbit Float -> CircularOrbit Float+pattern CircularOrbitF o = CircularOrbit o++pattern EllipticOrbitF :: Orbit Float -> EllipticOrbit Float+pattern EllipticOrbitF o = EllipticOrbit o++pattern ParabolicOrbitF :: Orbit Float -> ParabolicOrbit Float+pattern ParabolicOrbitF o = ParabolicOrbit o++pattern HyperbolicOrbitF :: Orbit Float -> HyperbolicOrbit Float+pattern HyperbolicOrbitF o = HyperbolicOrbit o+ -- | Use aerobreaking to shrink an orbit without expending fuel instance (Num a, Ord a, Random a, Arbitrary a) => Arbitrary (Orbit a) where   arbitrary = oneof@@ -96,15 +128,29 @@       pure . HyperbolicOrbit $ Orbit { .. }   shrink (HyperbolicOrbit o) = HyperbolicOrbit <$> shrinkOrbit o +instance (Floating a, Real a, Random a, Arbitrary a) => Arbitrary (CanonicalOrbit a) where+  arbitrary = do+    PositiveQuantity eccentricity <- arbitrary+    PositiveQuantity periapsis    <- arbitrary+    PositiveQuantity _Ω           <- arbitrary+    PositiveQuantity i            <- arbitrary+    let inclinationSpecifier =+          Inclined (_Ω `mod'` turn) (i `mod'` (halfTurn |/| 2))+    ω <- arbitrary+    let periapsisSpecifier = Eccentric (ω `mod'` turn)+    PositiveQuantity primaryGravitationalParameter <- arbitrary+    pure . CanonicalOrbit $ Orbit { .. }+  -- shrink (CanonicalOrbit o) = CanonicalOrbit <$> shrinkOrbit o+ instance Arbitrary a => Arbitrary (InclinationSpecifier a) where   arbitrary = oneof [pure NonInclined, Inclined <$> arbitrary <*> arbitrary]-  shrink Inclined { .. } = [NonInclined]-  shrink NonInclined = []+  shrink Inclined {..} = [NonInclined]+  shrink NonInclined   = []  -- | The instance of Arbitrary for PeriapsisSpecifier doesn't generate Circular instance (Eq a, Num a, Arbitrary a) => Arbitrary (PeriapsisSpecifier a) where   arbitrary = Eccentric <$> arbitrary-  shrink (Eccentric x) = if x == zero then [] else [Eccentric zero]+  shrink (Eccentric x) = [Eccentric zero | x /= zero]   shrink Circular = []  --------------------------------------------------------------------------------@@ -137,8 +183,8 @@   :: (Num a, Eq a)   => MkQu_ULN [si|m^3 s^-2|] 'DefaultLCSU a   -> [MkQu_ULN [si|m^3 s^-2|] 'DefaultLCSU a]-shrinkPrimaryGravitationalParameter μ | μ == (Qu 1) = []-                                      | otherwise   = [Qu 1]+shrinkPrimaryGravitationalParameter μ | μ == Qu 1 = []+                                      | otherwise = [Qu 1]   --------------------------------------------------------------------------------@@ -152,3 +198,34 @@                  , periapsisSpecifier = Circular                  , primaryGravitationalParameter = 1 % [si|m^3 s^-2|]                  }+++----------------------------------------------------------------+-- Constructing test trees+----------------------------------------------------------------++overAllClasses+  :: (Random a, Arbitrary a, Num a, Ord a, Show a, Testable t)+  => (Orbit a -> t)+  -> [TestTree]+overAllClasses t =+  [ testProperty "circular"   (\(CircularOrbit o) -> t o)+  , testProperty "elliptic"   (\(EllipticOrbit o) -> t o)+  , testProperty "parabolic"  (\(ParabolicOrbit o) -> t o)+  , testProperty "hyperbolic" (\(HyperbolicOrbit o) -> t o)+  ]+++----------------------------------------------------------------+-- StateVectors+----------------------------------------------------------------++instance (Num a, Eq a, Arbitrary a) => Arbitrary (StateVectors a) where+  arbitrary =+    do+        r <- V3 <$> arbitrary <*> arbitrary <*> arbitrary+        v <- V3 <$> arbitrary <*> arbitrary <*> arbitrary+        pure $ StateVectors r v+      `suchThat` (\(StateVectors r v) ->+                   r /= V3 zero zero zero && v /= V3 zero zero zero+                 )

@@ -9,47 +9,38 @@   ( main   ) where -import           Control.Applicative            ( (<|>) ) import           Data.CReal                     ( CReal ) import           Data.CReal.QuickCheck          ( ) import           Data.Coerce                    ( coerce ) import           Data.Constants.Mechanics.Extra-import           Data.Maybe                     ( fromJust )+import           Data.Maybe import           Data.Metrology          hiding ( (%) ) import           Data.Metrology.Extra-import           Data.Proxy                     ( Proxy(..) )-import           Data.Ratio                     ( (%) )-import           Data.Tagged                    ( Tagged(..) ) import           Data.Units.SI.Parser-import           Numeric                        ( readFloat ) import           Physics.Orbit import           Physics.Orbit.QuickCheck import           Test.QuickCheck.Arbitrary      ( Arbitrary )-import           Test.QuickCheck.Checkers       ( inverse )+import           Test.QuickCheck.Checkers       ( inverse+                                                , inverseL+                                                )+import           Test.QuickCheck.Extra          ( slowTest+                                                , slowTestQCRatio+                                                ) import           Test.Tasty                     ( TestTree-                                                , adjustOption-                                                , askOption                                                 , defaultIngredients                                                 , defaultMainWithIngredients                                                 , includingOptions                                                 , testGroup                                                 )-import           Test.Tasty.Options             ( IsOption(..)-                                                , OptionDescription(..)-                                                ) import           Test.Tasty.QuickCheck          ( (===)                                                 , (==>)-                                                , QuickCheckTests(..)                                                 , testProperty                                                 ) import           Test.Tasty.TH                  ( testGroupGenerator )-import           Text.ParserCombinators.ReadP   ( char-                                                , eof-                                                , readP_to_S-                                                , readS_to_P-                                                ) import           WrappedAngle                   ( WrappedAngle(..) ) +import qualified Test.StateVectors+ {-# ANN module ("HLint: ignore Reduce duplication" :: String) #-}  -- | The type used for tests which require exact arithmetic. They are compared@@ -57,42 +48,6 @@ type Exact = CReal 32  ----------------------------------------------------------------------------------- Disable some really slow tests by default-----------------------------------------------------------------------------------newtype SlowTestQCRatio = SlowTestQCRatio Rational--slowTestQCRatio :: OptionDescription-slowTestQCRatio = Option (Proxy :: Proxy SlowTestQCRatio)--readRational :: String -> Maybe Rational-readRational s = case readP_to_S readRationalP s of-                   [(r,"")] -> Just r-                   _ -> Nothing-  where readRationalP = readS_to_P readFloat <* eof-                    <|> do n <- readS_to_P reads-                           _ <- char '/'-                           d <- readS_to_P reads-                           eof-                           pure (n%d)--instance IsOption SlowTestQCRatio where-  defaultValue = SlowTestQCRatio (1%10)-  parseValue = fmap SlowTestQCRatio . readRational-  optionName = Tagged "slow-test-ratio"-  optionHelp = Tagged $-    unwords [ "Some of the slow tests can take a long time to run; set this"-            , "flag to change the number of slow test QuickCheck test cases as"-            , "a proportion of the non-slow test number."-            ]--slowTest :: TestTree -> TestTree-slowTest t = askOption (\(SlowTestQCRatio r) ->-                          adjustOption (qcRatio r) t)-  where qcRatio r (QuickCheckTests n) =-          QuickCheckTests (floor (fromIntegral n * r))---------------------------------------------------------------------------------- -- The tests -------------------------------------------------------------------------------- @@ -330,6 +285,9 @@ (.:) :: (a -> b) -> (c -> d -> a) -> c -> d -> b f .: g = \x y -> f (g x y) +(~>) :: Bool -> Bool -> Bool+a ~> b = not a || b+ test_conversions :: [TestTree] test_conversions = [ conversionToTime                    , conversionToMeanAnomaly@@ -347,6 +305,9 @@       , testGroup "from true anomaly"                   (anomalyTimeConversionTests (fromJust .: timeAtTrueAnomaly)                                               "true anomaly")+      , testProperty "from true anomaly out of bounds parabolic"+        (\ν (ParabolicOrbitF o) ->+          validTrueAnomaly o ν ~> isJust (timeAtTrueAnomaly o ν))       ]      conversionToMeanAnomaly = let s = "mean anomaly" in testGroup ("conversion to " ++ s)@@ -377,9 +338,9 @@      conversionToTrueAnomaly = let s = "true anomaly" in testGroup ("conversion to " ++ s)       [ testGroup "from time"-                  (timeAnomalyConversionTests (fromJust .: trueAnomalyAtTime) s)+                  (timeAnomalyConversionTests trueAnomalyAtTime s)       , testGroup "from mean anomaly"-                  (anomalyConversionTests (fromJust .: trueAnomalyAtMeanAnomaly)+                  (anomalyConversionTests trueAnomalyAtMeanAnomaly                                           "mean anomaly"                                           s)       , testGroup "from eccentric anomaly"@@ -398,27 +359,62 @@             inverse (coerce (fromJust . meanAnomalyAtEccentricAnomaly (o :: Orbit Exact)) :: WrappedAngle Exact -> WrappedAngle Exact)                     (coerce (fromJust . eccentricAnomalyAtMeanAnomaly o))) +      , slowTest $ testProperty "mean hyperbolic inverse"+          (\(HyperbolicOrbit o) ->+            inverseL (fromJust . meanAnomalyAtHyperbolicAnomaly @Exact o)+                     (fromJust . hyperbolicAnomalyAtMeanAnomaly o))+       , slowTest $ testProperty "mean true inverse"           (\(EllipticOrbit o) ->             inverse (fromJust . meanAnomalyAtTrueAnomaly (o :: Orbit Exact))-                    (fromJust . trueAnomalyAtMeanAnomaly o))+                    (trueAnomalyAtMeanAnomaly o)) -      , slowTest $ testProperty "time true inverse"+      , slowTest $ testProperty "time true inverse elliptic"           (\(EllipticOrbit o) ->             inverse (fromJust . timeAtTrueAnomaly (o :: Orbit Exact))-                    (fromJust . trueAnomalyAtTime o))+                    (trueAnomalyAtTime o)) +      , slowTest $ testProperty "true time inverse parabolic"+          (\(ParabolicOrbit o) ->+            -- Use inverseL because there doesn't exist a time for every true+            -- anomaly+            inverseL (fromJust . timeAtTrueAnomaly (o :: Orbit Exact))+                     (trueAnomalyAtTime o)+                    )+       , testProperty "time eccentric inverse"           (\(EllipticOrbit o) ->             inverse (fromJust . timeAtEccentricAnomaly (o :: Orbit Exact))                     (fromJust . eccentricAnomalyAtTime o)) +      -- , slowTest $ testProperty "time hyperbolic inverse"+      --     (\(HyperbolicOrbit o) ->+      --       inverseL (fromJust . timeAtHyperbolicAnomaly @Exact o)+      --                (fromJust . hyperbolicAnomalyAtTime o))+       , testProperty "eccentric true inverse"           (\(EllipticOrbit o) ->             inverse (coerce (fromJust . eccentricAnomalyAtTrueAnomaly (o:: Orbit Exact)) :: WrappedAngle Exact -> WrappedAngle Exact)                     (fromJust . coerce (trueAnomalyAtEccentricAnomaly o)))++      , testProperty "hyperbolic true inverse"+          (\(HyperbolicOrbit o) ->+            inverseL (fromJust . hyperbolicAnomalyAtTrueAnomaly o)+                     (fromJust . trueAnomalyAtHyperbolicAnomaly @Exact o))       ] +test_anomalies :: [TestTree]+test_anomalies =+  [ slowTest $ testProperty+      "hyperbolic true"+      (\(HyperbolicOrbit o) _M ->+        let Just _H = hyperbolicAnomalyAtMeanAnomaly @Exact o _M+            ν       = trueAnomalyAtMeanAnomaly o _M+            e       = eccentricity o+        in  qCosh _H === (qCos ν + e) / (1 + e * qCos ν)+      )+  ]+ -- TODO: Put parabolic and hyperbolic tests here test_areal :: [TestTree] test_areal = [ testProperty "elliptic areal area"@@ -429,8 +425,119 @@                                         in area === p |*| arealVelocity o)              ] +test_orbitalEnergy :: [TestTree]+test_orbitalEnergy =+  [ testProperty "negative elliptical energy"+                 (\(EllipticOrbitF o) -> specificOrbitalEnergy o < zero)+  , testProperty "zero parabolic energy"+                 (\(ParabolicOrbitF o) -> specificOrbitalEnergy o === zero)+  , testProperty "positive hyperbolic energy"+                 (\(HyperbolicOrbitF o) -> specificOrbitalEnergy o > zero)+  , testGroup+    "potential + kinetic"+    (overAllClasses+      (\o ν ->+        specificOrbitalEnergy @Exact o+          === specificPotentialEnergyAtTrueAnomaly o ν+          |+| specificKineticEnergyAtTrueAnomaly o ν+      )+    )+  ]++test_radius :: [TestTree]+test_radius =+  [ testGroup+    "periapsis when ν = 0"+    (overAllClasses (\o -> radiusAtTrueAnomaly @Exact o zero === periapsis o))+  , testProperty+    "constant on circular"+    (\(CircularOrbitF o) ν -> radiusAtTrueAnomaly o ν === periapsis o)+  , testProperty+    "apoapsis when ν == π for elliptic"+    (\(EllipticOrbit o) ->+      radiusAtTrueAnomaly @Exact o halfTurn === fromJust (apoapsis o)+    )+  , testGroup+    "l when ν == π/2"+    (overAllClasses+      (\o -> radiusAtTrueAnomaly @Exact o (halfTurn |*| (-0.5))+        === semiLatusRectum o+      )+    )+  , testGroup+    "l when ν == -π/2"+    (overAllClasses+      (\o -> radiusAtTrueAnomaly @Exact o (halfTurn |*| (-0.5))+        === semiLatusRectum o+      )+    )+  , testProperty+    "from E"+    (\(EllipticOrbit o) ν ->+      let Just _E = eccentricAnomalyAtTrueAnomaly @Exact o ν+      in  radiusAtTrueAnomaly o ν+            === fromJust (semiMajorAxis o)+            |*| (1 - eccentricity o |*| qCos _E)+    )+  ]++test_speed :: [TestTree]+test_speed =+  [ testProperty+    "constant on circular"+    (\(CircularOrbitF o) ν ν' ->+      speedAtTrueAnomaly o ν === speedAtTrueAnomaly o ν'+    )+  , testProperty+    "zero at apex"+    (\(ParabolicOrbitF o) -> speedAtTrueAnomaly o halfTurn === zero)+  , testProperty+    "below escape velocity for elliptical"+    (\(EllipticOrbitF o) ν -> speedAtTrueAnomaly o ν < escapeVelocityAtDistance+      (primaryGravitationalParameter o)+      (radiusAtTrueAnomaly o ν)+    )+  , testProperty+    "escape velocity for parabolic"+    (\(ParabolicOrbitF o) ν ->+      speedAtTrueAnomaly o ν === escapeVelocityAtDistance+        (primaryGravitationalParameter o)+        (radiusAtTrueAnomaly o ν)+    )+  , testProperty+    "above escape velocity for hyperbolic"+    (\(HyperbolicOrbitF o) _M ->+      let ν = trueAnomalyAtMeanAnomaly o _M+      in  speedAtTrueAnomaly o ν > escapeVelocityAtDistance+            (primaryGravitationalParameter o)+            (radiusAtTrueAnomaly o ν)+    )+  ]++test_angularMomentum :: [TestTree]+test_angularMomentum =+  [ testProperty "negative elliptical energy"+                 (\(EllipticOrbitF o) -> specificOrbitalEnergy o < zero)+  , testProperty "zero parabolic energy"+                 (\(ParabolicOrbitF o) -> specificOrbitalEnergy o === zero)+  , testProperty "positive hyperbolic energy"+                 (\(HyperbolicOrbitF o) -> specificOrbitalEnergy o > zero)+  ]++test_stateVectors :: [TestTree]+test_stateVectors = [Test.StateVectors.tests]+ main :: IO () main = do   let is = includingOptions [slowTestQCRatio] : defaultIngredients   defaultMainWithIngredients is $(testGroupGenerator)++----------------------------------------------------------------+-- Orbit utils+----------------------------------------------------------------++validTrueAnomaly :: (Floating a, Ord a) => Orbit a -> Angle a -> Bool+validTrueAnomaly o ν = case hyperbolicDepartureAngle o of+  Nothing -> True+  Just d  -> qAbs ν < d 

@@ -1,9 +1,31 @@ module Test.QuickCheck.Extra   ( (<=!)   , (>=!)+  , slowTest+  , slowTestQCRatio   ) where -import Test.QuickCheck (Property, counterexample)+import           Control.Applicative            ( (<|>) )+import           Data.Proxy                     ( Proxy(..) )+import           Data.Ratio                     ( (%) )+import           Data.Tagged                    ( Tagged(..) )+import           Numeric                        ( readFloat )+import           Test.QuickCheck                ( Property+                                                , counterexample+                                                )+import           Test.Tasty                     ( TestTree+                                                , adjustOption+                                                , askOption+                                                )+import           Test.Tasty.Options             ( IsOption(..)+                                                , OptionDescription(..)+                                                )+import           Test.Tasty.QuickCheck          ( QuickCheckTests(..) )+import           Text.ParserCombinators.ReadP   ( char+                                                , eof+                                                , readP_to_S+                                                , readS_to_P+                                                )  infix 4 <=! (<=!) :: (Ord a, Show a) => a -> a -> Property@@ -12,4 +34,40 @@ infix 4 >=! (>=!) :: (Ord a, Show a) => a -> a -> Property x >=! y = counterexample (show x ++ " ≱ " ++ show y) (x >= y)++--------------------------------------------------------------------------------+-- Reduce the number of slow tests+--------------------------------------------------------------------------------++newtype SlowTestQCRatio = SlowTestQCRatio Rational++slowTestQCRatio :: OptionDescription+slowTestQCRatio = Option (Proxy :: Proxy SlowTestQCRatio)++readRational :: String -> Maybe Rational+readRational s = case readP_to_S readRationalP s of+                   [(r,"")] -> Just r+                   _ -> Nothing+  where readRationalP = readS_to_P readFloat <* eof+                    <|> do n <- readS_to_P reads+                           _ <- char '/'+                           d <- readS_to_P reads+                           eof+                           pure (n%d)++instance IsOption SlowTestQCRatio where+  defaultValue = SlowTestQCRatio (1%10)+  parseValue = fmap SlowTestQCRatio . readRational+  optionName = Tagged "slow-test-ratio"+  optionHelp = Tagged $+    unwords [ "Some of the slow tests can take a long time to run; set this"+            , "flag to change the number of slow test QuickCheck test cases as"+            , "a proportion of the non-slow test number."+            ]++slowTest :: TestTree -> TestTree+slowTest t = askOption (\(SlowTestQCRatio r) ->+                          adjustOption (qcRatio r) t)+  where qcRatio r (QuickCheckTests n) =+          QuickCheckTests (floor (fromIntegral n * r)) 

@@ -0,0 +1,201 @@+{-# language QuasiQuotes #-}++module Test.StateVectors where++import           Control.Lens.Operators         ( (^.) )+import           Data.CReal                     ( CReal )+import           Data.CReal.QuickCheck          ( )+import           Data.Constants.Mechanics.Extra+import           Data.Metrology+import           Data.Metrology.Extra+import           Data.Metrology.QuickCheck+import           Data.Units.SI.Parser+import           Linear.Metric+import           Linear.QuickCheck              ( )+import           Linear.V3+import           Test.QuickCheck.Checkers+import           Test.QuickCheck.Extra+import           Test.Tasty+import           Test.Tasty.QuickCheck+import           Test.Tasty.TH                  ( testGroupGenerator )++import           Physics.Orbit+import           Physics.Orbit.QuickCheck+import           Physics.Orbit.StateVectors++-- | The type used for tests which require exact arithmetic. They are compared+-- at a resolution of 2^16+type Exact = CReal 16++test_planeRotation :: [TestTree]+test_planeRotation =+  [ testProperty+      "plane rotation inverse"+      (\o -> inverse @(Position Exact) (rotateToPlane o) (rotateFromPlane o))+  ]++test_stateVectorInverse :: [TestTree]+test_stateVectorInverse =+  [ testProperty+    "state vector elements inverse"+    (\(PositiveQuantity μ) sv ->+      let (o, ν) = elementsFromStateVectors @Exact μ sv+          sv'    = stateVectorsAtTrueAnomaly o ν+      in  sv' === sv+    )+  , slowTest $ testProperty+    "elements state vector inverse"+    (\(CanonicalOrbit o) (PositiveQuantity ((`mod'` turn) -> ν)) ->+      let μ        = primaryGravitationalParameter @Exact o+          sv       = stateVectorsAtTrueAnomaly o ν+          (o', ν') = elementsFromStateVectors μ sv+      in  validTrueAnomaly o ν ==> (o', ν') === (o, ν)+    )+  , slowTest $ testProperty+    "elements state vector inverse 2"+    (\(normalizeOrbit -> o) (PositiveQuantity ((`mod'` turn) -> ν)) ->+      let μ        = primaryGravitationalParameter @Exact o+          sv       = stateVectorsAtTrueAnomaly o ν+          (o', ν') = elementsFromStateVectors μ sv+      in  validTrueAnomaly o ν ==> normalizeν (o', ν') === normalizeν (o, ν)+    )+  ]+ where+  normalizeν (o, ν) = case periapsisSpecifier o of+    Eccentric ω | eccentricity o == 0 ->+      (o { periapsisSpecifier = Circular }, (ν |+| ω) `mod'` turn)+    _ -> (o, ν `mod'` turn)++test_normalize :: [TestTree]+test_normalize =+  [ testProperty+    "state vectors invariant over normalize"+    (\o ν ->+      let oN = normalizeOrbit @Exact o+      in  stateVectorsAtTrueAnomaly o ν === stateVectorsAtTrueAnomaly oN ν+    )+  , testProperty+    "plane quaternion invariant over normalize"+    (\o ->+      let q1 = orbitalPlaneQuaternion @Exact o+          q2 = orbitalPlaneQuaternion (normalizeOrbit o)+      in  q1 === q2 .||. q1 === negate q2+    )+  ]++test_positionVelocity :: [TestTree]+test_positionVelocity =+  [ testProperty+    "position magnitude"+    (\o ν ->+      let r1 = fmap (# [si|m|]) . positionAtTrueAnomaly @Exact o $ ν+          r2 = (# [si|m|]) . radiusAtTrueAnomaly o $ ν+      in  r2 * r2 === quadrance r1+    )+  , testProperty+    "position in plane z"+    (\o ν ->+      let r = positionInPlaneAtTrueAnomaly @Float o ν in r ^. _z === zero+    )+  , testProperty+    "velocity magnitude"+    (\o ν ->+      let r1 = fmap (# [si|m/s|]) . velocityAtTrueAnomaly @Exact o $ ν+          r2 = (# [si|m/s|]) . speedAtTrueAnomaly o $ ν+      in  r2 * r2 === quadrance r1+    )+  , testProperty+    "velocity in plane z"+    (\o ν ->+      let v = velocityInPlaneAtTrueAnomaly @Float o ν in v ^. _z === zero+    )+  , testProperty+    "velocity at ν=0"+    (\o ->+      let v     = velocityInPlaneAtTrueAnomaly @Exact o zero+          speed = speedAtTrueAnomaly o zero+      in  v === V3 zero speed zero+    )+  , testProperty+    "velocity in circular orbit"+    (\(CircularOrbit o) ν ->+      let v     = velocityInPlaneAtTrueAnomaly @Exact o ν+          speed = speedAtTrueAnomaly o zero+      in  qNorm v === speed+    )+  , testProperty+    "velocity perpendicular to radius in circular orbit"+    (\(CircularOrbit o) ν ->+      let v = velocityInPlaneAtTrueAnomaly @Exact o ν+          r = positionInPlaneAtTrueAnomaly o ν+      in  v `qDot` r === zero+    )+  ]++test_flightPathAngle :: [TestTree]+test_flightPathAngle =+  [ testProperty+    "fpa circular orbit "+    (\(CircularOrbit o) ν ->+      let φ = flightPathAngleAtTrueAnomaly @Exact o ν in φ === zero+    )+  , testProperty+    "fpa angular momentum"+    (\o ν ->+      let φ = flightPathAngleAtTrueAnomaly @Exact o ν+          h = specificAngularMomentum o+          r = radiusAtTrueAnomaly o ν+          v = speedAtTrueAnomaly o ν+      in  h === r |*| v |*| qCos φ+    )+  , testProperty+    "fpa velocity direction"+    (\o ν ->+      let φ = flightPathAngleAtTrueAnomaly @Exact o ν+          r = (# [si|m|]) <$> positionInPlaneAtTrueAnomaly o ν+          v = (# [si|m/s|]) <$> velocityInPlaneAtTrueAnomaly o ν+      in  validTrueAnomaly o ν+            ==>   sin (φ # [si|rad|])+            ===   normalize r+            `dot` normalize v+    )+  ]++test_specificAngularMomentum :: [TestTree]+test_specificAngularMomentum =+  [ testProperty+    "momentum from vectors"+    (\o -> specificAngularMomentum @Exact o === specificAngularMomentumSV o)+  , testProperty+    "momentum vector length"+    (\o ν ->+      let sv = stateVectorsAtTrueAnomaly @Exact o ν+          h1 = specificAngularMomentumVector sv+          h2 = specificAngularMomentum o+      in  qNorm h1 === h2+    )+  ]++prop_specificAngularMomentum :: Orbit Exact -> Property+prop_specificAngularMomentum o =+  specificAngularMomentum o === specificAngularMomentumSV o++specificAngularMomentumSV+  :: (Ord a, Floating a) => Orbit a -> Quantity [si|m^2 s^-1|] a+specificAngularMomentumSV o = rx |*| vy |-| ry |*| vx+   where+    ν          = zero+    V3 rx ry _ = positionInPlaneAtTrueAnomaly o ν+    V3 vx vy _ = velocityInPlaneAtTrueAnomaly o ν++tests :: TestTree+tests = $(testGroupGenerator)++----------------------------------------------------------------+-- Orbit utils+----------------------------------------------------------------++validTrueAnomaly :: (Floating a, Ord a) => Orbit a -> Angle a -> Bool+validTrueAnomaly o ν = case hyperbolicDepartureAngle o of+  Nothing -> True+  Just d  -> qAbs ν < d