{-# OPTIONS -Wall #-}
{-# LANGUAGE MultiParamTypeClasses #-}

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

Code from chapters 16, 17, and 18 of the book Learn Physics with Functional Programming
-}

module LPFPCore.Mechanics3D where

import LPFPCore.SimpleVec
    ( R, Vec, PosVec, (^+^), (^-^), (*^), (^*), (^/), (<.>), (><)
    , vec, sumV, magnitude, zeroV, xComp, yComp, zComp, iHat, jHat, kHat)
import LPFPCore.Mechanics1D
    ( RealVectorSpace(..), Diff(..), NumericalMethod
    , Time, TimeStep, rungeKutta4, solver )

-- | Data type for the state of a single particle in three-dimensional space.
data ParticleState = ParticleState { ParticleState -> R
mass     :: R
                                   , ParticleState -> R
charge   :: R
                                   , ParticleState -> R
time     :: R
                                   , ParticleState -> Vec
posVec   :: Vec
                                   , ParticleState -> Vec
velocity :: Vec }
                     deriving Int -> ParticleState -> ShowS
[ParticleState] -> ShowS
ParticleState -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [ParticleState] -> ShowS
$cshowList :: [ParticleState] -> ShowS
show :: ParticleState -> String
$cshow :: ParticleState -> String
showsPrec :: Int -> ParticleState -> ShowS
$cshowsPrec :: Int -> ParticleState -> ShowS
Show

-- | A default particle state.
defaultParticleState :: ParticleState
defaultParticleState :: ParticleState
defaultParticleState = ParticleState { mass :: R
mass     = R
1
                                     , charge :: R
charge   = R
0
                                     , time :: R
time     = R
0
                                     , posVec :: Vec
posVec   = Vec
zeroV
                                     , velocity :: Vec
velocity = Vec
zeroV }

rockState :: ParticleState
rockState :: ParticleState
rockState
    = ParticleState
defaultParticleState { mass :: R
mass     = R
2                        -- kg
                           , velocity :: Vec
velocity = R
3 R -> Vec -> Vec
*^ Vec
iHat Vec -> Vec -> Vec
^+^ R
4 R -> Vec -> Vec
*^ Vec
kHat  -- m/s
                           }

-- | Data type for a one-body force.
type OneBodyForce = ParticleState -> Vec

-- | Data type for the time-derivative of a particle state.
data DParticleState = DParticleState { DParticleState -> R
dmdt :: R
                                     , DParticleState -> R
dqdt :: R
                                     , DParticleState -> R
dtdt :: R
                                     , DParticleState -> Vec
drdt :: Vec
                                     , DParticleState -> Vec
dvdt :: Vec }
                      deriving Int -> DParticleState -> ShowS
[DParticleState] -> ShowS
DParticleState -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [DParticleState] -> ShowS
$cshowList :: [DParticleState] -> ShowS
show :: DParticleState -> String
$cshow :: DParticleState -> String
showsPrec :: Int -> DParticleState -> ShowS
$cshowsPrec :: Int -> DParticleState -> ShowS
Show

-- | Given a list of forces, return a differential equation
--   based on Newton's second law.
newtonSecondPS :: [OneBodyForce]
               -> ParticleState -> DParticleState  -- ^ a differential equation
newtonSecondPS :: [ParticleState -> Vec] -> ParticleState -> DParticleState
newtonSecondPS [ParticleState -> Vec]
fs ParticleState
st
    = let fNet :: Vec
fNet = [Vec] -> Vec
sumV [ParticleState -> Vec
f ParticleState
st | ParticleState -> Vec
f <- [ParticleState -> Vec]
fs]
          m :: R
m = ParticleState -> R
mass ParticleState
st
          v :: Vec
v = ParticleState -> Vec
velocity ParticleState
st
          acc :: Vec
acc = Vec
fNet Vec -> R -> Vec
^/ R
m
      in DParticleState { dmdt :: R
dmdt = R
0    -- dm/dt
                        , dqdt :: R
dqdt = R
0    -- dq/dt
                        , dtdt :: R
dtdt = R
1    -- dt/dt
                        , drdt :: Vec
drdt = Vec
v    -- dr/dt
                        , dvdt :: Vec
dvdt = Vec
acc  -- dv/dt
                        }

-- | The force of gravity near Earth's surface.
--   The z direction is toward the sky.
--   Assumes SI units.
earthSurfaceGravity :: OneBodyForce
earthSurfaceGravity :: ParticleState -> Vec
earthSurfaceGravity ParticleState
st
    = let g :: R
g = R
9.80665  -- m/s^2
      in (-ParticleState -> R
mass ParticleState
st forall a. Num a => a -> a -> a
* R
g) R -> Vec -> Vec
*^ Vec
kHat

-- | The force of the Sun's gravity on an object.
--   The origin is at center of the Sun.
--   Assumes SI units.
sunGravity :: OneBodyForce
sunGravity :: ParticleState -> Vec
sunGravity (ParticleState R
m R
_q R
_t Vec
r Vec
_v)
    = let bigG :: R
bigG = R
6.67408e-11  -- N m^2/kg^2
          sunMass :: R
sunMass = R
1.98848e30  -- kg
      in (-R
bigG forall a. Num a => a -> a -> a
* R
sunMass forall a. Num a => a -> a -> a
* R
m) R -> Vec -> Vec
*^ Vec
r Vec -> R -> Vec
^/ Vec -> R
magnitude Vec
r forall a. Floating a => a -> a -> a
** R
3

-- | The force of air resistance on an object.
airResistance :: R  -- ^ drag coefficient
              -> R  -- ^ air density
              -> R  -- ^ cross-sectional area of object
              -> OneBodyForce
airResistance :: R -> R -> R -> ParticleState -> Vec
airResistance R
drag R
rho R
area (ParticleState R
_m R
_q R
_t Vec
_r Vec
v)
    = (-R
0.5 forall a. Num a => a -> a -> a
* R
drag forall a. Num a => a -> a -> a
* R
rho forall a. Num a => a -> a -> a
* R
area forall a. Num a => a -> a -> a
* Vec -> R
magnitude Vec
v) R -> Vec -> Vec
*^ Vec
v

-- | The force of wind on an object.
windForce :: Vec  -- ^ wind velocity
          -> R    -- ^ drag coefficient
          -> R    -- ^ air density
          -> R    -- ^ cross-sectional area of object
          -> OneBodyForce
windForce :: Vec -> R -> R -> R -> ParticleState -> Vec
windForce Vec
vWind R
drag R
rho R
area (ParticleState R
_m R
_q R
_t Vec
_r Vec
v)
    = let vRel :: Vec
vRel = Vec
v Vec -> Vec -> Vec
^-^ Vec
vWind
      in (-R
0.5 forall a. Num a => a -> a -> a
* R
drag forall a. Num a => a -> a -> a
* R
rho forall a. Num a => a -> a -> a
* R
area forall a. Num a => a -> a -> a
* Vec -> R
magnitude Vec
vRel) R -> Vec -> Vec
*^ Vec
vRel

-- | The force of uniform electric and magnetic fields on an object.
uniformLorentzForce :: Vec  -- ^ E
                    -> Vec  -- ^ B
                    -> OneBodyForce
uniformLorentzForce :: Vec -> Vec -> ParticleState -> Vec
uniformLorentzForce Vec
vE Vec
vB (ParticleState R
_m R
q R
_t Vec
_r Vec
v)
    = R
q R -> Vec -> Vec
*^ (Vec
vE Vec -> Vec -> Vec
^+^ Vec
v Vec -> Vec -> Vec
>< Vec
vB)

-- | Euler-Cromer method for the 'ParticleState' data type.
eulerCromerPS :: TimeStep        -- dt for stepping
              -> NumericalMethod ParticleState DParticleState
eulerCromerPS :: R -> NumericalMethod ParticleState DParticleState
eulerCromerPS R
dt ParticleState -> DParticleState
deriv ParticleState
st
    = let t :: R
t   = ParticleState -> R
time     ParticleState
st
          r :: Vec
r   = ParticleState -> Vec
posVec   ParticleState
st
          v :: Vec
v   = ParticleState -> Vec
velocity ParticleState
st
          dst :: DParticleState
dst = ParticleState -> DParticleState
deriv ParticleState
st
          acc :: Vec
acc = DParticleState -> Vec
dvdt DParticleState
dst
          v' :: Vec
v'  = Vec
v Vec -> Vec -> Vec
^+^ Vec
acc Vec -> R -> Vec
^* R
dt
      in ParticleState
st { time :: R
time     = R
t  forall a. Num a => a -> a -> a
+         R
dt
            , posVec :: Vec
posVec   = Vec
r Vec -> Vec -> Vec
^+^ Vec
v'  Vec -> R -> Vec
^* R
dt
            , velocity :: Vec
velocity = Vec
v Vec -> Vec -> Vec
^+^ Vec
acc Vec -> R -> Vec
^* R
dt
            }

instance RealVectorSpace DParticleState where
    DParticleState
dst1 +++ :: DParticleState -> DParticleState -> DParticleState
+++ DParticleState
dst2
        = DParticleState { dmdt :: R
dmdt = DParticleState -> R
dmdt DParticleState
dst1  forall a. Num a => a -> a -> a
+  DParticleState -> R
dmdt DParticleState
dst2
                         , dqdt :: R
dqdt = DParticleState -> R
dqdt DParticleState
dst1  forall a. Num a => a -> a -> a
+  DParticleState -> R
dqdt DParticleState
dst2
                         , dtdt :: R
dtdt = DParticleState -> R
dtdt DParticleState
dst1  forall a. Num a => a -> a -> a
+  DParticleState -> R
dtdt DParticleState
dst2
                         , drdt :: Vec
drdt = DParticleState -> Vec
drdt DParticleState
dst1 Vec -> Vec -> Vec
^+^ DParticleState -> Vec
drdt DParticleState
dst2
                         , dvdt :: Vec
dvdt = DParticleState -> Vec
dvdt DParticleState
dst1 Vec -> Vec -> Vec
^+^ DParticleState -> Vec
dvdt DParticleState
dst2
                         }
    scale :: R -> DParticleState -> DParticleState
scale R
w DParticleState
dst
        = DParticleState { dmdt :: R
dmdt = R
w forall a. Num a => a -> a -> a
*  DParticleState -> R
dmdt DParticleState
dst
                         , dqdt :: R
dqdt = R
w forall a. Num a => a -> a -> a
*  DParticleState -> R
dqdt DParticleState
dst
                         , dtdt :: R
dtdt = R
w forall a. Num a => a -> a -> a
*  DParticleState -> R
dtdt DParticleState
dst
                         , drdt :: Vec
drdt = R
w R -> Vec -> Vec
*^ DParticleState -> Vec
drdt DParticleState
dst
                         , dvdt :: Vec
dvdt = R
w R -> Vec -> Vec
*^ DParticleState -> Vec
dvdt DParticleState
dst
                         }

instance Diff ParticleState DParticleState where
    shift :: R -> DParticleState -> ParticleState -> ParticleState
shift R
dt DParticleState
dps (ParticleState R
m R
q R
t Vec
r Vec
v)
        = R -> R -> R -> Vec -> Vec -> ParticleState
ParticleState (R
m  forall a. Num a => a -> a -> a
+  DParticleState -> R
dmdt DParticleState
dps  forall a. Num a => a -> a -> a
* R
dt)
                        (R
q  forall a. Num a => a -> a -> a
+  DParticleState -> R
dqdt DParticleState
dps  forall a. Num a => a -> a -> a
* R
dt)
                        (R
t  forall a. Num a => a -> a -> a
+  DParticleState -> R
dtdt DParticleState
dps  forall a. Num a => a -> a -> a
* R
dt)
                        (Vec
r Vec -> Vec -> Vec
^+^ DParticleState -> Vec
drdt DParticleState
dps Vec -> R -> Vec
^* R
dt)
                        (Vec
v Vec -> Vec -> Vec
^+^ DParticleState -> Vec
dvdt DParticleState
dps Vec -> R -> Vec
^* R
dt)

-- | Given a numerical method,
--   a list of one-body forces, and an initial state,
--   return a list of states describing how the particle
--   evolves in time.
statesPS :: NumericalMethod ParticleState DParticleState  -- ^ numerical method
         -> [OneBodyForce]  -- ^ list of force funcs
         -> ParticleState -> [ParticleState]  -- ^ evolver
statesPS :: NumericalMethod ParticleState DParticleState
-> [ParticleState -> Vec] -> ParticleState -> [ParticleState]
statesPS NumericalMethod ParticleState DParticleState
method = forall a. (a -> a) -> a -> [a]
iterate forall b c a. (b -> c) -> (a -> b) -> a -> c
. NumericalMethod ParticleState DParticleState
method forall b c a. (b -> c) -> (a -> b) -> a -> c
. [ParticleState -> Vec] -> ParticleState -> DParticleState
newtonSecondPS

-- | Given a numerical method and a list of one-body forces,
--   return a state-update function.
updatePS :: NumericalMethod ParticleState DParticleState
         -> [OneBodyForce]
         -> ParticleState -> ParticleState
updatePS :: NumericalMethod ParticleState DParticleState
-> [ParticleState -> Vec] -> ParticleState -> ParticleState
updatePS NumericalMethod ParticleState DParticleState
method = NumericalMethod ParticleState DParticleState
method forall b c a. (b -> c) -> (a -> b) -> a -> c
. [ParticleState -> Vec] -> ParticleState -> DParticleState
newtonSecondPS

-- | Given a numerical method,
--   a list of one-body forces, and an initial state,
--   return a position function describing how the particle
--   evolves in time.
positionPS :: NumericalMethod ParticleState DParticleState
           -> [OneBodyForce]  -- ^ list of force funcs
           -> ParticleState   -- ^ initial state
           -> Time -> PosVec  -- ^ position function
positionPS :: NumericalMethod ParticleState DParticleState
-> [ParticleState -> Vec] -> ParticleState -> R -> Vec
positionPS NumericalMethod ParticleState DParticleState
method [ParticleState -> Vec]
fs ParticleState
st R
t
    = let states :: [ParticleState]
states = NumericalMethod ParticleState DParticleState
-> [ParticleState -> Vec] -> ParticleState -> [ParticleState]
statesPS NumericalMethod ParticleState DParticleState
method [ParticleState -> Vec]
fs ParticleState
st
          dt :: R
dt = ParticleState -> R
time ([ParticleState]
states forall a. [a] -> Int -> a
!! Int
1) forall a. Num a => a -> a -> a
- ParticleState -> R
time ([ParticleState]
states forall a. [a] -> Int -> a
!! Int
0)
          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 (R
t forall a. Fractional a => a -> a -> a
/ R
dt)
          st1 :: ParticleState
st1 = forall s ds.
NumericalMethod s ds -> DifferentialEquation s ds -> s -> [s]
solver NumericalMethod ParticleState DParticleState
method ([ParticleState -> Vec] -> ParticleState -> DParticleState
newtonSecondPS [ParticleState -> Vec]
fs) ParticleState
st forall a. [a] -> Int -> a
!! Int
numSteps
      in ParticleState -> Vec
posVec ParticleState
st1

class HasTime s where
    timeOf :: s -> Time

instance HasTime ParticleState where
    timeOf :: ParticleState -> R
timeOf = ParticleState -> R
time

constantForce :: Vec -> OneBodyForce
constantForce :: Vec -> ParticleState -> Vec
constantForce Vec
f = forall a. HasCallStack => a
undefined Vec
f

moonSurfaceGravity :: OneBodyForce
moonSurfaceGravity :: ParticleState -> Vec
moonSurfaceGravity = forall a. HasCallStack => a
undefined

earthGravity :: OneBodyForce
earthGravity :: ParticleState -> Vec
earthGravity = forall a. HasCallStack => a
undefined

tvyPair :: ParticleState -> (R,R)
tvyPair :: ParticleState -> (R, R)
tvyPair ParticleState
st = forall a. HasCallStack => a
undefined ParticleState
st

tvyPairs :: [ParticleState] -> [(R,R)]
tvyPairs :: [ParticleState] -> [(R, R)]
tvyPairs [ParticleState]
sts = forall a. HasCallStack => a
undefined [ParticleState]
sts

tle1yr :: ParticleState -> Bool
tle1yr :: ParticleState -> Bool
tle1yr ParticleState
st = forall a. HasCallStack => a
undefined ParticleState
st

stateFunc :: [ParticleState]
          -> Time -> ParticleState
stateFunc :: [ParticleState] -> R -> ParticleState
stateFunc [ParticleState]
sts R
t
    = let t0 :: t
t0 = forall a. HasCallStack => a
undefined [ParticleState]
sts
          t1 :: t
t1 = forall a. HasCallStack => a
undefined [ParticleState]
sts
          dt :: t
dt = forall a. HasCallStack => a
undefined forall {t}. t
t0 forall {t}. t
t1
          numSteps :: t
numSteps = forall a. HasCallStack => a
undefined R
t forall {t}. t
dt
      in forall a. HasCallStack => a
undefined [ParticleState]
sts forall {t}. t
numSteps

airResAtAltitude :: R  -- ^ drag coefficient
                 -> R  -- ^ air density at sea level
                 -> R  -- ^ cross-sectional area of object
                 -> OneBodyForce
airResAtAltitude :: R -> R -> R -> ParticleState -> Vec
airResAtAltitude R
drag R
rho0 R
area (ParticleState R
_m R
_q R
_t Vec
r Vec
v)
    = forall a. HasCallStack => a
undefined R
drag R
rho0 R
area Vec
r Vec
v

projectileRangeComparison :: R -> R -> (R,R,R)
projectileRangeComparison :: R -> R -> (R, R, R)
projectileRangeComparison R
v0 R
thetaDeg
    = let vx0 :: R
vx0 = R
v0 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos (R
thetaDeg forall a. Fractional a => a -> a -> a
/ R
180 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi)
          vz0 :: R
vz0 = R
v0 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin (R
thetaDeg forall a. Fractional a => a -> a -> a
/ R
180 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi)
          drag :: R
drag = R
1
          ballRadius :: R
ballRadius = R
0.05    -- meters
          area :: R
area = forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* R
ballRadiusforall a. Floating a => a -> a -> a
**R
2
          airDensity :: R
airDensity  =     R
1.225  -- kg/m^3 @ sea level
          leadDensity :: R
leadDensity = R
11342      -- kg/m^3
          m :: R
m = R
leadDensity forall a. Num a => a -> a -> a
* R
4 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* R
ballRadiusforall a. Floating a => a -> a -> a
**R
3 forall a. Fractional a => a -> a -> a
/ R
3
          stateInitial :: t
stateInitial = forall a. HasCallStack => a
undefined R
m R
vx0 R
vz0
          aboveSeaLevel :: ParticleState -> Bool
          aboveSeaLevel :: ParticleState -> Bool
aboveSeaLevel ParticleState
st = Vec -> R
zComp (ParticleState -> Vec
posVec ParticleState
st) forall a. Ord a => a -> a -> Bool
>= R
0
          range :: [ParticleState] -> R
          range :: [ParticleState] -> R
range = Vec -> R
xComp forall b c a. (b -> c) -> (a -> b) -> a -> c
. ParticleState -> Vec
posVec forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. [a] -> a
last forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (a -> Bool) -> [a] -> [a]
takeWhile ParticleState -> Bool
aboveSeaLevel
          method :: NumericalMethod ParticleState DParticleState
method = forall s ds. Diff s ds => R -> (s -> ds) -> s -> s
rungeKutta4 R
0.01
          forcesNoAir :: [ParticleState -> Vec]
forcesNoAir
              = [ParticleState -> Vec
earthSurfaceGravity]
          forcesConstAir :: [ParticleState -> Vec]
forcesConstAir
              = [ParticleState -> Vec
earthSurfaceGravity, R -> R -> R -> ParticleState -> Vec
airResistance    R
drag R
airDensity R
area]
          forcesVarAir :: [ParticleState -> Vec]
forcesVarAir
              = [ParticleState -> Vec
earthSurfaceGravity, R -> R -> R -> ParticleState -> Vec
airResAtAltitude R
drag R
airDensity R
area]
          rangeNoAir :: R
rangeNoAir    = [ParticleState] -> R
range forall a b. (a -> b) -> a -> b
$ NumericalMethod ParticleState DParticleState
-> [ParticleState -> Vec] -> ParticleState -> [ParticleState]
statesPS NumericalMethod ParticleState DParticleState
method [ParticleState -> Vec]
forcesNoAir    forall {t}. t
stateInitial
          rangeConstAir :: R
rangeConstAir = [ParticleState] -> R
range forall a b. (a -> b) -> a -> b
$ NumericalMethod ParticleState DParticleState
-> [ParticleState -> Vec] -> ParticleState -> [ParticleState]
statesPS NumericalMethod ParticleState DParticleState
method [ParticleState -> Vec]
forcesConstAir forall {t}. t
stateInitial
          rangeVarAir :: R
rangeVarAir   = [ParticleState] -> R
range forall a b. (a -> b) -> a -> b
$ NumericalMethod ParticleState DParticleState
-> [ParticleState -> Vec] -> ParticleState -> [ParticleState]
statesPS NumericalMethod ParticleState DParticleState
method [ParticleState -> Vec]
forcesVarAir   forall {t}. t
stateInitial
      in forall a. HasCallStack => a
undefined R
rangeNoAir R
rangeConstAir R
rangeVarAir

halleyUpdate :: TimeStep
             -> ParticleState -> ParticleState
halleyUpdate :: R -> ParticleState -> ParticleState
halleyUpdate R
dt
    = NumericalMethod ParticleState DParticleState
-> [ParticleState -> Vec] -> ParticleState -> ParticleState
updatePS (R -> NumericalMethod ParticleState DParticleState
eulerCromerPS R
dt) [ParticleState -> Vec
sunGravity]

halleyInitial :: ParticleState
halleyInitial :: ParticleState
halleyInitial = ParticleState { mass :: R
mass     = R
2.2e14            -- kg
                              , charge :: R
charge   = R
0
                              , time :: R
time     = R
0
                              , posVec :: Vec
posVec   = R
8.766e10 R -> Vec -> Vec
*^ Vec
iHat  -- m
                              , velocity :: Vec
velocity = R
54569 R -> Vec -> Vec
*^ Vec
jHat }   -- m/s

baseballForces :: [OneBodyForce]
baseballForces :: [ParticleState -> Vec]
baseballForces
    = let area :: R
area = forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* (R
0.074 forall a. Fractional a => a -> a -> a
/ R
2) forall a. Floating a => a -> a -> a
** R
2
      in [ParticleState -> Vec
earthSurfaceGravity
         ,R -> R -> R -> ParticleState -> Vec
airResistance R
0.3 R
1.225 R
area]

baseballTrajectory :: R  -- time step
                   -> R  -- initial speed
                   -> R  -- launch angle in degrees
                   -> [(R,R)]  -- (y,z) pairs
baseballTrajectory :: R -> R -> R -> [(R, R)]
baseballTrajectory R
dt R
v0 R
thetaDeg
    = let thetaRad :: R
thetaRad = R
thetaDeg forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ R
180
          vy0 :: R
vy0 = R
v0 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos R
thetaRad
          vz0 :: R
vz0 = R
v0 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin R
thetaRad
          initialState :: ParticleState
initialState
              = ParticleState { mass :: R
mass     = R
0.145
                              , charge :: R
charge   = R
0
                              , time :: R
time     = R
0
                              , posVec :: Vec
posVec   = Vec
zeroV
                              , velocity :: Vec
velocity = R -> R -> R -> Vec
vec R
0 R
vy0 R
vz0 }
      in [ParticleState] -> [(R, R)]
trajectory forall a b. (a -> b) -> a -> b
$ [ParticleState] -> [ParticleState]
zGE0 forall a b. (a -> b) -> a -> b
$
         NumericalMethod ParticleState DParticleState
-> [ParticleState -> Vec] -> ParticleState -> [ParticleState]
statesPS (R -> NumericalMethod ParticleState DParticleState
eulerCromerPS R
dt) [ParticleState -> Vec]
baseballForces ParticleState
initialState

zGE0 :: [ParticleState] -> [ParticleState]
zGE0 :: [ParticleState] -> [ParticleState]
zGE0 = forall a. (a -> Bool) -> [a] -> [a]
takeWhile (\(ParticleState R
_ R
_ R
_ Vec
r Vec
_) -> Vec -> R
zComp Vec
r forall a. Ord a => a -> a -> Bool
>= R
0)

trajectory :: [ParticleState] -> [(R,R)]
trajectory :: [ParticleState] -> [(R, R)]
trajectory [ParticleState]
sts = [(Vec -> R
yComp Vec
r,Vec -> R
zComp Vec
r) | (ParticleState R
_ R
_ R
_ Vec
r Vec
_) <- [ParticleState]
sts]

baseballRange :: R  -- time step
              -> R  -- initial speed
              -> R  -- launch angle in degrees
              -> R  -- range
baseballRange :: R -> R -> R -> R
baseballRange R
dt R
v0 R
thetaDeg
    = let (R
y,R
_) = forall a. [a] -> a
last forall a b. (a -> b) -> a -> b
$ R -> R -> R -> [(R, R)]
baseballTrajectory R
dt R
v0 R
thetaDeg
      in R
y

bestAngle :: (R,R)
bestAngle :: (R, R)
bestAngle
    = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [(R -> R -> R -> R
baseballRange R
0.01 R
45 R
thetaDeg,R
thetaDeg) |
               R
thetaDeg <- [R
30,R
31..R
60]]

projectileUpdate :: TimeStep
                 -> ParticleState  -- old state
                 -> ParticleState  -- new state
projectileUpdate :: R -> ParticleState -> ParticleState
projectileUpdate R
dt
    = NumericalMethod ParticleState DParticleState
-> [ParticleState -> Vec] -> ParticleState -> ParticleState
updatePS (R -> NumericalMethod ParticleState DParticleState
eulerCromerPS R
dt) [ParticleState -> Vec]
baseballForces

projectileInitial :: [String] -> ParticleState
projectileInitial :: [String] -> ParticleState
projectileInitial []        = forall a. HasCallStack => String -> a
error String
"Please supply initial speed and angle."
projectileInitial [String
_]       = forall a. HasCallStack => String -> a
error String
"Please supply initial speed and angle."
projectileInitial (String
_:String
_:String
_:[String]
_)
    = forall a. HasCallStack => String -> a
error String
"First argument is speed.  Second is angle in degrees."
projectileInitial (String
arg1:String
arg2:[String]
_)
    = let v0 :: R
v0       = forall a. Read a => String -> a
read String
arg1 :: R       -- initial speed, m/s
          angleDeg :: R
angleDeg = forall a. Read a => String -> a
read String
arg2 :: R       -- initial angle, degrees
          theta :: R
theta    = R
angleDeg forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ R
180  -- in radians
      in ParticleState
defaultParticleState
             { mass :: R
mass     = R
0.145  -- kg
             , posVec :: Vec
posVec   = Vec
zeroV
             , velocity :: Vec
velocity = R -> R -> R -> Vec
vec R
0 (R
v0 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos R
theta) (R
v0 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin R
theta)
             }

protonUpdate :: TimeStep -> ParticleState -> ParticleState
protonUpdate :: R -> ParticleState -> ParticleState
protonUpdate R
dt
    = NumericalMethod ParticleState DParticleState
-> [ParticleState -> Vec] -> ParticleState -> ParticleState
updatePS (forall s ds. Diff s ds => R -> (s -> ds) -> s -> s
rungeKutta4 R
dt) [Vec -> Vec -> ParticleState -> Vec
uniformLorentzForce Vec
zeroV (R
3e-8 R -> Vec -> Vec
*^ Vec
kHat)]

protonInitial :: ParticleState
protonInitial :: ParticleState
protonInitial
    = ParticleState
defaultParticleState { mass :: R
mass     = R
1.672621898e-27  -- kg
                           , charge :: R
charge   = R
1.602176621e-19  -- C
                           , posVec :: Vec
posVec   = Vec
zeroV
                           , velocity :: Vec
velocity = R
1.5R -> Vec -> Vec
*^Vec
jHat Vec -> Vec -> Vec
^+^ R
0.3R -> Vec -> Vec
*^Vec
kHat  -- m/s
                           }

apR :: R
apR :: R
apR = R
0.04  -- meters

wallForce :: OneBodyForce
wallForce :: ParticleState -> Vec
wallForce ParticleState
ps
    = let m :: R
m = ParticleState -> R
mass ParticleState
ps
          r :: Vec
r = ParticleState -> Vec
posVec ParticleState
ps
          x :: R
x = Vec -> R
xComp Vec
r
          y :: R
y = Vec -> R
yComp Vec
r
          z :: R
z = Vec -> R
zComp Vec
r
          v :: Vec
v = ParticleState -> Vec
velocity ParticleState
ps
          timeStep :: R
timeStep = R
5e-4 forall a. Fractional a => a -> a -> a
/ R
60
      in if R
y forall a. Ord a => a -> a -> Bool
>= R
1 Bool -> Bool -> Bool
&& R
y forall a. Ord a => a -> a -> Bool
< R
1.1 Bool -> Bool -> Bool
&& forall a. Floating a => a -> a
sqrt (R
xforall a. Floating a => a -> a -> a
**R
2 forall a. Num a => a -> a -> a
+ R
zforall a. Floating a => a -> a -> a
**R
2) forall a. Ord a => a -> a -> Bool
> R
apR
         then (-R
m) R -> Vec -> Vec
*^ (Vec
v Vec -> R -> Vec
^/ R
timeStep)
         else Vec
zeroV

energy :: ParticleState -> R
energy :: ParticleState -> R
energy ParticleState
ps = forall a. HasCallStack => a
undefined ParticleState
ps

firstOrbit :: ParticleState -> Bool
firstOrbit :: ParticleState -> Bool
firstOrbit ParticleState
st
    = let year :: R
year = R
365.25 forall a. Num a => a -> a -> a
* R
24 forall a. Num a => a -> a -> a
* R
60 forall a. Num a => a -> a -> a
* R
60
      in ParticleState -> R
time ParticleState
st forall a. Ord a => a -> a -> Bool
< R
50 forall a. Num a => a -> a -> a
* R
year Bool -> Bool -> Bool
|| Vec -> R
yComp (ParticleState -> Vec
posVec ParticleState
st) forall a. Ord a => a -> a -> Bool
<= R
0

-- | Given a list of forces, return a differential equation
--   based on the theory of special relativity.
relativityPS :: [OneBodyForce]
             -> ParticleState -> DParticleState  -- a differential equation
relativityPS :: [ParticleState -> Vec] -> ParticleState -> DParticleState
relativityPS [ParticleState -> Vec]
fs ParticleState
st
    = let fNet :: Vec
fNet = [Vec] -> Vec
sumV [ParticleState -> Vec
f ParticleState
st | ParticleState -> Vec
f <- [ParticleState -> Vec]
fs]
          c :: R
c = R
299792458  -- m / s
          m :: R
m = ParticleState -> R
mass ParticleState
st
          v :: Vec
v = ParticleState -> Vec
velocity ParticleState
st
          u :: Vec
u = Vec
v Vec -> R -> Vec
^/ R
c
          acc :: Vec
acc = forall a. Floating a => a -> a
sqrt (R
1 forall a. Num a => a -> a -> a
- Vec
u Vec -> Vec -> R
<.> Vec
u) R -> Vec -> Vec
*^ (Vec
fNet Vec -> Vec -> Vec
^-^ (Vec
fNet Vec -> Vec -> R
<.> Vec
u) R -> Vec -> Vec
*^ Vec
u) Vec -> R -> Vec
^/ R
m
      in DParticleState { dmdt :: R
dmdt = R
0    -- dm/dt
                        , dqdt :: R
dqdt = R
0    -- dq/dt
                        , dtdt :: R
dtdt = R
1    -- dt/dt
                        , drdt :: Vec
drdt = Vec
v    -- dr/dt
                        , dvdt :: Vec
dvdt = Vec
acc  -- dv/vt
                        }

twoProtUpdate :: TimeStep
              -> (ParticleState,ParticleState)
              -> (ParticleState,ParticleState)
twoProtUpdate :: R
-> (ParticleState, ParticleState) -> (ParticleState, ParticleState)
twoProtUpdate R
dt (ParticleState
stN,ParticleState
stR)
    = let forces :: [ParticleState -> Vec]
forces = [Vec -> Vec -> ParticleState -> Vec
uniformLorentzForce Vec
zeroV Vec
kHat]
      in (forall s ds. Diff s ds => R -> (s -> ds) -> s -> s
rungeKutta4 R
dt ([ParticleState -> Vec] -> ParticleState -> DParticleState
newtonSecondPS [ParticleState -> Vec]
forces) ParticleState
stN
         ,forall s ds. Diff s ds => R -> (s -> ds) -> s -> s
rungeKutta4 R
dt ([ParticleState -> Vec] -> ParticleState -> DParticleState
relativityPS   [ParticleState -> Vec]
forces) ParticleState
stR)

twoProtInitial :: (ParticleState,ParticleState)
twoProtInitial :: (ParticleState, ParticleState)
twoProtInitial
    = let c :: R
c = R
299792458  -- m/s
          pInit :: ParticleState
pInit = ParticleState
protonInitial { velocity :: Vec
velocity = R
0.8 R -> Vec -> Vec
*^ R
c R -> Vec -> Vec
*^ Vec
jHat }
      in (ParticleState
pInit,ParticleState
pInit)

relativityPS' :: R  -- c
              -> [OneBodyForce]
              -> ParticleState -> DParticleState
relativityPS' :: R -> [ParticleState -> Vec] -> ParticleState -> DParticleState
relativityPS' R
c [ParticleState -> Vec]
fs ParticleState
st = forall a. HasCallStack => a
undefined R
c [ParticleState -> Vec]
fs ParticleState
st