{-# 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 :: R
mass   = R
2
           , posVec :: Vec
posVec = R
0.4 R -> Vec -> Vec
*^ Vec
jHat Vec -> Vec -> Vec
^-^ R
0.3 R -> Vec -> Vec
*^ Vec
kHat }
          ,ParticleState
defaultParticleState
           { mass :: R
mass   = R
3
           , posVec :: Vec
posVec = R
0.4 R -> Vec -> Vec
*^ Vec
jHat Vec -> Vec -> Vec
^-^ R
0.8 R -> Vec -> Vec
*^ Vec
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
1forall a. Fractional a => a -> a -> a
/R
2) forall a. Num a => a -> a -> a
* R
m forall a. Num a => a -> a -> a
* R
vforall a. Floating a => a -> a -> a
**R
2

systemKE :: MultiParticleState -> R
systemKE :: MultiParticleState -> R
systemKE (MPS [ParticleState]
sts) = 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 forall a. Num a => a -> a -> a
* (R
r21mag forall a. Num a => a -> a -> a
- R
re)forall a. Floating a => a -> a -> a
**R
2 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 forall a. Num a => a -> a -> a
* R
g 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 forall a. [a] -> Int -> a
!! Int
0)
      forall a. Num a => a -> a -> a
+ R -> R -> ParticleState -> ParticleState -> R
linearSpringPE R
100 R
0.5 ([ParticleState]
sts forall a. [a] -> Int -> a
!! Int
0) ([ParticleState]
sts forall a. [a] -> Int -> a
!! Int
1)
      forall a. Num a => a -> a -> a
+ ParticleState -> R
earthSurfaceGravityPE ([ParticleState]
sts forall a. [a] -> Int -> a
!! Int
0)
      forall a. Num a => a -> a -> a
+ ParticleState -> R
earthSurfaceGravityPE ([ParticleState]
sts forall a. [a] -> Int -> a
!! Int
1)

twoSpringsME :: MultiParticleState -> R
twoSpringsME :: MultiParticleState -> R
twoSpringsME MultiParticleState
mpst = MultiParticleState -> R
systemKE MultiParticleState
mpst 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
2forall 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 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 :: R
mass     = R
ballMass
                                   , posVec :: Vec
posVec   = Vec
zeroV
                                   , velocity :: Vec
velocity = R
0.2 R -> Vec -> Vec
*^ Vec
iHat }
             ,ParticleState
defaultParticleState { mass :: R
mass     = R
ballMass
                                   , posVec :: Vec
posVec   = Vec
iHat Vec -> Vec -> Vec
^+^ R
0.02 R -> Vec -> Vec
*^ Vec
jHat
                                   , velocity :: Vec
velocity = Vec
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
    = forall a. (a -> Bool) -> [a] -> [a]
takeWhile (\MultiParticleState
st -> forall s. HasTime s => s -> R
timeOf MultiParticleState
st 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 (forall a. [a] -> a
head [MultiParticleState]
mpsts)
          p1 :: Vec
p1 = MultiParticleState -> Vec
systemP (forall a. [a] -> a
last [MultiParticleState]
mpsts)
      in R
100 forall a. Num a => a -> a -> a
* Vec -> R
magnitude (Vec
p1 Vec -> Vec -> Vec
^-^ Vec
p0) 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 = forall a b. (RealFrac a, Integral b) => a -> b
floor (forall a. Floating a => a -> a -> a
logBase R
10 R
x) forall a. Num a => a -> a -> a
- Int
n forall a. Num a => a -> a -> a
+ Int
1
                  toInt :: R -> Int
                  toInt :: R -> Int
toInt = forall a b. (RealFrac a, Integral b) => a -> b
round
              in (Float
10forall a b. (Fractional a, Integral b) => a -> b -> a
^^Int
expon forall a. Num a => a -> a -> a
*) forall a b. (a -> b) -> a -> b
$ forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ R -> Int
toInt (R
10forall a b. (Fractional a, Integral b) => a -> b -> a
^^(-Int
expon) forall a. Num a => a -> a -> a
* R
x)

data Justification = LJ | RJ deriving Int -> Justification -> ShowS
[Justification] -> ShowS
Justification -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Justification] -> ShowS
$cshowList :: [Justification] -> ShowS
show :: Justification -> String
$cshow :: Justification -> String
showsPrec :: Int -> Justification -> ShowS
$cshowsPrec :: Int -> 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 = forall a. Show a => a -> String
show p
x in (String
str, forall (t :: * -> *) a. Foldable t => t a -> Int
length String
str)
              pairss :: [[(String, Int)]]
pairss = forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map forall {p}. Show p => p -> (String, Int)
pairWithLength) [[a]]
xss
              maxLength :: Int
