{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}

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

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

module LPFPCore.Mechanics1D where

import LPFPCore.Newton2 ( fAir )

import LPFPCore.SimpleVec ( R )

type Time     = R
type TimeStep = R
type Mass     = R
type Position = R
type Velocity = R
type Force    = R

type State1D = (Time,Position,Velocity)

newtonSecond1D :: Mass
               -> [State1D -> Force]  -- force funcs
               -> State1D             -- current state
               -> (R,R,R)             -- deriv of state
newtonSecond1D :: Force
-> [(Force, Force, Force) -> Force]
-> (Force, Force, Force)
-> (Force, Force, Force)
newtonSecond1D Force
m [(Force, Force, Force) -> Force]
fs (Force
t,Force
x0,Force
v0)
    = let fNet :: Force
fNet = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [(Force, Force, Force) -> Force
f (Force
t,Force
x0,Force
v0) | (Force, Force, Force) -> Force
f <- [(Force, Force, Force) -> Force]
fs]
          acc :: Force
acc = Force
fNet forall a. Fractional a => a -> a -> a
/ Force
m
      in (Force
1,Force
v0,Force
acc)

euler1D :: R                     -- time step dt
        -> (State1D -> (R,R,R))  -- differential equation
        -> State1D -> State1D    -- state-update function
euler1D :: Force
-> ((Force, Force, Force) -> (Force, Force, Force))
-> (Force, Force, Force)
-> (Force, Force, Force)
euler1D Force
dt (Force, Force, Force) -> (Force, Force, Force)
deriv (Force
t0,Force
x0,Force
v0)
    = let (Force
_, Force
_, Force
dvdt) = (Force, Force, Force) -> (Force, Force, Force)
deriv (Force
t0,Force
x0,Force
v0)
          t1 :: Force
t1 = Force
t0 forall a. Num a => a -> a -> a
+ Force
dt
          x1 :: Force
x1 = Force
x0 forall a. Num a => a -> a -> a
+ Force
v0 forall a. Num a => a -> a -> a
* Force
dt
          v1 :: Force
v1 = Force
v0 forall a. Num a => a -> a -> a
+ Force
dvdt forall a. Num a => a -> a -> a
* Force
dt
      in (Force
t1,Force
x1,Force
v1)

updateTXV :: R                   -- time interval dt
          -> Mass
          -> [State1D -> Force]  -- list of force funcs
          -> State1D -> State1D  -- state-update function
updateTXV :: Force
-> Force
-> [(Force, Force, Force) -> Force]
-> (Force, Force, Force)
-> (Force, Force, Force)
updateTXV Force
dt Force
m [(Force, Force, Force) -> Force]
fs = Force
-> ((Force, Force, Force) -> (Force, Force, Force))
-> (Force, Force, Force)
-> (Force, Force, Force)
euler1D Force
dt (Force
-> [(Force, Force, Force) -> Force]
-> (Force, Force, Force)
-> (Force, Force, Force)
newtonSecond1D Force
m [(Force, Force, Force) -> Force]
fs)

statesTXV :: R                   -- time step
          -> Mass
          -> State1D             -- initial state
          -> [State1D -> Force]  -- list of force funcs
          -> [State1D]           -- infinite list of states
statesTXV :: Force
-> Force
-> (Force, Force, Force)
-> [(Force, Force, Force) -> Force]
-> [(Force, Force, Force)]
statesTXV Force
dt Force
m (Force, Force, Force)
txv0 [(Force, Force, Force) -> Force]
fs = forall a. (a -> a) -> a -> [a]
iterate (Force
-> Force
-> [(Force, Force, Force) -> Force]
-> (Force, Force, Force)
-> (Force, Force, Force)
updateTXV Force
dt Force
m [(Force, Force, Force) -> Force]
fs) (Force, Force, Force)
txv0

-- assume that dt is the same between adjacent pairs
velocity1D :: [State1D]           -- infinite list
           -> Time -> Velocity    -- velocity function
velocity1D :: [(Force, Force, Force)] -> Force -> Force
velocity1D [(Force, Force, Force)]
sts Force
t
    = let (Force
t0,Force
_,Force
_) = [(Force, Force, Force)]
sts forall a. [a] -> Int -> a
!! Int
0
          (Force
t1,Force
_,Force
_) = [(Force, Force, Force)]
sts forall a. [a] -> Int -> a
!! Int
1
          dt :: Force
