{-# OPTIONS -Wall #-}

{- | 
Module      :  LPFPCore.MOExamples
Copyright   :  (c) Scott N. Walck 2023
License     :  BSD3 (see LICENSE)
Maintainer  :  Scott N. Walck <walck@lvc.edu>
Stability   :  stable

Code from chapter 20 of the book Learn Physics with Functional Programming
-}

module LPFPCore.MOExamples where

import LPFPCore.SimpleVec
    ( R, Vec, (^+^), (^-^), (*^), vec, zeroV, magnitude
    , sumV, iHat, jHat, kHat, zComp )
import LPFPCore.Mechanics1D ( TimeStep, NumericalMethod, euler, rungeKutta4 )
import LPFPCore.Mechanics3D
    ( ParticleState(..), HasTime(..), defaultParticleState
    , earthSurfaceGravity )
import LPFPCore.MultipleObjects
    ( MultiParticleState(..), DMultiParticleState, Force(..), TwoBodyForce
    , newtonSecondMPS, updateMPS, statesMPS, eulerCromerMPS
    , linearSpring, fixedLinearSpring, billiardForce )

twoSpringsForces :: [Force]
twoSpringsForces :: [Force]
twoSpringsForces
    = [Int -> OneBodyForce -> Force
ExternalForce Int
0 (R -> R -> Vec -> OneBodyForce
fixedLinearSpring R
100 R
0.5 Vec
zeroV)
      ,Int -> Int -> TwoBodyForce -> Force
InternalForce Int
0 Int
1 (R -> R -> TwoBodyForce
linearSpring R
100 R
0.5)
      ,Int -> OneBodyForce -> Force
ExternalForce Int
0 OneBodyForce
earthSurfaceGravity
      ,Int -> OneBodyForce -> Force
ExternalForce Int
1 OneBodyForce
earthSurfaceGravity
      ]

twoSpringsInitial :: MultiParticleState
twoSpringsInitial :: MultiParticleState
twoSpringsInitial
    = [ParticleState] -> MultiParticleState
MPS [ParticleState
defaultParticleState
           { mass   = 2
           , posVec = 0.4 *^ jHat ^-^ 0.3 *^ kHat }
          ,ParticleState
defaultParticleState
           { mass   = 3
           , posVec = 0.4 *^ jHat ^-^ 0.8 *^ kHat }
          ]

twoSpringsUpdate :: TimeStep
                 -> MultiParticleState  -- old state
                 -> MultiParticleState  -- new state
twoSpringsUpdate :: R -> MultiParticleState -> MultiParticleState
twoSpringsUpdate R
dt = NumericalMethod MultiParticleState DMultiParticleState
-> [Force] -> MultiParticleState -> MultiParticleState
updateMPS (R -> NumericalMethod MultiParticleState DMultiParticleState
eulerCromerMPS R
dt) [Force]
twoSpringsForces

kineticEnergy :: ParticleState -> R
kineticEnergy :: ParticleState -> R
kineticEnergy ParticleState
st = let m :: R
m = ParticleState -> R
mass ParticleState
st
                       v :: R
v = Vec -> R
magnitude (OneBodyForce
velocity ParticleState
st)
                   in (R
1R -> R -> R
forall a. Fractional a => a -> a -> a
/R
2) R -> R -> R
forall a. Num a => a -> a -> a
* R
m R -> R -> R
forall a. Num a => a -> a -> a
* R
vR -> R -> R
forall a. Floating a => a -> a -> a
**R
2

systemKE :: MultiParticleState -> R
systemKE :: MultiParticleState -> R
systemKE (MPS [ParticleState]
sts) = [R] -> R
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ParticleState -> R
kineticEnergy ParticleState
st | ParticleState
st <- [ParticleState]
sts]

linearSpringPE :: R              -- spring constant
               -> R              -- equilibrium length
               -> ParticleState  -- state of particle at one end of spring
               -> ParticleState  -- state of particle at other end of spring
               -> R              -- potential energy of the spring
linearSpringPE :: R -> R -> ParticleState -> ParticleState -> R
linearSpringPE R
k R
re ParticleState
st1 ParticleState
st2
    = let r1 :: Vec
r1 = OneBodyForce
posVec ParticleState
st1
          r2 :: Vec
