{-# LANGUAGE GeneralizedNewtypeDeriving #-}
module Physics where

import Control.Applicative (liftA2)

import Geometry
import Numerics.VectorAnalysis

data PhaseSpace = PhaseSpace
    { PhaseSpace -> Vec2
p :: Vec2
    , PhaseSpace -> Vec2
q :: Vec2
    } deriving (PhaseSpace -> PhaseSpace -> Bool
(PhaseSpace -> PhaseSpace -> Bool)
-> (PhaseSpace -> PhaseSpace -> Bool) -> Eq PhaseSpace
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: PhaseSpace -> PhaseSpace -> Bool
== :: PhaseSpace -> PhaseSpace -> Bool
$c/= :: PhaseSpace -> PhaseSpace -> Bool
/= :: PhaseSpace -> PhaseSpace -> Bool
Eq, Int -> PhaseSpace -> ShowS
[PhaseSpace] -> ShowS
PhaseSpace -> String
(Int -> PhaseSpace -> ShowS)
-> (PhaseSpace -> String)
-> ([PhaseSpace] -> ShowS)
-> Show PhaseSpace
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> PhaseSpace -> ShowS
showsPrec :: Int -> PhaseSpace -> ShowS
$cshow :: PhaseSpace -> String
show :: PhaseSpace -> String
$cshowList :: [PhaseSpace] -> ShowS
showList :: [PhaseSpace] -> ShowS
Show)

instance VectorSpace PhaseSpace where
    zero :: PhaseSpace
zero = Vec2 -> Vec2 -> PhaseSpace
PhaseSpace Vec2
forall v. VectorSpace v => v
zero Vec2
forall v. VectorSpace v => v
zero
    PhaseSpace Vec2
p1 Vec2
q1 +. :: PhaseSpace -> PhaseSpace -> PhaseSpace
+. PhaseSpace Vec2
p2 Vec2
q2 = Vec2 -> Vec2 -> PhaseSpace
PhaseSpace (Vec2
p1 Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Vec2
p2) (Vec2
q1 Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Vec2
q2)
    Double
c *. :: Double -> PhaseSpace -> PhaseSpace
*. PhaseSpace Vec2
p Vec2
q = Vec2 -> Vec2 -> PhaseSpace
PhaseSpace (Double
c Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
p) (Double
c Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
q)
    negateV :: PhaseSpace -> PhaseSpace
negateV (PhaseSpace Vec2
p Vec2
q) = Vec2 -> Vec2 -> PhaseSpace
PhaseSpace (Vec2 -> Vec2
forall v. VectorSpace v => v -> v
negateV Vec2
p) (Vec2 -> Vec2
forall v. VectorSpace v => v -> v
negateV Vec2
q)

particleInPotential :: Double -> (Vec2 -> Double) -> PhaseSpace -> PhaseSpace
particleInPotential :: Double -> (Vec2 -> Double) -> PhaseSpace -> PhaseSpace
particleInPotential Double
mass Vec2 -> Double
potential PhaseSpace{Vec2
p :: PhaseSpace -> Vec2
q :: PhaseSpace -> Vec2
p :: Vec2
q :: Vec2
..} = Vec2 -> Vec2 -> PhaseSpace
PhaseSpace Vec2
deltaP Vec2
deltaQ
  where
    deltaP :: Vec2
deltaP = Vec2 -> Vec2
forall v. VectorSpace v => v -> v
negateV ((Vec2 -> Double) -> Vec2 -> Vec2
grad Vec2 -> Double
potential Vec2
q)
    deltaQ :: Vec2
deltaQ = Vec2
p Vec2 -> Double -> Vec2
forall v. VectorSpace v => v -> Double -> v
/. Double
mass

--particleInPotentialWithFriction :: Double -> (Vec2 -> Double) -> Double -> PhaseSpace -> PhaseSpace
--particleInPotentialWithFriction mass potential friction = PhaseSpace deltaP deltaQ
--  where
--    deltaP = negateV (grad potential q)
--    deltaQ = p /. mass


harmonicPotential :: (Double, Double) -> Vec2 -> Vec2 -> Double
harmonicPotential :: (Double, Double) -> Vec2 -> Vec2 -> Double
harmonicPotential (Double
w, Double
h) (Vec2 Double
x0 Double
y0) (Vec2 Double
x Double
y) = ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
x0)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
w)Double -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
y0)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
h)Double -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2

coulombPotential :: Double -> Vec2 -> Vec2 -> Double
coulombPotential :: Double -> Vec2 -> Vec2 -> Double
coulombPotential Double
charge Vec2
center Vec2
particle = Double
charge Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Vec2 -> Double
norm (Vec2
particle Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
center)

