Physics.hs

{-# LANGUAGE
	UndecidableInstances,
	FlexibleInstances,
	MultiParamTypeClasses,
	EmptyDataDecls,
	TypeFamilies #-}
module Physics (Physicsable) where

import Data.Maybe


class NonNegative a
class NonPositive a

data Zero
instance NonNegative Zero
instance NonPositive Zero

data (NonNegative a) => Succ a = Succ a
data (NonPositive a) => Prev a = Prev a

instance (NonNegative a) => NonNegative (Succ a)
instance (NonPositive a) => NonPositive (Prev a)

class Add a b where
  type Sum a b
  (<+>) :: a -> b -> Sum a b
  _ <+> _ = undefined
  type Diff a b
  (<->) :: a -> b -> Diff a b
  _ <-> _ = undefined
  
class Negate a where
  type Neg a
  neg :: a -> Neg a
  neg _ = undefined

instance Negate Zero where type Neg Zero = Zero
instance Negate (Succ a) where type Neg (Succ a) = Prev (Neg a)
instance Negate (Prev a) where type Neg (Prev a) = Succ (Neg a)

instance Add Zero b where
  type Sum Zero b = b
  type Diff Zero b = Sum Zero (Neg b)
instance Add a Zero where
  type Sum a Zero = a
  type Diff a Zero = Sum a (Neg Zero)

instance Add (Succ a) (Succ b) where
  type Sum (Succ a) (Succ b) = Sum (Succ (Succ a)) b
  type Diff (Succ a) (Succ b) = Sum (Succ a) (Neg (Succ b))

instance Add (Prev a) (Prev b) where
  type Sum (Prev a) (Prev b) = Sum (Prev (Prev a)) b
  type Diff (Prev a) (Prev b) = Sum (Prev a) (Neg (Prev b))

instance Add (Succ a) (Prev b) where
  type Sum (Succ a) (Prev b) = Sum a b
  type Diff (Succ a) (Prev b) = Sum (Succ a) (Neg (Prev b))
instance Add (Prev a) (Succ b) where
  type Sum (Prev a) (Succ b) = Sum a b
  type Diff (Prev a) (Succ b) = Sum (Prev a) (Neg (Succ b))

data Physical s m g = Physical
			{ seconds :: s,
			  meters :: m,
			  kgrams :: g,
			  value :: Double }

instance Eq (Physical s m g) where
  a == b = (value a) == (value b)
instance Show (Physical s m g) where
  show = show . value

instance Add (Physical s m g) (Physical s m g) where
  type Sum (Physical s m g) (Physical s m g) = Physical s m g
  a <+> b = Physical
		{ seconds = undefined,
		  meters  = undefined,
		  kgrams  = undefined,
		  value = (value a) + (value b) }
  type Diff (Physical s m g) (Physical s m g) = Physical s m g
  a <-> b = Physical
		{ seconds = undefined,
		  meters  = undefined,
		  kgrams  = undefined,
		  value = (value a) - (value b) }

type KGrams = Physical Zero Zero (Succ Zero)
type Seconds = Physical (Succ Zero) Zero Zero
type Meters  = Physical Zero (Succ Zero) Zero
type MetersPerSecond = Div Meters Seconds
type MetersPerSecond2 = Div (Div Meters Seconds) Seconds
type Newtons = Prod KGrams MetersPerSecond2

class Mult a b where
  type Prod a b
  (<*>) :: a -> b -> Prod a b
  _ <*> _ = undefined
  type Div a b
  (</>) :: a -> b -> Div a b
  _ </> _ = undefined

class Invert a where
  type Inv a
  inv :: a -> Inv a
  inv _ = undefined

instance Invert (Physical s m g) where
  type Inv (Physical s m g) = Physical (Neg s) (Neg m) (Neg g)
  inv p = Physical
		{ seconds = undefined,
		  meters  = undefined,
		  kgrams  = undefined,
		  value = 1/(value p) }

instance Mult (Physical s m g) (Physical s' m' g') where
  type Prod (Physical s m g) (Physical s' m' g') =
	Physical (Sum s s') (Sum m m') (Sum g g')
  a <*> b = Physical
		{ seconds = undefined,
		  meters  = undefined,
		  kgrams  = undefined,
		  value = (value a) * (value b) }
  type Div (Physical s m g) (Physical s' m' g') =
	Physical (Diff s s') (Diff m m') (Diff g g')
  a </> b = Physical
		{ seconds = undefined,
		  meters  = undefined,
		  kgrams  = undefined,
		  value = (value a) / (value b) }

instance (Add a a', Add b b', Add c c') => Add (a, b, c) (a', b', c') where
  type Sum (a, b, c) (a', b', c') = (Sum a a', Sum b b', Sum c c')
  (a, b, c) <+> (a', b', c') = (a <+> a', b <+> b', c <+> c')
  type Diff (a, b, c) (a', b', c') = (Diff a a', Diff b b', Diff c c')
  (a, b, c) <-> (a', b', c') = (a <-> a', b <-> b', c <-> c')

instance (Mult s a, Mult s b, Mult s c) => Mult s (a, b, c) where
  type Prod s (a, b, c) = (Prod s a, Prod s b, Prod s c)
  s <*> (a, b, c) = (s <*> a, s <*> b, s <*> c)
  type Div s (a, b, c) = (Div s a, Div s b, Div s c)
  s </> (a, b, c) = (s </> a, s </> b, s </> c)

type Position = (Meters, Meters, Meters)
type Velocity = (MetersPerSecond, MetersPerSecond, MetersPerSecond)
type Acceleration = (MetersPerSecond2, MetersPerSecond2, MetersPerSecond2)
type Force = (Newtons, Newtons, Newtons)
type Mass = KGrams
type Time = Seconds

class Physicsable a where
  mass :: a -> Mass
  position :: a -> Position
  velocity :: a -> Velocity
  acceleration :: a -> Acceleration
  force :: a -> Force
  updatePhysics :: a -> (Position, Velocity, Acceleration) -> a
  doPhysics :: a -> Time -> a
  doPhysics a t = updatePhysics a (position', velocity', acceleration')
    where iterations = round $ toScalar (t </> h)
	  h = (physical 0.1 :: Time)
	  physics' = take iterations $ iterate (newton h)
		   $ (physical 0 :: Time, position a, velocity a, acceleration')
	  position' = (\(t1, p1, v1, a1) -> p1) $ last physics'
	  velocity' = (\(t1, p1, v1, a1) -> v1) $ last physics'
	  acceleration' = (inv $ mass a) <*> (force a)

newton :: Time -> (Time, Position, Velocity, Acceleration) -> (Time, Position, Velocity, Acceleration)
newton h (t0, p0, v0, a0) = (t1, p1, v1, a0)
  where (t1, v1) = rk4 (const $ const a0) v0 t0 h
	(_,  p1) = rk4 (const $ const v1) p0 t0 h

toScalar :: Physical Zero Zero Zero -> Double
toScalar = value

scalar :: Double -> Physical Zero Zero Zero
scalar = physical

physical :: Double -> Physical s m g
physical d = Physical
		{ seconds = undefined,
		  meters  = undefined,
		  kgrams  = undefined,
		  value = d }

-- Arguments:
-- rk4 f(t,y) y0 t0 h
-- idea: y' = f(t,y) and y(t0) = y0
-- computes (t1 = t0 + h) and y1
rk4 f y0 t0 h = (t0 <+> h,
		 y0 <+> ((scalar (1/6)) <*> h <*>
			 (k1
			  <+> ((scalar 2)<*>k2)
			  <+> ((scalar 2)<*>k3)
			  <+> k4)))
  where k1 = f t0 y0
	k2 = f (t0 <+> ((scalar (1/2)) <*> h))
	       (y0 <+> ((scalar (1/2)) <*> h <*> k1))
	k3 = f (t0 <+> ((scalar (1/2)) <*> h))
	       (y0 <+> ((scalar (1/2)) <*> h <*> k2))
	k4 = f (t0 <+> h) (y0 <+> (h <*> k3))


{-
Now for a quick test of the system.  We take constant force F = (1,2,3), mass m = 1,
initial values all zero.  Run it for 5 seconds.
-}


data DeadWeight = DW { m :: Mass, p :: Position, v :: Velocity, a :: Acceleration, f :: Force } deriving Show
instance Physicsable DeadWeight where
  mass = m
  position = p
  velocity = v
  acceleration = a
  force = f
  updatePhysics dw (p', v', a') = dw{ p=p', v=v', a=a' }

dw = DW {
	m = physical 1,
	p = (physical 0, physical 0, physical 0),
	v = (physical 0, physical 0, physical 0),
	a = (physical 0, physical 0, physical 0),
	f = (physical 1, physical 2, physical 3)
	}

dw' = doPhysics dw (physical 5 :: Time)