r2 = OneBodyForce
posVec ParticleState
st2
          r21 :: Vec
r21 = Vec
r2 Vec -> Vec -> Vec
^-^ Vec
r1
          r21mag :: R
r21mag = Vec -> R
magnitude Vec
r21
      in R
k R -> R -> R
forall a. Num a => a -> a -> a
* (R
r21mag R -> R -> R
forall a. Num a => a -> a -> a
- R
re)R -> R -> R
forall a. Floating a => a -> a -> a
**R
2 R -> R -> R
forall a. Fractional a => a -> a -> a
/ R
2

-- z direction is toward the sky
-- assumes SI units
earthSurfaceGravityPE :: ParticleState -> R
earthSurfaceGravityPE :: ParticleState -> R
earthSurfaceGravityPE ParticleState
st
    = let g :: R
g = R
9.80665  -- m/s^2
          m :: R
m = ParticleState -> R
mass ParticleState
st
          z :: R
z = Vec -> R
zComp (OneBodyForce
posVec ParticleState
st)
      in R
m R -> R -> R
forall a. Num a => a -> a -> a
* R
g R -> R -> R
forall a. Num a => a -> a -> a
* R
z

twoSpringsPE :: MultiParticleState -> R
twoSpringsPE :: MultiParticleState -> R
twoSpringsPE (MPS [ParticleState]
sts)
    = R -> R -> ParticleState -> ParticleState -> R
linearSpringPE R
100 R
0.5 ParticleState
defaultParticleState ([ParticleState]
sts [ParticleState] -> Int -> ParticleState
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)
      R -> R -> R
forall a. Num a => a -> a -> a
+ R -> R -> ParticleState -> ParticleState -> R
linearSpringPE R
100 R
0.5 ([ParticleState]
sts [ParticleState] -> Int -> ParticleState
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) ([ParticleState]
sts [ParticleState] -> Int -> ParticleState
forall a. HasCallStack => [a] -> Int -> a
!! Int
1)
      R -> R -> R
forall a. Num a => a -> a -> a
+ ParticleState -> R
earthSurfaceGravityPE ([ParticleState]
sts [ParticleState] -> Int -> ParticleState
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)
      R -> R -> R
forall a. Num a => a -> a -> a
+ ParticleState -> R
earthSurfaceGravityPE ([ParticleState]
sts [ParticleState] -> Int -> ParticleState
forall a. HasCallStack => [a] -> Int -> a
!! Int
1)

twoSpringsME :: MultiParticleState -> R
twoSpringsME :: MultiParticleState -> R
twoSpringsME MultiParticleState
mpst = MultiParticleState -> R
systemKE MultiParticleState
mpst R -> R -> R
forall a. Num a => a -> a -> a
+ MultiParticleState -> R
twoSpringsPE MultiParticleState
mpst

billiardForces :: R -> [Force]
billiardForces :: R -> [Force]
billiardForces R
k = [Int -> Int -> TwoBodyForce -> Force
InternalForce Int
0 Int
1 (R -> R -> TwoBodyForce
billiardForce R
k (R
2R -> R -> R
forall a. Num a => a -> a -> a
*R
ballRadius))]

ballRadius :: R
ballRadius :: R
ballRadius = R
0.03  -- 6cm diameter = 0.03m radius

billiardDiffEq :: R -> MultiParticleState -> DMultiParticleState
billiardDiffEq :: R -> MultiParticleState -> DMultiParticleState
billiardDiffEq R
k = [Force] -> MultiParticleState -> DMultiParticleState
newtonSecondMPS ([Force] -> MultiParticleState -> DMultiParticleState)
-> [Force] -> MultiParticleState -> DMultiParticleState
forall a b. (a -> b) -> a -> b
$ R -> [Force]
billiardForces R
k

billiardUpdate
    :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState)
    -> R         -- k
    -> TimeStep  -- dt
    -> MultiParticleState -> MultiParticleState
billiardUpdate :: (R -> NumericalMethod MultiParticleState DMultiParticleState)
-> R -> R -> MultiParticleState -> MultiParticleState
billiardUpdate R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod R
k R
dt = NumericalMethod MultiParticleState DMultiParticleState
-> [Force] -> MultiParticleState -> MultiParticleState
updateMPS (R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod R
dt) (R -> [Force]
billiardForces R
k)