dt = Force
t1 forall a. Num a => a -> a -> a
- Force
t0
          numSteps :: Int
numSteps = forall a. Num a => a -> a
abs forall a b. (a -> b) -> a -> b
$ forall a b. (RealFrac a, Integral b) => a -> b
round (Force
t forall a. Fractional a => a -> a -> a
/ Force
dt)
          (Force
_,Force
_,Force
v0) = [(Force, Force, Force)]
sts forall a. [a] -> Int -> a
!! Int
numSteps
      in Force
v0

velocityFtxv :: R                   -- time step
             -> Mass
             -> State1D             -- initial state
             -> [State1D -> Force]  -- list of force funcs
             -> Time -> Velocity    -- velocity function
velocityFtxv :: Force
-> Force
-> (Force, Force, Force)
-> [(Force, Force, Force) -> Force]
-> Force
-> Force
velocityFtxv Force
dt Force
m (Force, Force, Force)
txv0 [(Force, Force, Force) -> Force]
fs = [(Force, Force, Force)] -> Force -> Force
velocity1D (Force
-> Force
-> (Force, Force, Force)
-> [(Force, Force, Force) -> Force]
-> [(Force, Force, Force)]
statesTXV Force
dt Force
m (Force, Force, Force)
txv0 [(Force, Force, Force) -> Force]
fs)

-- assume that dt is the same between adjacent pairs
position1D :: [State1D]           -- infinite list
           -> Time -> Position    -- position function
position1D :: [(Force, Force, Force)] -> Force -> Force
position1D [(Force, Force, Force)]
sts Force
t
    = let (Force
t0,Force
_,Force
_) = [(Force, Force, Force)]
sts forall a. [a] -> Int -> a
!! Int
0
          (Force
t1,Force
_,Force
_) = [(Force, Force, Force)]
sts forall a. [a] -> Int -> a
!! Int
1
          dt :: Force
dt = Force
t1 forall a. Num a => a -> a -> a
- Force
t0
          numSteps :: Int
numSteps = forall a. Num a => a -> a
abs forall a b. (a -> b) -> a -> b
$ forall a b. (RealFrac a, Integral b) => a -> b
round (Force
t forall a. Fractional a => a -> a -> a
/ Force
dt)
          (Force
_,Force
x0,Force
_) = [(Force, Force, Force)]
sts forall a. [a] -> Int -> a
!! Int
numSteps
      in Force
x0

positionFtxv :: R                   -- time step
             -> Mass
             -> State1D             -- initial state
             -> [State1D -> Force]  -- list of force funcs
             -> Time -> Position    -- position function
positionFtxv :: Force
-> Force
-> (Force, Force, Force)
-> [(Force, Force, Force) -> Force]
-> Force
-> Force
positionFtxv Force
dt Force
m (Force, Force, Force)
txv0 [(Force, Force, Force) -> Force]
fs = [(Force, Force, Force)] -> Force -> Force
position1D (Force
-> Force
-> (Force, Force, Force)
-> [(Force, Force, Force) -> Force]
-> [(Force, Force, Force)]
statesTXV Force
dt Force
m (Force, Force, Force)
txv0 [(Force, Force, Force) -> Force]
fs)

springForce :: R -> State1D -> Force
springForce :: Force -> (Force, Force, Force) -> Force
springForce Force
k (Force
_,Force
x0,Force
_) = -Force
k forall a. Num a => a -> a -> a
* Force
x0

dampedHOForces :: [State1D -> Force]
dampedHOForces :: [(Force, Force, Force) -> Force]
dampedHOForces = [Force -> (Force, Force, Force) -> Force
springForce Force
0.8
                 ,\(Force
_,Force
_,Force
v0) -> Force -> Force -> Force -> Force -> Force
fAir Force
2 Force
1.225 (forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Force
0.02forall a. Floating a => a -> a -> a
**Force
2) Force
v0
                 ,\(Force, Force, Force)
_ -> -Force
0.0027 forall a. Num a => a -> a -> a
* Force
9.80665
                 ]

dampedHOStates :: [State1D]
dampedHOStates :: [(Force, Force, Force)]
dampedHOStates = Force
-> Force
-> (Force, Force, Force)
-> [(Force, Force, Force) -> Force]
-> [(Force, Force, Force)]
statesTXV Force
0.001 Force
0.0027 (Force
0.0,Force
0.1,Force
0.0) [(Force, Force, Force) -> Force]
dampedHOForces

