diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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` diff --git a/default.nix b/default.nix --- a/default.nix +++ b/default.nix @@ -1,4 +1,6 @@ -{ pkgs ? import { }, 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; diff --git a/orbits.cabal b/orbits.cabal --- a/orbits.cabal +++ b/orbits.cabal @@ -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 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 diff --git a/readme.md b/readme.md --- a/readme.md +++ b/readme.md @@ -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 diff --git a/src/Data/Metrology/Extra.hs b/src/Data/Metrology/Extra.hs --- a/src/Data/Metrology/Extra.hs +++ b/src/Data/Metrology/Extra.hs @@ -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 (|-|) diff --git a/src/Physics/Orbit.hs b/src/Physics/Orbit.hs --- a/src/Physics/Orbit.hs +++ b/src/Physics/Orbit.hs @@ -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 -- , 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 , -- 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) diff --git a/src/Physics/Orbit/Metrology.hs b/src/Physics/Orbit/Metrology.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Orbit/Metrology.hs @@ -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||] diff --git a/src/Physics/Orbit/Sol.hs b/src/Physics/Orbit/Sol.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Orbit/Sol.hs @@ -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 + } + + diff --git a/src/Physics/Orbit/StateVectors.hs b/src/Physics/Orbit/StateVectors.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Orbit/StateVectors.hs @@ -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 + +-} diff --git a/test/Data/CReal/QuickCheck.hs b/test/Data/CReal/QuickCheck.hs --- a/test/Data/CReal/QuickCheck.hs +++ b/test/Data/CReal/QuickCheck.hs @@ -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 diff --git a/test/Data/Metrology/Extra.hs b/test/Data/Metrology/Extra.hs new file mode 100644 --- /dev/null +++ b/test/Data/Metrology/Extra.hs @@ -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 (|-|) diff --git a/test/Data/Metrology/QuickCheck.hs b/test/Data/Metrology/QuickCheck.hs --- a/test/Data/Metrology/QuickCheck.hs +++ b/test/Data/Metrology/QuickCheck.hs @@ -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 diff --git a/test/Linear/QuickCheck.hs b/test/Linear/QuickCheck.hs new file mode 100644 --- /dev/null +++ b/test/Linear/QuickCheck.hs @@ -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 diff --git a/test/Physics/Orbit/QuickCheck.hs b/test/Physics/Orbit/QuickCheck.hs --- a/test/Physics/Orbit/QuickCheck.hs +++ b/test/Physics/Orbit/QuickCheck.hs @@ -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 + ) diff --git a/test/Test.hs b/test/Test.hs --- a/test/Test.hs +++ b/test/Test.hs @@ -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 diff --git a/test/Test/QuickCheck/Extra.hs b/test/Test/QuickCheck/Extra.hs --- a/test/Test/QuickCheck/Extra.hs +++ b/test/Test/QuickCheck/Extra.hs @@ -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)) diff --git a/test/Test/StateVectors.hs b/test/Test/StateVectors.hs new file mode 100644 --- /dev/null +++ b/test/Test/StateVectors.hs @@ -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