billiardEvolver
    :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState)
    -> R         -- k
    -> TimeStep  -- dt
    -> MultiParticleState -> [MultiParticleState]
billiardEvolver :: (R -> NumericalMethod MultiParticleState DMultiParticleState)
-> R -> R -> MultiParticleState -> [MultiParticleState]
billiardEvolver R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod R
k R
dt = NumericalMethod MultiParticleState DMultiParticleState
-> [Force] -> MultiParticleState -> [MultiParticleState]
statesMPS (R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod R
dt) (R -> [Force]
billiardForces R
k)

billiardInitial :: MultiParticleState
billiardInitial :: MultiParticleState
billiardInitial
    = let ballMass :: R
ballMass = R
0.160  -- 160g
      in [ParticleState] -> MultiParticleState
MPS [ParticleState
defaultParticleState { mass     = ballMass
                                   , posVec   = zeroV
                                   , velocity = 0.2 *^ iHat }
             ,ParticleState
defaultParticleState { mass     = ballMass
                                   , posVec   = iHat ^+^ 0.02 *^ jHat
                                   , velocity = zeroV }
             ]

billiardStates
    :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState)
    -> R         -- k
    -> TimeStep  -- dt
    -> [MultiParticleState]
billiardStates :: (R -> NumericalMethod MultiParticleState DMultiParticleState)
-> R -> R -> [MultiParticleState]
billiardStates R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod R
k R
dt
    = NumericalMethod MultiParticleState DMultiParticleState
-> [Force] -> MultiParticleState -> [MultiParticleState]
statesMPS (R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod R
dt) (R -> [Force]
billiardForces R
k) MultiParticleState
billiardInitial

billiardStatesFinite
    :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState)
    -> R         -- k
    -> TimeStep  -- dt
    -> [MultiParticleState]
billiardStatesFinite :: (R -> NumericalMethod MultiParticleState DMultiParticleState)
-> R -> R -> [MultiParticleState]
billiardStatesFinite R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod R
k R
dt
    = (MultiParticleState -> Bool)
-> [MultiParticleState] -> [MultiParticleState]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (\MultiParticleState
st -> MultiParticleState -> R
forall s. HasTime s => s -> R
timeOf MultiParticleState
st R -> R -> Bool
forall a. Ord a => a -> a -> Bool
<= R
10) ((R -> NumericalMethod MultiParticleState DMultiParticleState)
-> R -> R -> [MultiParticleState]
billiardStates R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod R
k R
dt)

momentum :: ParticleState -> Vec
momentum :: OneBodyForce
momentum ParticleState
st = let m :: R
m = ParticleState -> R
mass ParticleState
st
                  v :: Vec
v = OneBodyForce
velocity ParticleState
st
              in R
m R -> Vec -> Vec
*^ Vec
v

systemP :: MultiParticleState -> Vec
systemP :: MultiParticleState -> Vec
systemP (MPS [ParticleState]
sts) = [Vec] -> Vec
sumV [OneBodyForce
momentum ParticleState
st | ParticleState
st <- [ParticleState]
sts]

percentChangePMag :: [MultiParticleState] -> R
percentChangePMag :: [MultiParticleState] -> R
percentChangePMag [MultiParticleState]
mpsts
    = let p0 :: Vec
p0 = MultiParticleState -> Vec
systemP ([MultiParticleState] -> MultiParticleState
forall a. HasCallStack => [a] -> a
head [MultiParticleState]
mpsts)
          p1 :: Vec
p1 = MultiParticleState -> Vec
systemP ([MultiParticleState] -> MultiParticleState
forall a. HasCallStack => [a] -> a
last [MultiParticleState]
mpsts)
      in R
100 R -> R -> R
forall a. Num a => a -> a -> a
* Vec -> R
magnitude (Vec
p1 Vec -> Vec -> Vec
^-^ Vec
p0) R -> R -> R
forall a. Fractional a => a -> a -> a
/ Vec -> R
magnitude Vec
p0

sigFigs :: Int -> R -> Float
sigFigs :: Int -> R -> Float
sigFigs Int
n R
x = let expon :: Int
                  expon :: Int