twoBody :: (Vec2 -> Double) -> (Vec2 -> Vec2 -> Double) -> (Double, Double) -> (PhaseSpace, PhaseSpace) -> (PhaseSpace, PhaseSpace)
twoBody :: (Vec2 -> Double)
-> (Vec2 -> Vec2 -> Double)
-> (Double, Double)
-> (PhaseSpace, PhaseSpace)
-> (PhaseSpace, PhaseSpace)
twoBody Vec2 -> Double
potential Vec2 -> Vec2 -> Double
interaction (Double
m1, Double
m2) (PhaseSpace Vec2
p1 Vec2
q1, PhaseSpace Vec2
p2 Vec2
q2) = (Vec2 -> Vec2 -> PhaseSpace
PhaseSpace Vec2
deltaP1 Vec2
deltaQ1, Vec2 -> Vec2 -> PhaseSpace
PhaseSpace Vec2
deltaP2 Vec2
deltaQ2)
  where
    deltaP1 :: Vec2
deltaP1 = Vec2 -> Vec2
forall v. VectorSpace v => v -> v
negateV ((Vec2 -> Double) -> Vec2 -> Vec2
grad (Vec2 -> Double
potential (Vec2 -> Double) -> (Vec2 -> Double) -> Vec2 -> Double
forall v. VectorSpace v => v -> v -> v
+. Vec2 -> Vec2 -> Double
interaction Vec2
q2) Vec2
q1)
    deltaP2 :: Vec2
deltaP2 = Vec2 -> Vec2
forall v. VectorSpace v => v -> v
negateV ((Vec2 -> Double) -> Vec2 -> Vec2
grad (Vec2 -> Double
potential (Vec2 -> Double) -> (Vec2 -> Double) -> Vec2 -> Double
forall v. VectorSpace v => v -> v -> v
+. Vec2 -> Vec2 -> Double
interaction Vec2
q1) Vec2
q2)
    deltaQ1 :: Vec2
deltaQ1 = Vec2
p1 Vec2 -> Double -> Vec2
forall v. VectorSpace v => v -> Double -> v
/. Double
m1
    deltaQ2 :: Vec2
deltaQ2 = Vec2
p2 Vec2 -> Double -> Vec2
forall v. VectorSpace v => v -> Double -> v
/. Double
m2