pingpongPosition :: Time -> Velocity
pingpongPosition :: Force -> Force
pingpongPosition = Force
-> Force
-> (Force, Force, Force)
-> [(Force, Force, Force) -> Force]
-> Force
-> Force
positionFtxv Force
0.001 Force
0.0027 (Force
0,Force
0.1,Force
0) [(Force, Force, Force) -> Force]
dampedHOForces

pingpongVelocity :: Time -> Velocity
pingpongVelocity :: Force -> Force
pingpongVelocity = Force
-> Force
-> (Force, Force, Force)
-> [(Force, Force, Force) -> Force]
-> Force
-> Force
velocityFtxv Force
0.001 Force
0.0027 (Force
0,Force
0.1,Force
0) [(Force, Force, Force) -> Force]
dampedHOForces

eulerCromer1D :: R                     -- time step dt
              -> (State1D -> (R,R,R))  -- differential equation
              -> State1D -> State1D    -- state-update function
eulerCromer1D :: Force
-> ((Force, Force, Force) -> (Force, Force, Force))
-> (Force, Force, Force)
-> (Force, Force, Force)
eulerCromer1D Force
dt (Force, Force, Force) -> (Force, Force, Force)
deriv (Force
t0,Force
x0,Force
v0)
    = let (Force
_, Force
_, Force
dvdt) = (Force, Force, Force) -> (Force, Force, Force)
deriv (Force
t0,Force
x0,Force
v0)
          t1 :: Force
t1 = Force
t0 forall a. Num a => a -> a -> a
+ Force
dt
          x1 :: Force
x1 = Force
x0 forall a. Num a => a -> a -> a
+ Force
v1 forall a. Num a => a -> a -> a
* Force
dt
          v1 :: Force
v1 = Force
v0 forall a. Num a => a -> a -> a
+ Force
dvdt forall a. Num a => a -> a -> a
* Force
dt
      in (Force
t1,Force
x1,Force
v1)

updateTXVEC :: R                   -- time interval dt
            -> Mass
            -> [State1D -> Force]  -- list of force funcs
            -> State1D -> State1D  -- state-update function
updateTXVEC :: Force
-> Force
-> [(Force, Force, Force) -> Force]
-> (Force, Force, Force)
-> (Force, Force, Force)
updateTXVEC Force
dt Force
m [(Force, Force, Force) -> Force]
fs = Force
-> ((Force, Force, Force) -> (Force, Force, Force))
-> (Force, Force, Force)
-> (Force, Force, Force)
eulerCromer1D Force
dt (Force
-> [(Force, Force, Force) -> Force]
-> (Force, Force, Force)
-> (Force, Force, Force)
newtonSecond1D Force
m [(Force, Force, Force) -> Force]
fs)

-- | An update function takes a state as input and returns an updated state as output.
type UpdateFunction s = s -> s

-- | A differential equation takes a state as input and returns as output the rate at which
--   the state is changing.
type DifferentialEquation s ds = s -> ds

-- | A numerical method turns a differential equation into a state-update function.
type NumericalMethod s ds = DifferentialEquation s ds -> UpdateFunction s

-- | Given a numerical method, a differential equation, and an initial state,
--   return a list of states.
solver :: NumericalMethod s ds -> DifferentialEquation s ds -> s -> [s]
solver :: forall s ds.
NumericalMethod s ds -> DifferentialEquation s ds -> s -> [s]
solver NumericalMethod s ds
method = forall a. (a -> a) -> a -> [a]
iterate forall b c a. (b -> c) -> (a -> b) -> a -> c
. NumericalMethod s ds
method

-- | A real vector space allows vector addition and scalar multiplication by reals.
class RealVectorSpace ds where
      (+++) :: ds -> ds -> ds
      scale :: R -> ds -> ds

-- | A triple of real numbers is a real vector space.
instance RealVectorSpace (R,R,R) where
    (Force
dtdt0, Force
dxdt0, Force
dvdt0) +++ :: (Force, Force, Force)
-> (Force, Force, Force) -> (Force, Force, Force)
+++ (Force
dtdt1, Force
dxdt1, Force
dvdt1)
        = (Force
dtdt0 forall a. Num a => a -> a -> a
+ Force
dtdt1, Force
dxdt0 forall a. Num a => a -> a -> a
+ Force
dxdt1, Force
dvdt0 forall a. Num a => a -> a -> a
+ Force
dvdt1)
    scale :: Force -> (Force, Force, Force) -> (Force, Force, Force)