expon = R -> Int
forall b. Integral b => R -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (R -> R -> R
forall a. Floating a => a -> a -> a
logBase R
10 R
x) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
                  toInt :: R -> Int
                  toInt :: R -> Int
toInt = R -> Int
forall b. Integral b => R -> b
forall a b. (RealFrac a, Integral b) => a -> b
round
              in (Float
10Float -> Int -> Float
forall a b. (Fractional a, Integral b) => a -> b -> a
^^Int
expon Float -> Float -> Float
forall a. Num a => a -> a -> a
*) (Float -> Float) -> Float -> Float
forall a b. (a -> b) -> a -> b
$ Int -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Float) -> Int -> Float
forall a b. (a -> b) -> a -> b
$ R -> Int
toInt (R
10R -> Int -> R
forall a b. (Fractional a, Integral b) => a -> b -> a
^^(-Int
expon) R -> R -> R
forall a. Num a => a -> a -> a
* R
x)

data Justification = LJ | RJ deriving Int -> Justification -> ShowS
[Justification] -> ShowS
Justification -> String
(Int -> Justification -> ShowS)
-> (Justification -> String)
-> ([Justification] -> ShowS)
-> Show Justification
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Justification -> ShowS
showsPrec :: Int -> Justification -> ShowS
$cshow :: Justification -> String
show :: Justification -> String
$cshowList :: [Justification] -> ShowS
showList :: [Justification] -> ShowS
Show

data Table a = Table Justification [[a]]

instance Show a => Show (Table a) where
    show :: Table a -> String
show (Table Justification
j [[a]]
xss)
        = let pairWithLength :: p -> (String, Int)
pairWithLength p
x = let str :: String
str = p -> String
forall a. Show a => a -> String
show p
x in (String
str, String -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length String
str)
              pairss :: [[(String, Int)]]
pairss = ([a] -> [(String, Int)]) -> [[a]] -> [[(String, Int)]]
forall a b. (a -> b) -> [a] -> [b]
map ((a -> (String, Int)) -> [a] -> [(String, Int)]
forall a b. (a -> b) -> [a] -> [b]
map a -> (String, Int)
forall {p}. Show p => p -> (String, Int)
pairWithLength) [[a]]
xss
              maxLength :: Int
maxLength = [Int] -> Int
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (([Int] -> Int) -> [[Int]] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map [Int] -> Int
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (([(String, Int)] -> [Int]) -> [[(String, Int)]] -> [[Int]]
forall a b. (a -> b) -> [a] -> [b]
map (((String, Int) -> Int) -> [(String, Int)] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (String, Int) -> Int
forall a b. (a, b) -> b
snd) [[(String, Int)]]
pairss))
              showPair :: (String, Int) -> String
showPair (String
str,Int
len)
                  = case Justification
j of
                      Justification
LJ -> String
str String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> Char -> String
forall a. Int -> a -> [a]
replicate (Int
maxLength Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
len) Char
' '
                      Justification
RJ -> Int -> Char -> String
forall a. Int -> a -> [a]
replicate (Int
maxLength Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
len) Char
' ' String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
str
              showLine :: t (String, Int) -> String
showLine t (String, Int)
pairs = ((String, Int) -> String) -> t (String, Int) -> String
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (String, Int) -> String
showPair t (String, Int)
pairs String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"\n"
          in ShowS
forall a. HasCallStack => [a] -> [a]
init ShowS -> ShowS
forall a b. (a -> b) -> a -> b
$ ([(String, Int)] -> String) -> [[(String, Int)]] -> String
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap [(String, Int)] -> String
forall {t :: * -> *}. Foldable t => t (String, Int) -> String
showLine [[(String, Int)]]
pairss

pTable :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState)
       -> [R]         -- ks
       -> [TimeStep]  -- dts
       -> Table Float
pTable :: (R -> NumericalMethod MultiParticleState DMultiParticleState)
-> [R] -> [R] -> Table Float
pTable R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod [R]
ks [R]
dts
    = Justification -> [[Float]] -> Table Float
