orbiting.hs

{-# LANGUAGE FlexibleInstances,
             MultiParamTypeClasses,
             ScopedTypeVariables,
             TypeFamilies,
             TypeOperators,
             UndecidableInstances #-}
{-# OPTIONS_GHC -funbox-strict-fields #-}

module Main where


------------------------------------------------------------------------------
-- Type operators

infixr 8 :^:
infixl 7 :*:, :/:
infixl 6 :+:, :-:


type family a :+: b
type family Negate a
type family a :-: b
type family a :*: b
type family a :/: b
type family a :^: b


------------------------------------------------------------------------------
-- Type naturals

data Z
data S a


type instance Z :+: Z = Z
type instance Z :+: S a = S a
type instance S a :+: Z = S a
type instance S a :+: S b = S (S (a :+: b))

type instance Z :-: Z = Z
type instance S a :-: Z = S a
type instance S a :-: S b = a :-: b
-- No type instance Z :-: S a
-- Reducible only when the minuend is greater than the subtrahend

type instance Z :*: Z = Z
type instance Z :*: S a = Z
type instance S a :*: Z = Z
type instance S a :*: S b = S a :*: b :+: S a

type instance Z :/: S a = Z
type instance S a :/: S b = DivAux (S a) (S b) Z
-- No type instance Z :/: Z
-- No type instance S a :/: Z
-- Reducible only when the dividend is an integral multiple of the non-zero divisor

type instance Z :^: Z = S Z
type instance Z :^: S a = Z
type instance S a :^: Z = S Z
type instance S a :^: S b = S a :^: b :*: S a


type family DivAux a b x
type instance DivAux Z b x = x
type instance DivAux (S a) b x = DivAux (S a :-: b) b (S x)


class TypeNat a where
  reflectTypeNat :: (Num n) => a -> n

instance TypeNat Z where
  reflectTypeNat _ = 0
instance (TypeNat a) => TypeNat (S a) where
  reflectTypeNat _ = reflectTypeNat (undefined :: a) + 1


type ZeroNat = Z
type OneNat = S Z
type TwoNat = S (S Z)
type ThreeNat = S (S (S Z))


------------------------------------------------------------------------------
-- Type integers

data P a b


type instance P a b :+: P c d = Norm (P (a :+: c) (b :+: d))

type instance Negate (P a b) = P b a

type instance P a b :-: P c d = P a b :+: Negate (P c d)

type instance P a b :*: P c d = Norm (P ((a :*: c) :+: (b :*: d)) ((a :*: d) :+: (b :*: c)))

type instance P Z Z :/: P Z (S a) = P Z Z
type instance P Z Z :/: P (S a) Z = P Z Z
type instance P Z (S a) :/: P Z (S b) = P (S a :/: S b) Z
type instance P Z (S a) :/: P (S b) Z = P Z (S a :/: S b)
type instance P (S a) Z :/: P Z (S b) = P Z (S a :/: S b)
type instance P (S a) Z :/: P (S b) Z = P (S a :/: S b) Z
-- No type instance P Z Z :/: P Z Z
-- No type instance P Z (S a) :/: P Z Z
-- No type instance P (S a) Z :/: P Z Z
-- Reducible only when the dividend is an integral multiple of the non-zero divisor


type family Norm a
type instance Norm (P Z Z) = P Z Z
type instance Norm (P Z (S a)) = P Z (S a)
type instance Norm (P (S a) Z) = P (S a) Z
type instance Norm (P (S a) (S b)) = Norm (P a b)


class TypeInt a where
  reflectTypeInt :: (Num n) => a -> n

instance (TypeNat a, TypeNat b) => TypeInt (P a b) where
  reflectTypeInt _ = reflectTypeNat (undefined :: a) - reflectTypeNat (undefined :: b)


type Zero = P Z Z
type One = P (S Z) Z
type Two = P (S (S Z)) Z
type Three = P (S (S (S Z))) Z
type MinusOne = P Z (S Z)
type MinusTwo = P Z (S (S Z))
type MinusThree = P Z (S (S (S Z)))


------------------------------------------------------------------------------
-- Unit-aware operators

infixl 7 ***
infixl 7 ///


class Mulable a b where
  (***) :: a -> b -> a :*: b

class Divable a b where
  (///) :: a -> b -> a :/: b


------------------------------------------------------------------------------
-- Numbers

type instance Double :+: Double = Double
type instance Negate Double = Double
type instance Double :-: Double = Double
type instance Double :*: Double = Double
type instance Double :/: Double = Double
type instance Double :^: Double = Double

instance Mulable Double Double where
  (***) = (*)

instance Divable Double Double where
  (///) = (/)


------------------------------------------------------------------------------
-- Quantities

data Quantity m kg s = Quantity !Double
  deriving (Eq, Ord)


type instance Quantity m kg s :*: Quantity m' kg' s' = Quantity (m :+: m') (kg :+: kg') (s :+: s')

type instance Quantity m kg s :/: Quantity m' kg' s' = Quantity (m :-: m') (kg :-: kg') (s :-: s')

type instance Quantity m kg s :^: P Z Z = Quantity Zero Zero Zero
type instance Quantity m kg s :^: P (S a) Z = Quantity m kg s :^: P a Z :*: Quantity m kg s
type instance Quantity m kg s :^: P Z (S a) = Quantity m kg s :^: P Z a :/: Quantity m kg s


class TypeQuantity a where
  reflectTypeQuantity :: (Num n) => a -> (n, n, n)

instance (TypeInt m, TypeInt kg, TypeInt s, a ~ Quantity m kg s) => TypeQuantity a where
  reflectTypeQuantity _ = (reflectTypeInt (undefined :: m), reflectTypeInt (undefined :: kg), reflectTypeInt (undefined :: s))


instance (TypeInt m, TypeInt kg, TypeInt s) => Show (Quantity m kg s) where
  show q@(Quantity a) = show a ++ " m^" ++ show m ++ " kg^" ++ show kg ++ " s^" ++ show s
    where
      (m, kg, s) = reflectTypeQuantity q


instance Mulable (Quantity m kg s) (Quantity m' kg' s') where
  Quantity a *** Quantity b = Quantity (a * b)

instance Divable (Quantity m kg s) (Quantity m' kg' s') where
  Quantity a /// Quantity b = Quantity (a / b)


unq :: Quantity m kg s -> Double
unq (Quantity a) = a

qlift1 :: (Double -> Double) -> Quantity m kg s -> Quantity m kg s
qlift1 f (Quantity a) = Quantity (f a)

qlift2 :: (Double -> Double -> Double) -> Quantity m kg s -> Quantity m kg s -> Quantity m kg s
qlift2 f (Quantity a) (Quantity b) = Quantity (f a b)


instance Num (Quantity m kg s) where
  (+) = qlift2 (+)
  (-) = qlift2 (-)
  negate = qlift1 negate
  (*) = undefined
  abs = qlift1 abs
  signum = qlift1 signum
  fromInteger = Quantity . fromInteger

instance Real (Quantity m kg s) where
  toRational = toRational . unq

instance Fractional (Quantity m kg s) where
  (/) = undefined
  fromRational = Quantity . fromRational


type Meter = Quantity One Zero Zero
type Kilogram = Quantity Zero One Zero
type Second = Quantity Zero Zero One
type Newton = Quantity One One MinusTwo


------------------------------------------------------------------------------
-- Physics

gravF :: Kilogram -> Kilogram -> Meter -> Newton
gravF m1 m2 r = gravC *** m1 *** m2 /// (r *** r)

gravC :: Newton :*: Meter:^:Two :/: Kilogram:^:Two
gravC = 6.67428e-11

earthM :: Kilogram
earthM = 6.0e24

earthR :: Meter
earthR = 6.357e6