orbits (empty) → 0.3
raw patch · 17 files changed
+1550/−0 lines, 17 filesdep +QuickCheckdep +addep +basebuild-type:Customsetup-changed
Dependencies added: QuickCheck, ad, base, checkers, doctest, exact-real, orbits, random, tagged, tasty, tasty-quickcheck, tasty-th, units, units-defs
Files
- .gitignore +17/−0
- LICENSE +24/−0
- Setup.hs +6/−0
- changelog.md +6/−0
- default.nix +45/−0
- orbits.cabal +101/−0
- readme.md +51/−0
- src/Data/Constants/Mechanics/Extra.hs +35/−0
- src/Data/Metrology/Extra.hs +33/−0
- src/Physics/Orbit.hs +547/−0
- test/Data/CReal/QuickCheck.hs +15/−0
- test/Data/Metrology/QuickCheck.hs +31/−0
- test/Physics/Orbit/QuickCheck.hs +154/−0
- test/Test.hs +436/−0
- test/Test/QuickCheck/Extra.hs +15/−0
- test/WrappedAngle.hs +24/−0
- test/doctest/Doctests.hs +10/−0
+ .gitignore view
@@ -0,0 +1,17 @@+dist+cabal-dev+*.o+*.hi+*.chi+*.chs.h+*.dyn_o+*.dyn_hi+.hpc+.hsenv+.cabal-sandbox/+cabal.sandbox.config+*.prof+*.aux+*.hp+.stack-work/+result
+ LICENSE view
@@ -0,0 +1,24 @@+Copyright (c) 2020, Joe Hermaszewski+All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:+ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.+ * Redistributions in binary form must reproduce the above copyright+ notice, this list of conditions and the following disclaimer in the+ documentation and/or other materials provided with the distribution.+ * Neither the name of the <organization> nor the+ names of its contributors may be used to endorse or promote products+ derived from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE+DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,6 @@+module Main where++import Distribution.Extra.Doctest (defaultMainWithDoctests)++main :: IO ()+main = defaultMainWithDoctests "doctests"
+ changelog.md view
@@ -0,0 +1,6 @@+# Change Log++## WIP++## 0.3+ - Switch to `units` from `uom-plugin`
+ default.nix view
@@ -0,0 +1,45 @@+{ pkgs ? import <nixpkgs> { }, compiler ? null, hoogle ? true }:++let+ src = pkgs.nix-gitignore.gitignoreSource [ ] ./.;++ compiler' = if compiler != null then+ compiler+ else+ "ghc" + pkgs.lib.concatStrings+ (pkgs.lib.splitVersion pkgs.haskellPackages.ghc.version);++ # Any overrides we require to the specified haskell package set+ haskellPackages = with pkgs.haskell.lib;+ pkgs.haskell.packages.${compiler'}.override {+ 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")+ { };+ } // pkgs.lib.optionalAttrs hoogle {+ ghc = super.ghc // { withPackages = super.ghc.withHoogle; };+ ghcWithPackages = self.ghc.withPackages;+ };+ };++ # Any packages to appear in the environment provisioned by nix-shell+ extraEnvPackages = with haskellPackages; [ ];++ # Generate a haskell derivation using the cabal2nix tool on `package.yaml`+ drv = let old = haskellPackages.callCabal2nix "" src { };+ in old // {+ # Insert the extra environment packages into the environment generated by+ # cabal2nix+ env = pkgs.lib.overrideDerivation old.env (attrs:+ {+ buildInputs = attrs.buildInputs ++ extraEnvPackages;+ } // pkgs.lib.optionalAttrs hoogle {+ shellHook = attrs.shellHook + ''+ export HIE_HOOGLE_DATABASE="$(cat $(${pkgs.which}/bin/which hoogle) | sed -n -e 's|.*--database \(.*\.hoo\).*|\1|p')"+ '';+ });+ };++in if pkgs.lib.inNixShell then drv.env else drv
+ orbits.cabal view
@@ -0,0 +1,101 @@+cabal-version: 1.24++-- This file has been generated from package.yaml by hpack version 0.33.0.+--+-- see: https://github.com/sol/hpack+--+-- hash: 19ff836740d3a6d31f4b16ef06879a924c919fabab962046e35730122d362dff++name: orbits+version: 0.3+synopsis: Types and functions for Kepler orbits.+category: Physics+homepage: https://github.com/expipiplus1/orbit#readme+bug-reports: https://github.com/expipiplus1/orbit/issues+author: Joe Hermaszewski+maintainer: Joe Hermaszewski <keep.it.real@monoid.al>+copyright: 2020 Joe Hermaszewski+license: BSD3+license-file: LICENSE+build-type: Custom+extra-source-files:+ .gitignore+ readme.md+ default.nix+ changelog.md++source-repository head+ type: git+ location: https://github.com/expipiplus1/orbit++custom-setup+ setup-depends:+ Cabal+ , base+ , cabal-doctest >=1 && <1.1++library+ exposed-modules:+ Data.Constants.Mechanics.Extra+ Physics.Orbit+ other-modules:+ Data.Metrology.Extra+ hs-source-dirs:+ src+ default-extensions: DataKinds GeneralizedNewtypeDeriving QuasiQuotes ScopedTypeVariables TypeOperators+ ghc-options: -Wall -O2+ build-depends:+ ad >=4.3.2+ , base >=4.8 && <5+ , exact-real >=0.12+ , units+ , units-defs >=2.2+ default-language: Haskell2010++test-suite doctests+ type: exitcode-stdio-1.0+ main-is: Doctests.hs+ other-modules:+ + hs-source-dirs:+ test/doctest+ default-extensions: DataKinds GeneralizedNewtypeDeriving QuasiQuotes ScopedTypeVariables TypeOperators+ ghc-options: -Wall+ build-depends:+ base+ , doctest+ default-language: Haskell2010++test-suite test+ type: exitcode-stdio-1.0+ main-is: Test.hs+ other-modules:+ Data.CReal.QuickCheck+ Data.Metrology.QuickCheck+ Physics.Orbit.QuickCheck+ Test.QuickCheck.Extra+ 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+ build-depends:+ QuickCheck+ , ad+ , base+ , checkers+ , exact-real+ , orbits+ , random+ , tagged+ , tasty+ , tasty-quickcheck+ , tasty-th+ , units+ , units-defs+ default-language: Haskell2010
+ readme.md view
@@ -0,0 +1,51 @@+orbit+=====++*For my uncle Zbys who watched the planets and stars.*++-----++Types and functions for dealing with Kepler orbits.++The main data type is `Orbit`, which describes the path of a body in orbit.++Nomenclature+------------++| Symbol | Meaning | Notes |+|--------|----------------------------------|--------------------------------|+| a | Semi-major axis | Negative for hyperbolic orbits |+| b | Semi-minor axis | Negative for hyperbolic orbits |+| e | Eccentricity | |+| q | Periapsis | |+| i | Inclination | |+| μ | Standard gravitational parameter | |+| Ω | Longitude of the ascending node | |+| l | Semi-latus Rectum | |+| n | Mean motion | |+| p | Period | |+| t | Time since periapse | |+| M | Mean anomaly | |+| E | Eccentric anomaly | Only for elliptic orbits |+| ν | True anomaly | |+++Note that in the Haskell source uppercase symbols such as Ω and M are written+with a leading underscore.++Implementation+--------------++This package makes use of the+[`units`](https://hackage.haskell.org/package/units) package to ensure that the+implementation is correct regarding units of measure.++Contributing+------------++Contributions and bug reports are welcome!++Please feel free to contact me on GitHub or as "jophish" on freenode.++-Joe+
+ src/Data/Constants/Mechanics/Extra.hs view
@@ -0,0 +1,35 @@+module Data.Constants.Mechanics.Extra+ ( turn+ , halfTurn+ , addRad+ , delRad+ ) where++import Data.Coerce+import qualified Data.Dimensions.SI as D+import Data.Metrology+import Data.Metrology.SI ( )+import Data.Metrology.Unsafe+import Data.Units.SI++type PlaneAngle = MkQu_DLN D.PlaneAngle 'DefaultLCSU++-- | One complete revolution in radians+turn :: Floating a => PlaneAngle a+turn = 2 |*| halfTurn++-- | π radians+halfTurn :: Floating a => PlaneAngle a+halfTurn = pi % Radian++-- | Multiply by 1 radian+addRad+ :: Qu b 'DefaultLCSU a+ -> Qu (Normalize ('[ 'F D.PlaneAngle One] @+ b)) 'DefaultLCSU a+addRad = coerce++-- | Divide by 1 radian+delRad+ :: Qu u 'DefaultLCSU a+ -> Qu (Normalize (u @- '[ 'F D.PlaneAngle One])) 'DefaultLCSU a+delRad = coerce
+ src/Data/Metrology/Extra.hs view
@@ -0,0 +1,33 @@+module Data.Metrology.Extra+ ( mod'+ , div'+ , divMod'+ ) where++import Data.Coerce ( coerce )+import qualified Data.Fixed as F+ ( div'+ , divMod'+ , mod'+ )+import Data.Metrology+import Data.Metrology.Unsafe ( Qu(..) )++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))
+ src/Physics/Orbit.hs view
@@ -0,0 +1,547 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE QuasiQuotes #-}+{-# LANGUAGE TypeFamilies #-}+{-# OPTIONS_GHC -fno-warn-orphans #-}++-- | Types and functions for dealing with Kepler orbits.+module Physics.Orbit+ ( -- * The Orbit data type and dependencies+ Orbit(..)+ , InclinationSpecifier(..)+ , PeriapsisSpecifier(..)+ , Classification(..)++ -- * Functions for dealing with orbits+ -- ** Utilities+ , isValid+ , classify+ -- ** Orbital elements+ , apoapsis+ , meanMotion+ , period+ , arealVelocity+ -- *** Geometry+ , semiMajorAxis+ , semiMinorAxis+ , semiLatusRectum+ , hyperbolicApproachAngle+ , hyperbolicDepartureAngle+ -- ** Conversions++ -- *** To time since periapse+ , timeAtMeanAnomaly+ , timeAtEccentricAnomaly+ , timeAtTrueAnomaly++ -- *** To mean anomaly+ , meanAnomalyAtTime+ , meanAnomalyAtEccentricAnomaly+ , meanAnomalyAtTrueAnomaly++ -- *** To eccentric anomaly+ , eccentricAnomalyAtTime+ , eccentricAnomalyAtMeanAnomaly+ , eccentricAnomalyAtMeanAnomalyFloat+ , eccentricAnomalyAtTrueAnomaly++ -- *** To true anomaly+ , trueAnomalyAtTime+ , trueAnomalyAtMeanAnomaly+ , trueAnomalyAtEccentricAnomaly++ -- * Unit synonyms+ , Time+ , Distance+ , Speed+ , Mass+ , Angle+ , Unitless++ -- * Reexported from 'Data.CReal'+ , Converge+ ) where++import Control.Monad ( (<=<) )+import Data.Bifunctor ( bimap+ , second+ )+import Data.CReal.Converge ( Converge+ , convergeErr+ )+import Data.Constants.Mechanics.Extra+import Data.Metrology+import Data.Metrology.Extra+import Data.Metrology.Show ( )+import Data.Metrology.Unsafe ( UnsafeQu(..) )+import Data.Units.SI.Parser+import Numeric.AD ( Mode+ , Scalar+ , auto+ )+import Numeric.AD.Halley ( findZero+ , findZeroNoEq+ )+import Numeric.AD.Internal.Identity ( Id(..) )++--------------------------------------------------------------------------------+-- 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.+ --+ -- 'eccentricity' must be non-negative.+ --+ -- An eccentricity of 0 describes a circular orbit.+ --+ -- An eccentricity of less than 1 describes an elliptic+ -- orbit.+ --+ -- An eccentricity equal to 1 describes a parabolic orbit.+ --+ -- An eccentricity greater than 1 describes a hyperbolic+ -- orbit.+ eccentricity :: !(Unitless a)+ -- | The orbit's periapsis, q.+ --+ -- 'periapsis' must be positive.+ --+ -- The periapsis is the distance between the bodies at+ -- their closest approach.+ , periapsis :: !(Distance a)+ -- | The 'inclinationSpecifier' describes the angle+ -- between the obtital plane and the reference plane.+ , inclinationSpecifier :: !(InclinationSpecifier a)+ -- | 'periapsisSpecifier' is 'Circular' iff+ -- 'eccentricity' is 0+ --+ -- The periapsis specifier describes any rotation of+ -- the orbit relative to the reference direction in the+ -- orbital plane.+ , periapsisSpecifier :: !(PeriapsisSpecifier a)+ -- | The gravitational parameter of the system's+ -- primary, μ.+ --+ -- μ is equal to the mass of the primary times+ -- <https://en.wikipedia.org/wiki/Gravitational_constant+ -- G>.+ --+ -- 'primaryGravitationalParameter' must be positive.+ , primaryGravitationalParameter :: !(Quantity [si| m^3 s^-2 |] a)+ }+ deriving (Show, Eq)++-- | Along with 'PeriapsisSpecifier' the 'InclinationSpecifier' describes+-- orbital elements extra to its geometry.+data InclinationSpecifier a = -- | The orbit does not lie exactly in the+ -- reference plane+ Inclined { -- | The longitude of the ascending+ -- node, Ω.+ --+ -- The angle between the reference+ -- direction and the point where the+ -- orbiting body crosses the reference+ -- plane in the positive z direction.+ longitudeOfAscendingNode :: !(Angle a)+ -- | The orbit's inclination, i.+ --+ -- The angle between the reference+ -- plane and the orbital plane+ , inclination :: !(Angle a)+ }+ -- | The orbit lies in the reference plane+ | NonInclined+ deriving (Show, Eq)++-- | Along with 'InclinationSpecifier' the 'PeriapsisSpecifier' describes+-- orbital elements extra to its geometry.+data PeriapsisSpecifier a = -- | The orbit is not circular+ Eccentric { -- | The argument of periapsis, ω.+ --+ -- The 'argumentOfPeriapsis' is the+ -- angle of the periapsis relative to+ -- the reference direction in the+ -- orbital plane.+ 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+-- 'eccentricity', e >= 0, of the orbit.+data Classification = -- | 0 <= e < 1+ --+ -- This includes circular orbits.+ Elliptic+ -- | e == 1+ | Parabolic+ -- | e > 1+ | Hyperbolic+ deriving (Show, Read, Eq)++-- TODO, use the neat "UnsafeQu" newtype for unsafe instances+unsafeMapUnit :: (a -> b) -> Qu u l a -> Qu u l b+unsafeMapUnit f = qu . fmap f . UnsafeQu++unsafeMapOrbit :: (a -> b) -> Orbit a -> Orbit b+unsafeMapOrbit f (Orbit e q i p μ) = Orbit (unsafeMapUnit f e)+ (unsafeMapUnit f q)+ (unsafeMapInclinationSpecifier f i)+ (unsafeMapPeriapsisSpecifier f p)+ (unsafeMapUnit f μ)++unsafeMapInclinationSpecifier :: (a -> b)+ -> InclinationSpecifier a -> InclinationSpecifier b+unsafeMapInclinationSpecifier f s = case s of+ Inclined _Ω i -> Inclined (unsafeMapUnit f _Ω) (unsafeMapUnit f i)+ NonInclined -> NonInclined++unsafeMapPeriapsisSpecifier :: (a -> b)+ -> PeriapsisSpecifier a -> PeriapsisSpecifier b+unsafeMapPeriapsisSpecifier f p = case p of+ Circular -> Circular+ Eccentric a -> Eccentric (unsafeMapUnit f a)++--------------------------------------------------------------------------------+-- 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.+isValid :: (Ord a, Num a) => Orbit a -> Bool+isValid o = e >= 0 &&+ ((e == 0) `iff` (periapsisSpecifier o == Circular)) &&+ q > zero &&+ μ > zero+ where+ iff = (==) :: Bool -> Bool -> Bool+ e = eccentricity o+ q = periapsis o+ μ = primaryGravitationalParameter o++-- | 'classify' is a funciton which returns the orbit's class.+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++-- | Calculate the semi-major axis, a, of the 'Orbit'. Returns 'Nothing' when+-- given a parabolic orbit for which there is no semi-major axis. Note that the+-- semi-major axis of a hyperbolic orbit is negative.+semiMajorAxis :: (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)+semiMajorAxis o =+ case classify o of+ Parabolic -> Nothing+ _ -> Just $ q |/| (1 |-| e)+ where+ q = periapsis o+ e = eccentricity o++-- | Calculate the semi-minor axis, b, of the 'Orbit'. Like 'semiMajorAxis'+-- @\'semiMinorAxis\' o@ is negative when @o@ is a hyperbolic orbit. In the+-- case of a parabolic orbit 'semiMinorAxis' returns 0m.+semiMinorAxis :: (Floating a, Ord a) => Orbit a -> Distance a+semiMinorAxis o =+ case classify o of+ Elliptic -> a |*| qSqrt (1 |-| e ^ (2::Int))+ Parabolic -> zero+ Hyperbolic -> a |*| qSqrt (e ^ (2::Int) |-| 1)+ where+ e = eccentricity o+ Just a = semiMajorAxis o++-- | Calculate the semiLatusRectum, l, of the 'Orbit'+semiLatusRectum :: (Num a) => Orbit a -> Distance a+semiLatusRectum orbit = e |*| q |+| q+ where q = periapsis orbit+ e = eccentricity orbit++-- | Calculate the distance between the bodies when they are at their most+-- distant. 'apoapsis' returns 'Nothing' when given a parabolic or hyperbolic+-- orbit.+apoapsis :: (Fractional a, Ord a) => Orbit a -> Maybe (Distance a)+apoapsis o =+ case classify o of+ Elliptic -> Just $ a |*| (1 |+| e)+ _ -> Nothing+ where+ Just a = semiMajorAxis o+ e = eccentricity o++-- | Calculate the mean motion, n, of an orbit+--+-- This is the rate of change of the mean anomaly with respect to time.+meanMotion :: (Floating a, Ord a) => Orbit a -> Quantity [si|rad/s|] a+meanMotion o =+ case classify o of+ Elliptic -> addRad $ qSqrt (μ |/| qCube a)+ Hyperbolic -> addRad $ qSqrt (μ |/| qNegate (qCube a))+ Parabolic -> addRad $ 2 |*| qSqrt (μ |/| qCube l)+ where+ Just a = semiMajorAxis o+ μ = primaryGravitationalParameter o+ l = semiLatusRectum o++-- | Calculate the orbital period, p, of an elliptic orbit.+--+-- 'period' returns Nothing if given a parabolic or hyperbolic orbit.+period :: (Floating a, Ord a) => Orbit a -> Maybe (Time a)+period o =+ case classify o of+ Elliptic -> Just p+ _ -> Nothing+ where+ n = meanMotion o+ p = turn |/| n+++-- | Calculate the areal velocity, A, of the orbit.+--+-- The areal velocity is the area <https://xkcd.com/21/ swept out> by the line+-- between the orbiting body and the primary per second.+arealVelocity :: (Ord a, Floating a) => Orbit a -> Quantity [si|m^2/s|] a+arealVelocity o = qSqrt (l |*| μ) |/| 2+ where l = semiLatusRectum o+ μ = primaryGravitationalParameter o++-- | Calculate the angle at which a body leaves the system when on an escape+-- trajectory relative to the argument of periapsis. This is the limit of the+-- true anomaly as time tends towards infinity minus the argument of periapsis.+-- The departure angle is in the closed range (π/2..π).+--+-- This is the negation of the approach angle.+--+-- 'hyperbolicDepartureAngle' returns Nothing when given an elliptic orbit and+-- π when given a parabolic orbit.+hyperbolicDepartureAngle :: (Floating a, Ord a) => Orbit a -> Maybe (Angle a)+hyperbolicDepartureAngle o =+ case classify o of+ Hyperbolic ->+ let e = eccentricity o+ θ = addRad $ acos (-1 / e)+ in Just θ+ Parabolic -> Just (turn |/| 2)+ _ -> Nothing++-- | Calculate the angle at which a body leaves the system when on a hyperbolic+-- trajectory relative to the argument of periapsis. This is the limit of the+-- true anomaly as time tends towards -infinity minus the argument of+-- periapsis. The approach angle is in the closed range (-π..π/2).+--+-- This is the negation of the departure angle.+--+-- 'hyperbolicApproachAngle' returns Nothing when given a non-hyperbolic orbit+-- and -π when given a parabolic orbit.+hyperbolicApproachAngle :: (Floating a, Ord a) => Orbit a -> Maybe (Angle a)+hyperbolicApproachAngle = fmap qNegate . hyperbolicDepartureAngle++-- | Calculate the time since periapse, t, when the body has the given+-- <https://en.wikipedia.org/wiki/Mean_anomaly mean anomaly>, M. M may be+-- negative, indicating that the orbiting body has yet to reach periapse.+--+-- The sign of the time at mean anomaly M is the same as the sign of M.+--+-- The returned time is unbounded.+timeAtMeanAnomaly :: (Floating a, Ord a) => Orbit a -> Angle a -> Time a+timeAtMeanAnomaly o _M = _M |/| n+ where n = meanMotion o++-- | Calculate the time since periapse, t, of an elliptic orbit when at+-- eccentric anomaly E.+--+-- 'timeAtEccentricAnomaly' returns Nothing if given a parabolic or hyperbolic+-- orbit.+timeAtEccentricAnomaly :: (Floating a, Ord a) => Orbit a -> Angle a -> Maybe (Time a)+timeAtEccentricAnomaly o = fmap (timeAtMeanAnomaly o) . meanAnomalyAtEccentricAnomaly 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++-- | Calculate the <https://en.wikipedia.org/wiki/Mean_anomaly mean anomaly>,+-- M, at the given time since periapse, t. t may be negative, indicating that+-- the orbiting body has yet to reach periapse.+--+-- The sign of the mean anomaly at time t is the same as the sign of t.+--+-- The returned mean anomaly is unbounded.+meanAnomalyAtTime :: (Floating a, Ord a) => Orbit a -> Time a -> Angle a+meanAnomalyAtTime o t = t |*| n+ where n = meanMotion o++-- | Calculate the mean anomaly, M, of an elliptic orbit when at eccentric+-- anomaly E+--+-- 'meanAnomalyAtEccentricAnomaly' returns Nothing if given a parabolic or+-- hyperbolic orbit.+--+-- The number of orbits represented by the anomalies is preserved;+-- i.e. M `div` 2π = E `div` 2π+meanAnomalyAtEccentricAnomaly :: (Floating a, Ord a) => Orbit a -> Angle a -> Maybe (Angle a)+meanAnomalyAtEccentricAnomaly o _E = case classify o of+ Elliptic -> Just _M+ _ -> Nothing+ where e = eccentricity o+ untypedE = delRad _E+ _M = addRad (untypedE |-| e |*| sin untypedE)++-- | 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.+meanAnomalyAtTrueAnomaly :: (Real a, Floating a)+ => Orbit a -> Angle a -> Maybe (Angle a)+meanAnomalyAtTrueAnomaly o = case classify o of+ Elliptic -> meanAnomalyAtEccentricAnomaly o <=<+ eccentricAnomalyAtTrueAnomaly o+ _ -> error "TODO: meanAnomalyAtTrueAnomaly"++-- | Calculate the eccentric anomaly, E, of an elliptic orbit at time t.+--+-- 'eccentricAnomalyAtTime' returns Nothing when given a parabolic or+-- hyperbolic orbit.+--+-- The number of orbits represented by the time is preserved;+-- i.e. t `div` p = E `div` 2π+eccentricAnomalyAtTime :: (Converge [a], Floating a, Real a)+ => Orbit a -> Time a -> Maybe (Angle a)+eccentricAnomalyAtTime o t = case classify o of+ Elliptic -> eccentricAnomalyAtMeanAnomaly o . meanAnomalyAtTime o $ t+ _ -> Nothing++-- | Calculate the eccentric anomaly, E, of an elliptic orbit when at mean+-- anomaly M. This function is considerably slower than most other conversion+-- functions as it uses an iterative method as no closed form solution exists.+--+-- The number of orbits represented by the anomalies is preserved;+-- i.e. M `div` 2π = E `div` 2π+--+-- 'eccentricAnomalyAtMeanAnomaly' returns Nothing when given a parabolic or+-- hyperbolic orbit.+eccentricAnomalyAtMeanAnomaly :: forall a. (Converge [a], Floating a, Real a)+ => Orbit a -> Angle a -> Maybe (Angle a)+eccentricAnomalyAtMeanAnomaly o _M = case classify o of+ Elliptic -> _E+ _ -> Nothing+ where (n, wrappedM) = second (# [si|rad|]) (_M `divMod'` turn)+ e = eccentricity o # [si||]+ _MFloat = rad . realToFrac $ wrappedM+ oFloat = unsafeMapOrbit realToFrac o+ initialGuessFloat :: Angle Float+ Just initialGuessFloat = eccentricAnomalyAtMeanAnomalyFloat oFloat _MFloat+ initialGuess = realToFrac . (# [si|rad|]) $ initialGuessFloat+ err :: (Mode b, Floating b, Scalar b ~ a) => b -> b+ err _E = auto wrappedM - (_E - auto e * sin _E)+ wrappedE = fmap rad . convergeErr (runId . abs . err . Id) $+ findZeroNoEq err initialGuess+ _E = (|+| (unsafeMapUnit fromInteger n |*| turn)) <$> wrappedE++-- | 'eccentricAnomalyAtMeanAnomaly' specialized to 'Float'.+--+-- This function is used to calculate the initial guess for+-- 'eccentricAnomalyAtMeanAnomaly'.+eccentricAnomalyAtMeanAnomalyFloat :: Orbit Float -> Angle Float -> Maybe (Angle Float)+eccentricAnomalyAtMeanAnomalyFloat o _M = case classify o of+ Elliptic -> Just _E+ _ -> Nothing+ where wrappedM = (_M `mod'` turn) # [si|rad|]+ e = eccentricity o # [si||]+ sinM = sin wrappedM+ cosM = cos wrappedM+ -- Use a better initial guess+ -- http://alpheratz.net/dynamics/twobody/KeplerIterations_summary.pdf+ initialGuess = wrappedM ++ e * sinM ++ e * e * sinM * cosM ++ 0.5 * e * e * e * sinM * (3 * cosM * cosM - 1)+ _E :: Angle Float+ _E = rad . last . take 5 $+ findZero (\_E -> auto wrappedM - (_E - auto e * sin _E))+ initialGuess++-- | Calculate the eccentric anomaly, E, of an orbiting body when it has true+-- anomaly, ν.+--+-- The number of orbits represented by the anomalies is preserved;+-- i.e. ν `div` 2π = E `div` 2π+--+-- Returns Nothing if given a parabolic or hyperbolic orbit.+eccentricAnomalyAtTrueAnomaly :: (Floating a, Real a)+ => Orbit a -> Angle a -> Maybe (Angle a)+eccentricAnomalyAtTrueAnomaly o ν = case classify o of+ Elliptic -> Just _E+ _ -> Nothing+ where (n, wrappedν) = ν `divMod'` turn+ cosν = cos (ν # [si|rad|])+ -- sinν = sin (wrappedν # [si|rad|])+ e = eccentricity o # [si||]+ wrappedE = rad $ acos ((e + cosν) / (1 + e * cosν))+ -- wrappedE = rad $ atan2 (sqrt (1 - e*e) * sinν) (e + cosν)+ _E = if wrappedν < halfTurn+ then (unsafeMapUnit fromInteger n |*| turn) |+| wrappedE+ else (unsafeMapUnit fromInteger (n+1) |*| turn) |-| wrappedE++-- | 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++-- | 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++-- | Calculate the true anomaly, ν, of an orbiting body when it has the given+-- eccentric anomaly, _E.+--+-- The number of orbits represented by the anomalies is preserved;+-- i.e. ν `div` 2π = E `div` 2π+trueAnomalyAtEccentricAnomaly :: RealFloat a+ => Orbit a -- ^ An elliptic orbit+ -> Angle a -- ^ The eccentric anomaly _E+ -> Maybe (Angle a) -- ^ The true anomaly, ν+trueAnomalyAtEccentricAnomaly o _E = case classify o of+ Elliptic -> Just ν+ _ -> Nothing+ where (n, wrappedE) = bimap (unsafeMapUnit fromInteger) (# [si|rad|]) $+ _E `divMod'` turn+ e = eccentricity o # [si||]+ wrappedν = rad $ 2 * atan2 (sqrt (1 + e) * sin (wrappedE / 2))+ (sqrt (1 - e) * cos (wrappedE / 2))+ ν = turn |*| n |+| wrappedν++----------------------------------------------------------------+-- Utils+----------------------------------------------------------------++rad :: Fractional a => a -> Angle a+rad = (% [si|rad|])
+ test/Data/CReal/QuickCheck.hs view
@@ -0,0 +1,15 @@+{-# OPTIONS_GHC -fno-warn-orphans #-}++module Data.CReal.QuickCheck+ ( module Data.CReal+ ) where++import Data.CReal+import GHC.TypeLits+import Test.QuickCheck.Arbitrary (Arbitrary (..),+ arbitrarySizedFractional,+ shrinkRealFrac)++instance KnownNat n => Arbitrary (CReal n) where+ arbitrary = arbitrarySizedFractional+ shrink = shrinkRealFrac
+ test/Data/Metrology/QuickCheck.hs view
@@ -0,0 +1,31 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# OPTIONS_GHC -fno-warn-orphans #-}++module Data.Metrology.QuickCheck+ ( PositiveQuantity(..)+ ) where++import Data.Metrology.Unsafe+import System.Random ( Random )+import Test.QuickCheck ( Arbitrary(..)+ , Positive(..)+ )+import Test.QuickCheck.Checkers ( EqProp(..)+ , eq+ )++newtype PositiveQuantity a = PositiveQuantity { getPositiveQuantity :: a }++deriving instance Arbitrary a => Arbitrary (Qu u l a)++deriving instance Random a => Random (Qu u l a)++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++instance (Eq a) => EqProp (Qu u l a) where+ (=-=) = eq+
+ test/Physics/Orbit/QuickCheck.hs view
@@ -0,0 +1,154 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE QuasiQuotes #-}+{-# LANGUAGE RecordWildCards #-}+{-# OPTIONS_GHC -fno-warn-orphans #-}++module Physics.Orbit.QuickCheck+ ( CircularOrbit(..)+ , EllipticOrbit(..)+ , ParabolicOrbit(..)+ , HyperbolicOrbit(..)+ , unitOrbit+ ) where++import Data.Metrology+import Data.Metrology.Unsafe+import Data.Metrology.QuickCheck+import Data.Units.SI.Parser+import Physics.Orbit ( Distance+ , InclinationSpecifier(..)+ , Orbit(..)+ , PeriapsisSpecifier(..)+ , Unitless+ )+import System.Random ( Random )+import Test.QuickCheck ( Arbitrary(..)+ , choose+ , oneof+ , suchThat+ )++{-# ANN module ("HLint: ignore Reduce duplication" :: String) #-}++newtype CircularOrbit a = CircularOrbit {getCircularOrbit :: Orbit a}+ deriving(Show, Eq)++newtype EllipticOrbit a = EllipticOrbit {getEllipticOrbit :: Orbit a}+ deriving(Show, Eq)++newtype ParabolicOrbit a = ParabolicOrbit {getParabolicOrbit :: Orbit a}+ deriving(Show, Eq)++newtype HyperbolicOrbit a = HyperbolicOrbit {getHyperbolicOrbit :: Orbit a}+ deriving(Show, Eq)++-- | Use aerobreaking to shrink an orbit without expending fuel+instance (Num a, Ord a, Random a, Arbitrary a) => Arbitrary (Orbit a) where+ arbitrary = oneof+ [ getCircularOrbit <$> arbitrary+ , getEllipticOrbit <$> arbitrary+ , getParabolicOrbit <$> arbitrary+ , getHyperbolicOrbit <$> arbitrary+ ]+ shrink = shrinkOrbit++instance (Num a, Ord a, Arbitrary a) => Arbitrary (CircularOrbit a) where+ arbitrary =+ do+ let eccentricity = 0+ PositiveQuantity periapsis <- arbitrary+ inclinationSpecifier <- arbitrary+ let periapsisSpecifier = Circular+ PositiveQuantity primaryGravitationalParameter <- arbitrary+ pure . CircularOrbit $ Orbit { .. }+ shrink (CircularOrbit o) = CircularOrbit <$> shrinkOrbit o++instance (Num a, Ord a, Random a, Arbitrary a) => Arbitrary (EllipticOrbit a) where+ arbitrary =+ do+ eccentricity <- choose (0, 1) `suchThat` (/= 1)+ PositiveQuantity periapsis <- arbitrary+ inclinationSpecifier <- arbitrary+ periapsisSpecifier <- arbitrary+ PositiveQuantity primaryGravitationalParameter <- arbitrary+ pure . EllipticOrbit $ Orbit { .. }+ shrink (EllipticOrbit o) = EllipticOrbit <$> shrinkOrbit o++instance (Num a, Ord a, Random a, Arbitrary a) => Arbitrary (ParabolicOrbit a) where+ arbitrary =+ do+ let eccentricity = 1+ PositiveQuantity periapsis <- arbitrary+ inclinationSpecifier <- arbitrary+ periapsisSpecifier <- arbitrary+ PositiveQuantity primaryGravitationalParameter <- arbitrary+ pure . ParabolicOrbit $ Orbit { .. }+ shrink (ParabolicOrbit o) = ParabolicOrbit <$> shrinkOrbit o++instance (Num a, Ord a, Random a, Arbitrary a) => Arbitrary (HyperbolicOrbit a) where+ arbitrary =+ do+ eccentricity <- arbitrary `suchThat` (> 1)+ PositiveQuantity periapsis <- arbitrary+ inclinationSpecifier <- arbitrary+ periapsisSpecifier <- arbitrary+ PositiveQuantity primaryGravitationalParameter <- arbitrary+ pure . HyperbolicOrbit $ Orbit { .. }+ shrink (HyperbolicOrbit o) = HyperbolicOrbit <$> shrinkOrbit o++instance Arbitrary a => Arbitrary (InclinationSpecifier a) where+ arbitrary = oneof [pure NonInclined, Inclined <$> arbitrary <*> arbitrary]+ 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 Circular = []++--------------------------------------------------------------------------------+-- Shrinking+--------------------------------------------------------------------------------++-- | Note, this doesn't just lower the altitude, ho ho+shrinkOrbit :: (Arbitrary a, Num a, Ord a) => Orbit a -> [Orbit a]+shrinkOrbit o = [o{eccentricity = e} | e <- shrinkEccentricity (eccentricity o)] +++ [o{periapsis = q} | q <- shrinkPeriapsis (periapsis o)] +++ [o{inclinationSpecifier = i} | i <- shrink (inclinationSpecifier o)] +++ [o{periapsisSpecifier = ω} | ω <- shrink (periapsisSpecifier o)] +++ [o{primaryGravitationalParameter = μ} | μ <-+ shrinkPrimaryGravitationalParameter (primaryGravitationalParameter o)]++-- The semantics for shrinking lots of these values isn't to necessrily to+-- get a smaller value, but a more simple integral value could make+-- debugging easier. Try and skrink to the integers 0, 1, and 2+shrinkEccentricity :: (Num a, Ord a) => Unitless a -> [Unitless a]+shrinkEccentricity e | e == 0 || e == 1 || e == 2 = []+ | e < 1 = [0]+ | e > 1 = [2]+ | otherwise = error "shrinkEccentricity"++shrinkPeriapsis :: (Num a, Eq a) => Distance a -> [Distance a]+shrinkPeriapsis a | a == Qu 1 = []+ | otherwise = [Qu 1]++shrinkPrimaryGravitationalParameter+ :: (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]+++--------------------------------------------------------------------------------+-- Extras+--------------------------------------------------------------------------------++unitOrbit :: Fractional a => Orbit a+unitOrbit = Orbit{ eccentricity = 0+ , periapsis = 1 % [si|m|]+ , inclinationSpecifier = NonInclined+ , periapsisSpecifier = Circular+ , primaryGravitationalParameter = 1 % [si|m^3 s^-2|]+ }
+ test/Test.hs view
@@ -0,0 +1,436 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE QuasiQuotes #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TemplateHaskell #-}++module Main+ ( 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.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.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(..) )++{-# ANN module ("HLint: ignore Reduce duplication" :: String) #-}++-- | The type used for tests which require exact arithmetic. They are compared+-- at a resolution of 2^32+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+--------------------------------------------------------------------------------++test_sanity :: [TestTree]+test_sanity = [ testProperty "circular isValid"+ (\(CircularOrbit o) -> isValid (o :: Orbit Double))+ , testProperty "elliptic isValid"+ (\(EllipticOrbit o) -> isValid (o :: Orbit Double))+ , testProperty "parabolic isValid"+ (\(ParabolicOrbit o) -> isValid (o :: Orbit Double))+ , testProperty "hyperbolic isValid"+ (\(HyperbolicOrbit o) -> isValid (o :: Orbit Double))+ ]++test_classify :: [TestTree]+test_classify = [ testProperty "circular"+ (\(CircularOrbit o) ->+ classify (o :: Orbit Double) === Elliptic)+ , testProperty "elliptic"+ (\(EllipticOrbit o) ->+ classify (o :: Orbit Double) === Elliptic)+ , testProperty "parabolic"+ (\(ParabolicOrbit o) ->+ classify (o :: Orbit Double) === Parabolic)+ , testProperty "hyperbolic"+ (\(HyperbolicOrbit o) ->+ classify (o :: Orbit Double) === Hyperbolic)+ ]++test_semiMajorAxis :: [TestTree]+test_semiMajorAxis = [ testProperty "circular"+ (\(CircularOrbit o) ->+ fromJust (semiMajorAxis (o :: Orbit Double)) === periapsis o)+ , testProperty "elliptic"+ (\(EllipticOrbit o) ->+ fromJust (semiMajorAxis (o :: Orbit Double)) > zero)+ , testProperty "parabolic"+ (\(ParabolicOrbit o) ->+ semiMajorAxis (o :: Orbit Double) === Nothing)+ , testProperty "hyperbolic"+ (\(HyperbolicOrbit o) ->+ fromJust (semiMajorAxis (o :: Orbit Double)) < zero)+ ]++test_semiMinorAxis :: [TestTree]+test_semiMinorAxis = [ testGroup "range"+ [ testProperty "elliptic: b > 0"+ (\(EllipticOrbit o) ->+ semiMinorAxis (o :: Orbit Double) > zero)+ , testProperty "parabolic: b = 0"+ (\(ParabolicOrbit o) ->+ semiMinorAxis (o :: Orbit Double) === zero)+ , testProperty "hyperbolic: b < 0"+ (\(HyperbolicOrbit o) ->+ semiMinorAxis (o :: Orbit Double) < zero)+ ]+ , testProperty "semiMinorAxis circular = q"+ (\(CircularOrbit o) ->+ semiMinorAxis (o :: Orbit Double) === periapsis o)+ , testGroup "b^2 = a * l"+ [ testProperty "elliptic"+ (\(EllipticOrbit o) -> let a = fromJust (semiMajorAxis (o :: Orbit Exact))+ b = semiMinorAxis o+ l = semiLatusRectum o+ in b |*| b === a |*| l)+ , testProperty "hyperbolic"+ (\(HyperbolicOrbit o) -> let a = fromJust (semiMajorAxis (o :: Orbit Exact))+ b = semiMinorAxis o+ l = semiLatusRectum o+ in b |*| b === qNegate (a |*| l))+ ]+ ]++test_apoapsis :: [TestTree]+test_apoapsis = [ testProperty "ap > q"+ (\(EllipticOrbit o) ->+ eccentricity (o :: Orbit Double) /= 0+ ==> fromJust (apoapsis o) > periapsis o)+ , testProperty "circular: ap = q"+ (\(CircularOrbit o) ->+ fromJust (apoapsis (o :: Orbit Double)) === periapsis o)+ , testProperty "parabolic: no ap"+ (\(ParabolicOrbit o) ->+ apoapsis (o :: Orbit Double) === Nothing)+ , testProperty "hyperbolic: no ap"+ (\(HyperbolicOrbit o) ->+ apoapsis (o :: Orbit Double) === Nothing)+ ]++test_meanMotion :: [TestTree]+test_meanMotion = [ testProperty "n > 0"+ (\o -> meanMotion (o :: Orbit Double) > zero)+ ]++test_period :: [TestTree]+test_period = [ testProperty "p > 0"+ (\(EllipticOrbit o) ->+ fromJust (period (o :: Orbit Double)) > zero)+ , testProperty "4 π a^3 / p^2 = μ"+ (\(EllipticOrbit o) ->+ let Just p = period (o :: Orbit Exact)+ Just a = semiMajorAxis o+ μ = primaryGravitationalParameter o+ in (4 * qSq pi) |*| qCube a |/| qSq p === μ)+ , testProperty "parabolic: no p"+ (\(ParabolicOrbit o) ->+ period (o :: Orbit Double) === Nothing)+ , testProperty "hyperbolic: no p"+ (\(HyperbolicOrbit o) ->+ period (o :: Orbit Double) === Nothing)+ ]++-- TODO: Put converge test here+test_hyperbolicAngles :: [TestTree]+test_hyperbolicAngles = [ testProperty "parabolic approach"+ (\(ParabolicOrbit o) ->+ fromJust+ (hyperbolicApproachAngle (o :: Orbit Double))+ === qNegate halfTurn)+ , testProperty "parabolic departure"+ (\(ParabolicOrbit o) ->+ fromJust+ (hyperbolicDepartureAngle (o :: Orbit Double))+ === halfTurn)+ , testProperty "hyperbolic symmetry"+ (\(HyperbolicOrbit o) ->+ fromJust (hyperbolicDepartureAngle (o :: Orbit Double))+ === qNegate (fromJust (hyperbolicApproachAngle o)))+ , testProperty "elliptic: no approach"+ (\(EllipticOrbit o) ->+ hyperbolicApproachAngle (o :: Orbit Double) === Nothing)+ , testProperty "elliptic: no departure"+ (\(EllipticOrbit o) ->+ hyperbolicDepartureAngle (o :: Orbit Double) === Nothing)+ ]++anomalyConversionTests :: (forall a. (RealFloat a, Show a, Arbitrary a, Converge [a])+ => Orbit a -> Angle a -> Angle a)+ -> String -> String -> [TestTree]+anomalyConversionTests convertAnomaly fromName toName =+ [ testProperty (toName ++ " when " ++ fromName ++ " = 0")+ (\(EllipticOrbit o) ->+ let to = convertAnomaly (o :: Orbit Double) zero+ in to === zero)++ , testProperty (toName ++ " when " ++ fromName ++ " = π")+ (\(EllipticOrbit o) ->+ let to = convertAnomaly (o :: Orbit Double) halfTurn+ in to === halfTurn)++ , testProperty (toName ++ " when " ++ fromName ++ " = 2π")+ (\(EllipticOrbit o) ->+ let to = convertAnomaly (o :: Orbit Double) turn+ in to === turn)++ , testProperty "identity on circular orbits"+ (\(CircularOrbit o) from ->+ let to = convertAnomaly (o :: Orbit Exact) from+ in from === to)++ , testProperty "orbit number preservation"+ (\(EllipticOrbit o) from ->+ let to = convertAnomaly (o :: Orbit Double) from+ in from `div'` turn === (to `div'` turn :: Unitless Integer))+ ]++timeAnomalyConversionTests :: (forall a. (RealFloat a, Show a, Arbitrary a, Converge [a])+ => Orbit a -> Time a -> Angle a)+ -> String -> [TestTree]+timeAnomalyConversionTests timeToAnomaly toName =+ [ testProperty (toName ++ " when time = 0")+ (\(EllipticOrbit o) ->+ let to = timeToAnomaly (o :: Orbit Double) zero+ in to === zero)++ , testProperty (toName ++ " when time = p/2")+ (\(EllipticOrbit o) ->+ let to = timeToAnomaly (o :: Orbit Exact) (p|/|2)+ Just p = period o+ in to === halfTurn)++ , testProperty (toName ++ " when time = p")+ (\(EllipticOrbit o) ->+ let to = timeToAnomaly (o :: Orbit Exact) p+ Just p = period o+ in to === turn)++ , testProperty "identity on the unit orbit (modulo units!)"+ (\time ->+ let o = unitOrbit+ to = timeToAnomaly (o :: Orbit Exact) time+ in time # [si|s|] === to # [si|rad|])++ , testProperty "orbit number preservation"+ (\(EllipticOrbit o) time ->+ let to = timeToAnomaly (o :: Orbit Double) time+ Just p = period o+ in time `div'` p === (to `div'` turn :: Unitless Integer))+ ]++anomalyTimeConversionTests :: (forall a. (RealFloat a, Show a, Arbitrary a, Converge [a])+ => Orbit a -> Angle a -> Time a)+ -> String -> [TestTree]+anomalyTimeConversionTests anomalyToTime fromName =+ [ testProperty ("time when " ++ fromName ++ " = 0")+ (\(EllipticOrbit o) ->+ let t = anomalyToTime (o :: Orbit Double) zero+ in t === zero)++ , testProperty ("time when " ++ fromName ++ " = π")+ (\(EllipticOrbit o) ->+ let t = anomalyToTime (o :: Orbit Double) halfTurn+ Just p = period o+ in t === p |/| 2)++ , testProperty ("time when " ++ fromName ++ " = 2π")+ (\(EllipticOrbit o) ->+ let t = anomalyToTime (o :: Orbit Double) turn+ Just p = period o+ in t === p)++ , testProperty "identity on the unit orbit (modulo units!)"+ (\from ->+ let o = unitOrbit+ t = anomalyToTime (o :: Orbit Exact) from+ in from # [si|rad|] === t # [si|s|])++ , testProperty "orbit number preservation"+ (\(EllipticOrbit o) from ->+ let t = anomalyToTime (o :: Orbit Double) from+ Just p = period o+ in from `div'` turn === (t `div'` p :: Unitless Integer))+ ]++(.:) :: (a -> b) -> (c -> d -> a) -> c -> d -> b+f .: g = \x y -> f (g x y)++test_conversions :: [TestTree]+test_conversions = [ conversionToTime+ , conversionToMeanAnomaly+ , conversionToEccentricAnomaly+ , conversionToTrueAnomaly+ , conversionInverses+ ]+ where+ conversionToTime = testGroup "conversion to time"+ [ testGroup "from mean anomaly"+ (anomalyTimeConversionTests timeAtMeanAnomaly "mean anomaly")+ , testGroup "from eccentric anomaly"+ (anomalyTimeConversionTests (fromJust .: timeAtEccentricAnomaly)+ "eccentric anomaly")+ , testGroup "from true anomaly"+ (anomalyTimeConversionTests (fromJust .: timeAtTrueAnomaly)+ "true anomaly")+ ]++ conversionToMeanAnomaly = let s = "mean anomaly" in testGroup ("conversion to " ++ s)+ [ testGroup "from time"+ (timeAnomalyConversionTests meanAnomalyAtTime s)+ , testGroup "from eccentric anomaly"+ (anomalyConversionTests (fromJust .: meanAnomalyAtEccentricAnomaly)+ "eccentric anomaly"+ s)+ , testGroup "from true anomaly"+ (anomalyConversionTests (fromJust .: meanAnomalyAtTrueAnomaly)+ "true anomaly"+ s)+ ]++ conversionToEccentricAnomaly = let s = "eccentric anomaly" in testGroup ("conversion to " ++ s)+ [ testGroup "from time"+ (timeAnomalyConversionTests (fromJust .: eccentricAnomalyAtTime) s)+ , testGroup "from mean anomaly"+ (anomalyConversionTests (fromJust .: eccentricAnomalyAtMeanAnomaly)+ "mean anomaly"+ s)+ , testGroup "from true anomaly"+ (anomalyConversionTests (fromJust .: eccentricAnomalyAtTrueAnomaly)+ "true anomaly"+ s)+ ]++ conversionToTrueAnomaly = let s = "true anomaly" in testGroup ("conversion to " ++ s)+ [ testGroup "from time"+ (timeAnomalyConversionTests (fromJust .: trueAnomalyAtTime) s)+ , testGroup "from mean anomaly"+ (anomalyConversionTests (fromJust .: trueAnomalyAtMeanAnomaly)+ "mean anomaly"+ s)+ , testGroup "from eccentric anomaly"+ (anomalyConversionTests (fromJust .: trueAnomalyAtEccentricAnomaly)+ "eccentric anomaly"+ s)+ ]++ conversionInverses = testGroup "conversionInverses"+ [ testProperty "mean time inverse"+ (\o -> inverse (meanAnomalyAtTime (o :: Orbit Exact))+ (timeAtMeanAnomaly o))++ , slowTest $ testProperty "mean eccentric inverse"+ (\(EllipticOrbit o) ->+ inverse (coerce (fromJust . meanAnomalyAtEccentricAnomaly (o :: Orbit Exact)) :: WrappedAngle Exact -> WrappedAngle Exact)+ (coerce (fromJust . eccentricAnomalyAtMeanAnomaly o)))++ , slowTest $ testProperty "mean true inverse"+ (\(EllipticOrbit o) ->+ inverse (fromJust . meanAnomalyAtTrueAnomaly (o :: Orbit Exact))+ (fromJust . trueAnomalyAtMeanAnomaly o))++ , slowTest $ testProperty "time true inverse"+ (\(EllipticOrbit o) ->+ inverse (fromJust . timeAtTrueAnomaly (o :: Orbit Exact))+ (fromJust . trueAnomalyAtTime o))++ , testProperty "time eccentric inverse"+ (\(EllipticOrbit o) ->+ inverse (fromJust . timeAtEccentricAnomaly (o :: Orbit Exact))+ (fromJust . eccentricAnomalyAtTime o))++ , testProperty "eccentric true inverse"+ (\(EllipticOrbit o) ->+ inverse (coerce (fromJust . eccentricAnomalyAtTrueAnomaly (o:: Orbit Exact)) :: WrappedAngle Exact -> WrappedAngle Exact)+ (fromJust . coerce (trueAnomalyAtEccentricAnomaly o)))+ ]++-- TODO: Put parabolic and hyperbolic tests here+test_areal :: [TestTree]+test_areal = [ testProperty "elliptic areal area"+ (\(EllipticOrbit o) -> let Just a = semiMajorAxis (o :: Orbit Exact)+ b = semiMinorAxis o+ area = pi |*| a |*| b+ Just p = period o+ in area === p |*| arealVelocity o)+ ]++main :: IO ()+main = do+ let is = includingOptions [slowTestQCRatio] : defaultIngredients+ defaultMainWithIngredients is $(testGroupGenerator)+
+ test/Test/QuickCheck/Extra.hs view
@@ -0,0 +1,15 @@+module Test.QuickCheck.Extra+ ( (<=!)+ , (>=!)+ ) where++import Test.QuickCheck (Property, counterexample)++infix 4 <=!+(<=!) :: (Ord a, Show a) => a -> a -> Property+x <=! y = counterexample (show x ++ " ≰ " ++ show y) (x <= y)++infix 4 >=!+(>=!) :: (Ord a, Show a) => a -> a -> Property+x >=! y = counterexample (show x ++ " ≱ " ++ show y) (x >= y)+
+ test/WrappedAngle.hs view
@@ -0,0 +1,24 @@+module WrappedAngle+ ( WrappedAngle(..)+ ) where++import Data.Constants.Mechanics.Extra+import Data.Metrology+import Data.Metrology.Extra+import Data.Metrology.QuickCheck ( )+import Data.Metrology.Show ( )+import Data.Units.SI.Parser+import Test.QuickCheck.Checkers ( EqProp(..)+ , eq+ )+import Test.Tasty.QuickCheck ( Arbitrary )++-- A wrapper which compares angles for equality modulo 2π+newtype WrappedAngle a = WrappedAngle (MkQu_ULN [si|rad|] 'DefaultLCSU a)+ deriving (Show, Arbitrary)++instance (Floating a, Real a) => Eq (WrappedAngle a) where+ WrappedAngle x == WrappedAngle y = (x `mod'` turn) == (y `mod'` turn)++instance (Floating a, Real a) => EqProp (WrappedAngle a) where+ (=-=) = eq
+ test/doctest/Doctests.hs view
@@ -0,0 +1,10 @@+module Main where++import Build_doctests ( flags+ , module_sources+ , pkgs+ )+import Test.DocTest ( doctest )++main :: IO ()+main = doctest $ flags ++ pkgs ++ module_sources