forall a. Justification -> [[a]] -> Table a
Table Justification
LJ [[Int -> R -> Float
sigFigs Int
2 (R -> Float) -> R -> Float
forall a b. (a -> b) -> a -> b
$
                 [MultiParticleState] -> R
percentChangePMag ((R -> NumericalMethod MultiParticleState DMultiParticleState)
-> R -> R -> [MultiParticleState]
billiardStatesFinite R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod R
k R
dt)
                     | R
dt <- [R]
dts] | R
k <- [R]
ks]

pTableEu :: [R]         -- ks
         -> [TimeStep]  -- dts
         -> Table Float
pTableEu :: [R] -> [R] -> Table Float
pTableEu = (R -> NumericalMethod MultiParticleState DMultiParticleState)
-> [R] -> [R] -> Table Float
pTable R -> NumericalMethod MultiParticleState DMultiParticleState
forall s ds. Diff s ds => R -> (s -> ds) -> s -> s
euler

percentChangeKE :: [MultiParticleState] -> R
percentChangeKE :: [MultiParticleState] -> R
percentChangeKE [MultiParticleState]
mpsts
    = let ke0 :: R
ke0 = MultiParticleState -> R
systemKE ([MultiParticleState] -> MultiParticleState
forall a. HasCallStack => [a] -> a
head [MultiParticleState]
mpsts)
          ke1 :: R
ke1 = MultiParticleState -> R
systemKE ([MultiParticleState] -> MultiParticleState
forall a. HasCallStack => [a] -> a
last [MultiParticleState]
mpsts)
      in R
100 R -> R -> R
forall a. Num a => a -> a -> a
* (R
ke1 R -> R -> R
forall a. Num a => a -> a -> a
- R
ke0) R -> R -> R
forall a. Fractional a => a -> a -> a
/ R
ke0

tenths :: R -> Float
tenths :: R -> Float
tenths = let toInt :: R -> Int
             toInt :: R -> Int
toInt = R -> Int
forall b. Integral b => R -> b
forall a b. (RealFrac a, Integral b) => a -> b
round
         in (Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/ Float
10) (Float -> Float) -> (R -> Float) -> R -> Float
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Float) -> (R -> Int) -> R -> Float
forall b c a. (b -> c) -> (a -> b) -> a -> c
. R -> Int
toInt (R -> Int) -> (R -> R) -> R -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (R -> R -> R
forall a. Num a => a -> a -> a
* R
10)

keTable
    :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState)
    -> [R]         -- ks
    -> [TimeStep]  -- dts
    -> Table Float
keTable :: (R -> NumericalMethod MultiParticleState DMultiParticleState)
-> [R] -> [R] -> Table Float
keTable R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod [R]
ks [R]
dts
    = Justification -> [[Float]] -> Table Float
forall a. Justification -> [[a]] -> Table a
Table Justification
RJ [[R -> Float
tenths (R -> Float) -> R -> Float
forall a b. (a -> b) -> a -> b
$
                 [MultiParticleState] -> R
percentChangeKE ((R -> NumericalMethod MultiParticleState DMultiParticleState)
-> R -> R -> [MultiParticleState]
billiardStatesFinite R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod R
k R
dt)
                     | R
dt <- [R]
dts] | R
k <- [R]
ks]

contactSteps :: [MultiParticleState] -> Int
contactSteps :: [MultiParticleState] -> Int
contactSteps = [MultiParticleState] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([MultiParticleState] -> Int)
-> ([MultiParticleState] -> [MultiParticleState])
-> [MultiParticleState]
-> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (MultiParticleState -> Bool)
-> [MultiParticleState] -> [MultiParticleState]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile MultiParticleState -> Bool
inContact ([MultiParticleState] -> [MultiParticleState])
-> ([MultiParticleState] -> [MultiParticleState])
-> [MultiParticleState]
-> [MultiParticleState]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (MultiParticleState -> Bool)
-> [MultiParticleState] -> [MultiParticleState]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Bool -> Bool
not (Bool -> Bool)
-> (MultiParticleState -> Bool) -> MultiParticleState -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. MultiParticleState -> Bool
inContact)

inContact :: MultiParticleState -> Bool
inContact :: MultiParticleState -> Bool
inContact (MPS [ParticleState]
sts)
    = let r :: R
