diff --git a/.gitignore b/.gitignore new file mode 100644 --- /dev/null +++ b/.gitignore @@ -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 diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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 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 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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,6 @@ +module Main where + +import Distribution.Extra.Doctest (defaultMainWithDoctests) + +main :: IO () +main = defaultMainWithDoctests "doctests" diff --git a/changelog.md b/changelog.md new file mode 100644 --- /dev/null +++ b/changelog.md @@ -0,0 +1,6 @@ +# Change Log + +## WIP + +## 0.3 + - Switch to `units` from `uom-plugin` diff --git a/default.nix b/default.nix new file mode 100644 --- /dev/null +++ b/default.nix @@ -0,0 +1,45 @@ +{ pkgs ? import { }, 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 diff --git a/orbits.cabal b/orbits.cabal new file mode 100644 --- /dev/null +++ b/orbits.cabal @@ -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 +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 diff --git a/readme.md b/readme.md new file mode 100644 --- /dev/null +++ b/readme.md @@ -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 + diff --git a/src/Data/Constants/Mechanics/Extra.hs b/src/Data/Constants/Mechanics/Extra.hs new file mode 100644 --- /dev/null +++ b/src/Data/Constants/Mechanics/Extra.hs @@ -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 diff --git a/src/Data/Metrology/Extra.hs b/src/Data/Metrology/Extra.hs new file mode 100644 --- /dev/null +++ b/src/Data/Metrology/Extra.hs @@ -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)) diff --git a/src/Physics/Orbit.hs b/src/Physics/Orbit.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Orbit.hs @@ -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 + -- . + -- + -- '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 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 +-- , 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 , +-- 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|]) diff --git a/test/Data/CReal/QuickCheck.hs b/test/Data/CReal/QuickCheck.hs new file mode 100644 --- /dev/null +++ b/test/Data/CReal/QuickCheck.hs @@ -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 diff --git a/test/Data/Metrology/QuickCheck.hs b/test/Data/Metrology/QuickCheck.hs new file mode 100644 --- /dev/null +++ b/test/Data/Metrology/QuickCheck.hs @@ -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 + diff --git a/test/Physics/Orbit/QuickCheck.hs b/test/Physics/Orbit/QuickCheck.hs new file mode 100644 --- /dev/null +++ b/test/Physics/Orbit/QuickCheck.hs @@ -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|] + } diff --git a/test/Test.hs b/test/Test.hs new file mode 100644 --- /dev/null +++ b/test/Test.hs @@ -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) + diff --git a/test/Test/QuickCheck/Extra.hs b/test/Test/QuickCheck/Extra.hs new file mode 100644 --- /dev/null +++ b/test/Test/QuickCheck/Extra.hs @@ -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) + diff --git a/test/WrappedAngle.hs b/test/WrappedAngle.hs new file mode 100644 --- /dev/null +++ b/test/WrappedAngle.hs @@ -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 diff --git a/test/doctest/Doctests.hs b/test/doctest/Doctests.hs new file mode 100644 --- /dev/null +++ b/test/doctest/Doctests.hs @@ -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