scale Force
w (Force
dtdt0, Force
dxdt0, Force
dvdt0) = (Force
w forall a. Num a => a -> a -> a
* Force
dtdt0, Force
w forall a. Num a => a -> a -> a
* Force
dxdt0, Force
w forall a. Num a => a -> a -> a
* Force
dvdt0)

-- | A type class that expresses a relationship between a state space
--   and a time-derivative-state space.
class RealVectorSpace ds => Diff s ds where
    shift :: R -> ds -> s -> s

-- | A triple of real numbers can serve as the time derivative of a 'State1D'.
instance Diff State1D (R,R,R) where
    shift :: Force
-> (Force, Force, Force)
-> (Force, Force, Force)
-> (Force, Force, Force)
shift Force
dt (Force
dtdt,Force
dxdt,Force
dvdt) (Force
t,Force
x,Force
v)
        = (Force
t forall a. Num a => a -> a -> a
+ Force
dtdt forall a. Num a => a -> a -> a
* Force
dt, Force
x forall a. Num a => a -> a -> a
+ Force
dxdt forall a. Num a => a -> a -> a
* Force
dt, Force
v forall a. Num a => a -> a -> a
+ Force
dvdt forall a. Num a => a -> a -> a
* Force
dt)

-- | Given a step size, return the numerical method that uses the Euler
--   method with that step size.
euler :: Diff s ds => R -> (s -> ds) -> s -> s
euler :: forall s ds. Diff s ds => Force -> (s -> ds) -> s -> s
euler Force
dt s -> ds
deriv s
st0 = forall s ds. Diff s ds => Force -> ds -> s -> s
shift Force
dt (s -> ds
deriv s
st0) s
st0

-- | Given a step size, return the numerical method that uses the 4th order Runge Kutta
--   method with that step size.
rungeKutta4 :: Diff s ds => R -> (s -> ds) -> s -> s
rungeKutta4 :: forall s ds. Diff s ds => Force -> (s -> ds) -> s -> s
rungeKutta4 Force
dt s -> ds
deriv s
st0
    = let m0 :: ds
m0 = s -> ds
deriv                  s
st0
          m1 :: ds
m1 = s -> ds
deriv (forall s ds. Diff s ds => Force -> ds -> s -> s
shift (Force
dtforall a. Fractional a => a -> a -> a
/Force
2) ds
m0 s
st0)
          m2 :: ds
m2 = s -> ds
deriv (forall s ds. Diff s ds => Force -> ds -> s -> s
shift (Force
dtforall a. Fractional a => a -> a -> a
/Force
2) ds
m1 s
st0)
          m3 :: ds
m3 = s -> ds
deriv (forall s ds. Diff s ds => Force -> ds -> s -> s
shift  Force
dt    ds
m2 s
st0)
      in forall s ds. Diff s ds => Force -> ds -> s -> s
shift (Force
dtforall a. Fractional a => a -> a -> a
/Force
6) (ds
m0 forall ds. RealVectorSpace ds => ds -> ds -> ds
+++ ds
m1 forall ds. RealVectorSpace ds => ds -> ds -> ds
+++ ds
m1 forall ds. RealVectorSpace ds => ds -> ds -> ds
+++ ds
m2 forall ds. RealVectorSpace ds => ds -> ds -> ds
+++ ds
m2 forall ds. RealVectorSpace ds => ds -> ds -> ds
+++ ds
m3) s
st0

exponential :: DifferentialEquation (R,R,R) (R,R,R)
exponential :: (Force, Force, Force) -> (Force, Force, Force)
exponential (Force
_,Force
x0,Force
v0) = (Force
1,Force
v0,Force
x0)

update2 :: (R,R,R)  -- starting state
        -> (R,R,R)  -- ending state
update2 :: (Force, Force, Force) -> (Force, Force, Force)
update2 = forall a. HasCallStack => a
undefined

earthGravity :: Mass -> State1D -> Force
earthGravity :: Force -> (Force, Force, Force) -> Force
earthGravity Force
m (Force, Force, Force)
_ = let g :: Force
g = Force
9.80665
                   in -Force