r = Vec -> R
magnitude (Vec -> R) -> Vec -> R
forall a b. (a -> b) -> a -> b
$ OneBodyForce
posVec ([ParticleState]
sts [ParticleState] -> Int -> ParticleState
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) Vec -> Vec -> Vec
^-^ OneBodyForce
posVec ([ParticleState]
sts [ParticleState] -> Int -> ParticleState
forall a. HasCallStack => [a] -> Int -> a
!! Int
1)
      in R
r R -> R -> Bool
forall a. Ord a => a -> a -> Bool
< R
2 R -> R -> R
forall a. Num a => a -> a -> a
* R
ballRadius

contactTable
    :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState)
    -> [R]         -- ks
    -> [TimeStep]  -- dts
    -> Table Int
contactTable :: (R -> NumericalMethod MultiParticleState DMultiParticleState)
-> [R] -> [R] -> Table Int
contactTable R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod [R]
ks [R]
dts
    = Justification -> [[Int]] -> Table Int
forall a. Justification -> [[a]] -> Table a
Table Justification
RJ [[[MultiParticleState] -> Int
contactSteps ((R -> NumericalMethod MultiParticleState DMultiParticleState)
-> R -> R -> [MultiParticleState]
billiardStatesFinite R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod R
k R
dt)
                     | R
dt <- [R]
dts] | R
k <- [R]
ks]

closest :: [MultiParticleState] -> R
closest :: [MultiParticleState] -> R
closest = [R] -> R
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum ([R] -> R)
-> ([MultiParticleState] -> [R]) -> [MultiParticleState] -> R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (MultiParticleState -> R) -> [MultiParticleState] -> [R]
forall a b. (a -> b) -> [a] -> [b]
map MultiParticleState -> R
separation

separation :: MultiParticleState -> R
separation :: MultiParticleState -> R
separation (MPS [ParticleState]
sts)
    = Vec -> R
magnitude (Vec -> R) -> Vec -> R
forall a b. (a -> b) -> a -> b
$ OneBodyForce
posVec ([ParticleState]
sts [ParticleState] -> Int -> ParticleState
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) Vec -> Vec -> Vec
^-^ OneBodyForce
posVec ([ParticleState]
sts [ParticleState] -> Int -> ParticleState
forall a. HasCallStack => [a] -> Int -> a
!! Int
1)

closestTable
    :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState)
    -> [R]         -- ks
    -> [TimeStep]  -- dts
    -> Table Float
closestTable :: (R -> NumericalMethod MultiParticleState DMultiParticleState)
-> [R] -> [R] -> Table Float
closestTable R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod [R]
ks [R]
dts
    = Justification -> [[Float]] -> Table Float
forall a. Justification -> [[a]] -> Table a
Table Justification
RJ [[R -> Float
tenths (R -> Float) -> R -> Float
forall a b. (a -> b) -> a -> b
$ (R
100R -> R -> R
forall a. Num a => a -> a -> a
*) (R -> R) -> R -> R
forall a b. (a -> b) -> a -> b
$
                 [MultiParticleState] -> R
closest ((R -> NumericalMethod MultiParticleState DMultiParticleState)
-> R -> R -> [MultiParticleState]
billiardStatesFinite R -> NumericalMethod MultiParticleState DMultiParticleState
nMethod R
k R
dt)
                     | R
dt <- [R]
dts] | R
k <- [R]
ks]