maxLength = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (forall a b. (a -> b) -> [a] -> [b]
map forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map 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 forall a. [a] -> [a] -> [a]
++ forall a. Int -> a -> [a]
replicate (Int
maxLength forall a. Num a => a -> a -> a
+ Int
1 forall a. Num a => a -> a -> a
- Int
len) Char
' '
                      Justification
RJ -> forall a. Int -> a -> [a]
replicate (Int
maxLength forall a. Num a => a -> a -> a
+ Int
1 forall a. Num a => a -> a -> a
- Int
len) Char
' ' forall a. [a] -> [a] -> [a]
++ String
str
              showLine :: t (String, Int) -> String
showLine t (String, Int)
pairs = forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (String, Int) -> String
showPair t (String, Int)
pairs forall a. [a] -> [a] -> [a]
++ String
"\n"
          in forall a. [a] -> [a]
init forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap 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
    = forall a. Justification -> [[a]] -> Table a
Table Justification
LJ [[Int -> R -> Float
sigFigs Int
2 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 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 (forall a. [a] -> a
head [MultiParticleState]
mpsts)
          ke1 :: R
ke1 = MultiParticleState -> R
systemKE (forall a. [a] -> a
last [MultiParticleState]
mpsts)
      in R
100 forall a. Num a => a -> a -> a
* (R
ke1 forall a. Num a => a -> a -> a
- R
ke0) forall a. Fractional a => a -> a -> a
/ R
ke0

tenths :: R -> Float
tenths :: R -> Float
tenths = let toInt :: R -> Int
             toInt :: R -> Int
toInt = forall a b. (RealFrac a, Integral b) => a -> b
round
         in (forall a. Fractional a => a -> a -> a
/ Float
10) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (Integral a, Num b) => a -> b
fromIntegral forall b c a. (b -> c) -> (a -> b) -> a -> c
. R -> Int
toInt forall b c a. (b -> c) -> (a -> b) -> a -> c
. (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
    = forall a. Justification -> [[a]] -> Table a
Table Justification
RJ [[R -> Float
tenths 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 = forall (t :: * -> *) a. Foldable t => t a -> Int
length forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (a -> Bool) -> [a] -> [a]
takeWhile MultiParticleState -> Bool
inContact forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Bool -> Bool
not 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 forall a b. (a -> b) -> a -> b
$ OneBodyForce
posVec ([ParticleState]
sts forall a. [a] -> Int -> a
!! Int
0) Vec -> Vec -> Vec
^-^ OneBodyForce
posVec ([ParticleState]
sts forall a. [a] -> Int -> a
!! Int
1)
      in R
r forall a. Ord a => a -> a -> Bool
< R
2 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
    = 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 = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> [a] -> [b]
map MultiParticleState -> R
separation

separation :: MultiParticleState -> R
separation :: MultiParticleState -> R
separation (MPS [ParticleState]
sts)
    = Vec -> R
magnitude forall a b. (a -> b) -> a -> b
$ OneBodyForce
posVec ([ParticleState]
sts forall a. [a] -> Int -> a
!! Int
0) Vec -> Vec -> Vec
^-^ OneBodyForce
posVec ([ParticleState]
sts forall a. [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
    = forall a. Justification -> [[a]] -> Table a
Table Justification
RJ [[R -> Float
tenths forall a b. (a -> b) -> a -> b
$ (R
100forall a. Num a => a -> a -> a
*) 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))] forall a. [a] -> [a] -> [a]
++
      [Int -> Int -> TwoBodyForce -> Force
InternalForce Int
n (Int
nforall 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 (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 :: R
mass     = R
0.8293e-3 forall a. Num a => a -> a -> a
* R
0.65 forall a. Fractional a => a -> a -> a
/ R
64
           , posVec :: Vec
posVec   = R
x R -> Vec -> Vec
*^ Vec
iHat Vec -> Vec -> Vec
^+^ R
y R -> Vec -> Vec
*^ Vec
jHat
           , velocity :: Vec
velocity = Vec
zeroV
           } | R
x <- [R
0.01, R
0.02 .. R
0.64],
           let y :: R
y = R
0.005 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* R
x forall a. Fractional a => a -> a -> a
/ R
0.65)]

stringInitialPluck :: MultiParticleState
stringInitialPluck :: MultiParticleState
stringInitialPluck = [ParticleState] -> MultiParticleState
MPS [ParticleState
defaultParticleState
             { mass :: R
mass     = R
0.8293e-3 forall a. Num a => a -> a -> a
* R
0.65 forall a. Fractional a => a -> a -> a
/ R
64
             , posVec :: Vec
posVec   = R
x R -> Vec -> Vec
*^ Vec
iHat Vec -> Vec -> Vec
^+^ R
y R -> Vec -> Vec
*^ Vec
jHat
             , velocity :: Vec
velocity = Vec
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 forall a. Ord a => a -> a -> Bool
<= R
0.51  = R
0.005 forall a. Fractional a => a -> a -> a
/ (R
0.51 forall a. Num a => a -> a -> a
- R
0.00) forall a. Num a => a -> a -> a
* (R
x forall a. Num a => a -> a -> a
- R
0.00)
          | Bool
otherwise  = R
0.005 forall a. Fractional a => a -> a -> a
/ (R
0.51 forall a. Num a => a -> a -> a
- R
0.65) forall a. Num a => a -> a -> a
* (R
x forall a. Num a => a -> a -> a
- R
0.65)

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

mpsVel :: MultiParticleState -> IO ()
mpsVel :: MultiParticleState -> IO ()
mpsVel = 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 forall a. Ord a => a -> a -> Bool
>= R
re
         then Vec
zeroV
         else (-R
b) R -> Vec -> Vec
*^ Vec
v21