m forall a. Num a => a -> a -> a
* Force
g

type MState = (Time,Mass,Position,Velocity)

earthGravity2 :: MState -> Force
earthGravity2 :: MState -> Force
earthGravity2 (Force
_,Force
m,Force
_,Force
_) = let g :: Force
g = Force
9.80665
                          in -Force
m forall a. Num a => a -> a -> a
* Force
g

positionFtxv2 :: R                  -- time step
              -> MState             -- initial state
              -> [MState -> Force]  -- list of force funcs
              -> Time -> Position   -- position function
positionFtxv2 :: Force -> MState -> [MState -> Force] -> Force -> Force
positionFtxv2 = forall a. HasCallStack => a
undefined

statesTXV2 :: R                 -- time step
          -> MState             -- initial state
          -> [MState -> Force]  -- list of force funcs
          -> [MState]           -- infinite list of states
statesTXV2 :: Force -> MState -> [MState -> Force] -> [MState]
statesTXV2 = forall a. HasCallStack => a
undefined

updateTXV2 :: R                  -- dt for stepping
           -> [MState -> Force]  -- list of force funcs
           -> MState             -- current state
           -> MState             -- new state
updateTXV2 :: Force -> [MState -> Force] -> MState -> MState
updateTXV2 = forall a. HasCallStack => a
undefined

instance RealVectorSpace (R,R) where
    (Force
dtdt0, Force
dvdt0) +++ :: (Force, Force) -> (Force, Force) -> (Force, Force)
+++ (Force
dtdt1, Force
dvdt1) = (Force
dtdt0 forall a. Num a => a -> a -> a
+ Force
dtdt1, Force
dvdt0 forall a. Num a => a -> a -> a
+ Force
dvdt1)
    scale :: Force -> (Force, Force) -> (Force, Force)
scale Force
w (Force
dtdt0, Force
dvdt0) = (Force
w forall a. Num a => a -> a -> a
* Force
dtdt0, Force
w forall a. Num a => a -> a -> a
* Force
dvdt0)

instance Diff (Time,Velocity) (R,R) where
    shift :: Force -> (Force, Force) -> (Force, Force) -> (Force, Force)
shift Force
dt (Force
dtdt,Force
dvdt) (Force
t,Force
v)
        = (Force
t forall a. Num a => a -> a -> a
+ Force
dtdt forall a. Num a => a -> a -> a
* Force
dt, Force
v forall a. Num a => a -> a -> a
+ Force
dvdt forall a. Num a => a -> a -> a
* Force
dt)

updateTV' :: R                           -- dt for stepping
          -> Mass
          -> [(Time,Velocity) -> Force]  -- list of force funcs
          -> (Time,Velocity)             -- current state
          -> (Time,Velocity)             -- new state
updateTV' :: Force
-> Force
-> [(Force, Force) -> Force]
-> (Force, Force)
-> (Force, Force)
updateTV' = forall a. HasCallStack => a
undefined

forces :: R -> [State1D -> R]
forces :: Force -> [(Force, Force, Force) -> Force]
forces Force
mu = [\(Force
_t,Force
x,Force
_v) -> forall a. HasCallStack => a
undefined Force
x
            ,\(Force
_t,Force
x, Force
v) -> forall a. HasCallStack => a
undefined Force
mu Force
x Force
v]

vdp :: R -> [(R,R)]
vdp :: Force -> [(Force, Force)]
vdp Force
mu = forall a b. (a -> b) -> [a] -> [b]
map (\(Force
_,Force
x,Force
v) -> (Force
x,Force
v)) forall a b. (a -> b) -> a -> b
$ forall a. Int -> [a] -> [a]
take Int
10000 forall a b. (a -> b) -> a -> b
$
         forall s ds.
NumericalMethod s ds -> DifferentialEquation s ds -> s -> [s]
solver (forall s ds. Diff s ds => Force -> (s -> ds) -> s -> s
rungeKutta4 Force
0.01) (Force
-> [(Force, Force, Force) -> Force]
-> (Force, Force, Force)
-> (Force, Force, Force)
newtonSecond1D Force
1 forall a b. (a -> b) -> a -> b
$ Force -> [(Force, Force, Force) -> Force]
forces Force
mu) (Force
0,Force
2,Force
0)