-- 64 masses (0 to 63)
-- There are 63 internal springs, 2 external springs
forcesString :: [Force]
forcesString :: [Force]
forcesString
    = [Int -> OneBodyForce -> Force
ExternalForce  Int
0 (R -> R -> Vec -> OneBodyForce
fixedLinearSpring R
5384 R
0 (R -> R -> R -> Vec
vec    R
0 R
0 R
0))
      ,Int -> OneBodyForce -> Force
ExternalForce Int
63 (R -> R -> Vec -> OneBodyForce
fixedLinearSpring R
5384 R
0 (R -> R -> R -> Vec
vec R
0.65 R
0 R
0))] [Force] -> [Force] -> [Force]
forall a. [a] -> [a] -> [a]
++
      [Int -> Int -> TwoBodyForce -> Force
InternalForce Int
n (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (R -> R -> TwoBodyForce
linearSpring R
5384 R
0) | Int
n <- [Int
0..Int
62]]

stringUpdate :: TimeStep
             -> MultiParticleState  -- old state
             -> MultiParticleState  -- new state
stringUpdate :: R -> MultiParticleState -> MultiParticleState
stringUpdate R
dt = NumericalMethod MultiParticleState DMultiParticleState
-> [Force] -> MultiParticleState -> MultiParticleState
updateMPS (R -> NumericalMethod MultiParticleState DMultiParticleState
forall s ds. Diff s ds => R -> (s -> ds) -> s -> s
rungeKutta4 R
dt) [Force]
forcesString

stringInitialOvertone :: Int -> MultiParticleState
stringInitialOvertone :: Int -> MultiParticleState
stringInitialOvertone Int
n
    = [ParticleState] -> MultiParticleState
MPS [ParticleState
defaultParticleState
           { mass     = 0.8293e-3 * 0.65 / 64
           , posVec   = x *^ iHat ^+^ y *^ jHat
           , velocity = zeroV
           } | R
x <- [R
0.01, R
0.02 .. R
0.64],
           let y :: R
y = R
0.005 R -> R -> R
forall a. Num a => a -> a -> a
* R -> R
forall a. Floating a => a -> a
sin (Int -> R
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n R -> R -> R
forall a. Num a => a -> a -> a
* R
forall a. Floating a => a
pi R -> R -> R
forall a. Num a => a -> a -> a
* R
x R -> R -> R
forall a. Fractional a => a -> a -> a
/ R
0.65)]

stringInitialPluck :: MultiParticleState
stringInitialPluck :: MultiParticleState
stringInitialPluck = [ParticleState] -> MultiParticleState
MPS [ParticleState
defaultParticleState
             { mass     = 0.8293e-3 * 0.65 / 64
             , posVec   = x *^ iHat ^+^ y *^ jHat
             , velocity = zeroV
             } | R
x <- [R
0.01, R
0.02 .. R
0.64], let y :: R
y = R -> R
pluckEq R
x]
    where
      pluckEq :: R -> R
      pluckEq :: R -> R
pluckEq R
x
          | R
x R -> R -> Bool
forall a. Ord a => a -> a -> Bool
<= R
0.51  = R
0.005 R -> R -> R
forall a. Fractional a => a -> a -> a
/ (R
0.51 R -> R -> R
forall a. Num a => a -> a -> a
- R
0.00) R -> R -> R
forall a. Num a => a -> a -> a
* (R
x R -> R -> R
forall a. Num a => a -> a -> a
- R
0.00)
          | Bool
otherwise  = R
0.005 R -> R -> R
forall a. Fractional a => a -> a -> a
/ (R
0.51 R -> R -> R
forall a. Num a => a -> a -> a
- R
0.65) R -> R -> R
forall a. Num a => a -> a -> a
* (R
x R -> R -> R
forall a. Num a => a -> a -> a
- R
0.65)

mpsPos :: MultiParticleState -> IO ()
mpsPos :: MultiParticleState -> IO ()
mpsPos = MultiParticleState -> IO ()
forall a. HasCallStack => a
undefined

mpsVel :: MultiParticleState -> IO ()
mpsVel :: MultiParticleState -> IO ()
mpsVel = MultiParticleState -> IO ()
forall a. HasCallStack => a
undefined

dissipation :: R  -- damping constant
            -> R  -- threshold center separation
            -> TwoBodyForce
dissipation :: R -> R -> TwoBodyForce
dissipation R
b R
re ParticleState
st1 ParticleState
st2
    = let r1 :: Vec
r1 = OneBodyForce
posVec ParticleState
st1
          r2 :: Vec
r2 = OneBodyForce
posVec ParticleState
st2
          v1 :: Vec
v1 = OneBodyForce
velocity ParticleState
st1
          v2 :: Vec
v2 = OneBodyForce
velocity ParticleState
st2
          r21 :: Vec
r21 = Vec
r2 Vec -> Vec -> Vec
^-^ Vec
r1
          v21 :: Vec
v21 = Vec
v2 Vec -> Vec -> Vec
^-^ Vec
v1
      in if Vec -> R
magnitude Vec
r21 R -> R -> Bool
forall a. Ord a => a -> a -> Bool
>= R
re
         then Vec
zeroV
         else (-R
b) R -> Vec -> Vec
*^ Vec
v21