newtype NBody a = NBody { forall a. NBody a -> [a]
getNBody :: [a] } deriving (NBody a -> NBody a -> Bool
(NBody a -> NBody a -> Bool)
-> (NBody a -> NBody a -> Bool) -> Eq (NBody a)
forall a. Eq a => NBody a -> NBody a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => NBody a -> NBody a -> Bool
== :: NBody a -> NBody a -> Bool
$c/= :: forall a. Eq a => NBody a -> NBody a -> Bool
/= :: NBody a -> NBody a -> Bool
Eq, Int -> NBody a -> ShowS
[NBody a] -> ShowS
NBody a -> String
(Int -> NBody a -> ShowS)
-> (NBody a -> String) -> ([NBody a] -> ShowS) -> Show (NBody a)
forall a. Show a => Int -> NBody a -> ShowS
forall a. Show a => [NBody a] -> ShowS
forall a. Show a => NBody a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> NBody a -> ShowS
showsPrec :: Int -> NBody a -> ShowS
$cshow :: forall a. Show a => NBody a -> String
show :: NBody a -> String
$cshowList :: forall a. Show a => [NBody a] -> ShowS
showList :: [NBody a] -> ShowS
Show, (forall a b. (a -> b) -> NBody a -> NBody b)
-> (forall a b. a -> NBody b -> NBody a) -> Functor NBody
forall a b. a -> NBody b -> NBody a
forall a b. (a -> b) -> NBody a -> NBody b
forall (f :: * -> *).
(forall a b. (a -> b) -> f a -> f b)
-> (forall a b. a -> f b -> f a) -> Functor f
$cfmap :: forall a b. (a -> b) -> NBody a -> NBody b
fmap :: forall a b. (a -> b) -> NBody a -> NBody b
$c<$ :: forall a b. a -> NBody b -> NBody a
<$ :: forall a b. a -> NBody b -> NBody a
Functor, (forall m. Monoid m => NBody m -> m)
-> (forall m a. Monoid m => (a -> m) -> NBody a -> m)
-> (forall m a. Monoid m => (a -> m) -> NBody a -> m)
-> (forall a b. (a -> b -> b) -> b -> NBody a -> b)
-> (forall a b. (a -> b -> b) -> b -> NBody a -> b)
-> (forall b a. (b -> a -> b) -> b -> NBody a -> b)
-> (forall b a. (b -> a -> b) -> b -> NBody a -> b)
-> (forall a. (a -> a -> a) -> NBody a -> a)
-> (forall a. (a -> a -> a) -> NBody a -> a)
-> (forall a. NBody a -> [a])
-> (forall a. NBody a -> Bool)
-> (forall a. NBody a -> Int)
-> (forall a. Eq a => a -> NBody a -> Bool)
-> (forall a. Ord a => NBody a -> a)
-> (forall a. Ord a => NBody a -> a)
-> (forall a. Num a => NBody a -> a)
-> (forall a. Num a => NBody a -> a)
-> Foldable NBody
forall a. Eq a => a -> NBody a -> Bool
forall a. Num a => NBody a -> a
forall a. Ord a => NBody a -> a
forall m. Monoid m => NBody m -> m
forall a. NBody a -> Bool
forall a. NBody a -> Int
forall a. NBody a -> [a]
forall a. (a -> a -> a) -> NBody a -> a
forall m a. Monoid m => (a -> m) -> NBody a -> m
forall b a. (b -> a -> b) -> b -> NBody a -> b
forall a b. (a -> b -> b) -> b -> NBody a -> b
forall (t :: * -> *).
(forall m. Monoid m => t m -> m)
-> (forall m a. Monoid m => (a -> m) -> t a -> m)
-> (forall m a. Monoid m => (a -> m) -> t a -> m)
-> (forall a b. (a -> b -> b) -> b -> t a -> b)
-> (forall a b. (a -> b -> b) -> b -> t a -> b)
-> (forall b a. (b -> a -> b) -> b -> t a -> b)
-> (forall b a. (b -> a -> b) -> b -> t a -> b)
-> (forall a. (a -> a -> a) -> t a -> a)
-> (forall a. (a -> a -> a) -> t a -> a)
-> (forall a. t a -> [a])
-> (forall a. t a -> Bool)
-> (forall a. t a -> Int)
-> (forall a. Eq a => a -> t a -> Bool)
-> (forall a. Ord a => t a -> a)
-> (forall a. Ord a => t a -> a)
-> (forall a. Num a => t a -> a)
-> (forall a. Num a => t a -> a)
-> Foldable t
$cfold :: forall m. Monoid m => NBody m -> m
fold :: forall m. Monoid m => NBody m -> m
$cfoldMap :: forall m a. Monoid m => (a -> m) -> NBody a -> m
foldMap :: forall m a. Monoid m => (a -> m) -> NBody a -> m
$cfoldMap' :: forall m a. Monoid m => (a -> m) -> NBody a -> m
foldMap' :: forall m a. Monoid m => (a -> m) -> NBody a -> m
$cfoldr :: forall a b. (a -> b -> b) -> b -> NBody a -> b
foldr :: forall a b. (a -> b -> b) -> b -> NBody a -> b
$cfoldr' :: forall a b. (a -> b -> b) -> b -> NBody a -> b
foldr' :: forall a b. (a -> b -> b) -> b -> NBody a -> b
$cfoldl :: forall b a. (b -> a -> b) -> b -> NBody a -> b
foldl :: forall b a. (b -> a -> b) -> b -> NBody a -> b
$cfoldl' :: forall b a. (b -> a -> b) -> b -> NBody a -> b
foldl' :: forall b a. (b -> a -> b) -> b -> NBody a -> b
$cfoldr1 :: forall a. (a -> a -> a) -> NBody a -> a
foldr1 :: forall a. (a -> a -> a) -> NBody a -> a
$cfoldl1 :: forall a. (a -> a -> a) -> NBody a -> a
foldl1 :: forall a. (a -> a -> a) -> NBody a -> a
$ctoList :: forall a. NBody a -> [a]
toList :: forall a. NBody a -> [a]
$cnull :: forall a. NBody a -> Bool
null :: forall a. NBody a -> Bool
$clength :: forall a. NBody a -> Int
length :: forall a. NBody a -> Int
$celem :: forall a. Eq a => a -> NBody a -> Bool
elem :: forall a. Eq a => a -> NBody a -> Bool
$cmaximum :: forall a. Ord a => NBody a -> a
maximum :: forall a. Ord a => NBody a -> a
$cminimum :: forall a. Ord a => NBody a -> a
minimum :: forall a. Ord a => NBody a -> a
$csum :: forall a. Num a => NBody a -> a
sum :: forall a. Num a => NBody a -> a
$cproduct :: forall a. Num a => NBody a -> a
product :: forall a. Num a => NBody a -> a
Foldable)

instance Applicative NBody where
    pure :: forall a. a -> NBody a
pure a
a = [a] -> NBody a
forall a. [a] -> NBody a
NBody (a -> [a]
forall a. a -> [a]
repeat a
a)
    liftA2 :: forall a b c. (a -> b -> c) -> NBody a -> NBody b -> NBody c
liftA2 a -> b -> c
f (NBody [a]
xs) (NBody [b]
ys) = [c] -> NBody c
forall a. [a] -> NBody a
NBody ((a -> b -> c) -> [a] -> [b] -> [c]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> b -> c
f [a]
xs [b]
ys)

instance Traversable NBody where
    traverse :: forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> NBody a -> f (NBody b)
traverse a -> f b
f (NBody [a]
xs) = [b] -> NBody b
forall a. [a] -> NBody a
NBody ([b] -> NBody b) -> f [b] -> f (NBody b)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (a -> f b) -> [a] -> f [b]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse a -> f b
f [a]
xs

instance VectorSpace a => VectorSpace (NBody a) where
    +. :: NBody a -> NBody a -> NBody a
(+.) = (a -> a -> a) -> NBody a -> NBody a -> NBody a
forall a b c. (a -> b -> c) -> NBody a -> NBody b -> NBody c
forall (f :: * -> *) a b c.
Applicative f =>
(a -> b -> c) -> f a -> f b -> f c
liftA2 a -> a -> a
forall v. VectorSpace v => v -> v -> v
(+.)
    Double
a *. :: Double -> NBody a -> NBody a
*. NBody a
xs = (a -> a) -> NBody a -> NBody a
forall a b. (a -> b) -> NBody a -> NBody b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Double
a Double -> a -> a
forall v. VectorSpace v => Double -> v -> v
*.) NBody a
xs
    -. :: NBody a -> NBody a -> NBody a
(-.) = (a -> a -> a) -> NBody a -> NBody a -> NBody a
forall a b c. (a -> b -> c) -> NBody a -> NBody b -> NBody c
forall (f :: * -> *) a b c.
Applicative f =>
(a -> b -> c) -> f a -> f b -> f c
liftA2 a -> a -> a
forall v. VectorSpace v => v -> v -> v
(-.)
    zero :: NBody a
zero = [a] -> NBody a
forall a. [a] -> NBody a
NBody (a -> [a]
forall a. a -> [a]
repeat a
forall v. VectorSpace v => v
zero)

nBody
    :: (Vec2 -> Double)
    -- ^ External potential, affecting all particles. If you don't require an external potential, use 'zero'.
    -> (Vec2 -> Vec2 -> Double)
    -- ^ Interaction term between two particles. Should be symmetric.
    -> NBody Double
    -- ^ Particle masses
    -> NBody PhaseSpace -> NBody PhaseSpace
nBody :: (Vec2 -> Double)
-> (Vec2 -> Vec2 -> Double)
-> NBody Double
-> NBody PhaseSpace
-> NBody PhaseSpace
nBody Vec2 -> Double
potential Vec2 -> Vec2 -> Double
interaction NBody Double
masses NBody PhaseSpace
particles = [PhaseSpace] -> NBody PhaseSpace
forall a. [a] -> NBody a
NBody
    [ Vec2 -> Vec2 -> PhaseSpace
PhaseSpace Vec2
deltaP Vec2
deltaQ
    | (Integer
i, Double
m1, PhaseSpace Vec2
p1 Vec2
q1) <- [Integer]
-> [Double] -> [PhaseSpace] -> [(Integer, Double, PhaseSpace)]
forall a b c. [a] -> [b] -> [c] -> [(a, b, c)]
zip3 [Integer
0..] (NBody Double -> [Double]
forall a. NBody a -> [a]
getNBody NBody Double
masses) (NBody PhaseSpace -> [PhaseSpace]
forall a. NBody a -> [a]
getNBody NBody PhaseSpace
particles)
    , let deltaQ :: Vec2
deltaQ = Vec2
p1 Vec2 -> Double -> Vec2
forall v. VectorSpace v => v -> Double -> v
/. Double
m1
    , let deltaP :: Vec2
deltaP = Vec2 -> Vec2
forall v. VectorSpace v => v -> v
negateV (Vec2 -> Vec2) -> Vec2 -> Vec2
forall a b. (a -> b) -> a -> b
$ (Vec2 -> Double) -> Vec2 -> Vec2
grad Vec2 -> Double
potential Vec2
q1 Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. [Vec2] -> Vec2
forall a. VectorSpace a => [a] -> a
vsum
            [ (Vec2 -> Double) -> Vec2 -> Vec2
grad (Vec2 -> Vec2 -> Double
interaction Vec2
q2) Vec2
q1
            | (Integer
j, PhaseSpace Vec2
_ Vec2
q2) <- [Integer] -> [PhaseSpace] -> [(Integer, PhaseSpace)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Integer
0..] (NBody PhaseSpace -> [PhaseSpace]
forall a. NBody a -> [a]
getNBody NBody PhaseSpace
particles)
            , Integer
i Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
j
            ]
    ]