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

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

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

module LPFPCore.Lorentz where

import LPFPCore.SimpleVec ( R, Vec, (^+^), (*^), (^*), (^/), (><), zeroV )
import LPFPCore.Mechanics1D ( RealVectorSpace(..), Diff(..), rungeKutta4 )
import LPFPCore.Mechanics3D ( HasTime(..) )
import LPFPCore.CoordinateSystems ( Position(..), VectorField, origin
                         , shiftPosition, addVectorFields )

data ParticleFieldState = ParticleFieldState { ParticleFieldState -> R
mass          :: R
                                             , ParticleFieldState -> R
charge        :: R
                                             , ParticleFieldState -> R
time          :: R
                                             , ParticleFieldState -> Position
position      :: Position
                                             , ParticleFieldState -> Vec
velocity      :: Vec
                                             , ParticleFieldState -> VectorField
electricField :: VectorField
                                             , ParticleFieldState -> VectorField
magneticField :: VectorField }

data DParticleFieldState = DParticleFieldState { DParticleFieldState -> R
dmdt :: R
                                               , DParticleFieldState -> R
dqdt :: R
                                               , DParticleFieldState -> R
dtdt :: R
                                               , DParticleFieldState -> Vec
drdt :: Vec
                                               , DParticleFieldState -> Vec
dvdt :: Vec
                                               , DParticleFieldState -> VectorField
dEdt :: VectorField
                                               , DParticleFieldState -> VectorField
dBdt :: VectorField }

instance RealVectorSpace DParticleFieldState where
    DParticleFieldState
dst1 +++ :: DParticleFieldState -> DParticleFieldState -> DParticleFieldState
+++ DParticleFieldState
dst2
        = DParticleFieldState { dmdt :: R
dmdt = DParticleFieldState -> R
dmdt DParticleFieldState
dst1  forall a. Num a => a -> a -> a
+  DParticleFieldState -> R
dmdt DParticleFieldState
dst2
                              , dqdt :: R
dqdt = DParticleFieldState -> R
dqdt DParticleFieldState
dst1  forall a. Num a => a -> a -> a
+  DParticleFieldState -> R
dqdt DParticleFieldState
dst2
                              , dtdt :: R
dtdt = DParticleFieldState -> R
dtdt DParticleFieldState
dst1  forall a. Num a => a -> a -> a
+  DParticleFieldState -> R
dtdt DParticleFieldState
dst2
                              , drdt :: Vec
drdt = DParticleFieldState -> Vec
drdt DParticleFieldState
dst1 Vec -> Vec -> Vec
^+^ DParticleFieldState -> Vec
drdt DParticleFieldState
dst2
                              , dvdt :: Vec
dvdt = DParticleFieldState -> Vec
dvdt DParticleFieldState
dst1 Vec -> Vec -> Vec
^+^ DParticleFieldState -> Vec
dvdt DParticleFieldState
dst2
                              , dEdt :: VectorField
dEdt = [VectorField] -> VectorField
addVectorFields [DParticleFieldState -> VectorField
dEdt DParticleFieldState
dst1, DParticleFieldState -> VectorField
dEdt DParticleFieldState
dst2]
                              , dBdt :: VectorField
dBdt = [VectorField] -> VectorField
addVectorFields [DParticleFieldState -> VectorField
dBdt DParticleFieldState
dst1, DParticleFieldState -> VectorField
dBdt DParticleFieldState
dst2]
                              }
    scale :: R -> DParticleFieldState -> DParticleFieldState
scale R
w DParticleFieldState
dst
        = DParticleFieldState { dmdt :: R
dmdt = R
w forall a. Num a => a -> a -> a
*  DParticleFieldState -> R
dmdt DParticleFieldState
dst
                              , dqdt :: R
dqdt = R
w forall a. Num a => a -> a -> a
*  DParticleFieldState -> R
dqdt DParticleFieldState
dst
                              , dtdt :: R
dtdt = R
w forall a. Num a => a -> a -> a
*  DParticleFieldState -> R
dtdt DParticleFieldState
dst
                              , drdt :: Vec
drdt = R
w R -> Vec -> Vec
*^ DParticleFieldState -> Vec
drdt DParticleFieldState
dst
                              , dvdt :: Vec
dvdt = R
w R -> Vec -> Vec
*^ DParticleFieldState -> Vec
dvdt DParticleFieldState
dst
                              , dEdt :: VectorField
dEdt = (R
w R -> Vec -> Vec
*^) forall b c a. (b -> c) -> (a -> b) -> a -> c
. (DParticleFieldState -> VectorField
dEdt DParticleFieldState
dst)
                              , dBdt :: VectorField
dBdt = (R
w R -> Vec -> Vec
*^) forall b c a. (b -> c) -> (a -> b) -> a -> c
. (DParticleFieldState -> VectorField
dBdt DParticleFieldState
dst)
                              }

instance Diff ParticleFieldState DParticleFieldState where
    shift :: R
-> DParticleFieldState -> ParticleFieldState -> ParticleFieldState
shift R
dt DParticleFieldState
dst ParticleFieldState
st
        = ParticleFieldState
          { mass :: R
mass          = ParticleFieldState -> R
mass     ParticleFieldState
st  forall a. Num a => a -> a -> a
+  DParticleFieldState -> R
dmdt DParticleFieldState
dst  forall a. Num a => a -> a -> a
* R
dt
          , charge :: R
charge        = ParticleFieldState -> R
charge   ParticleFieldState
st  forall a. Num a => a -> a -> a
+  DParticleFieldState -> R
dqdt DParticleFieldState
dst  forall a. Num a => a -> a -> a
* R
dt
          , time :: R
time          = ParticleFieldState -> R
time     ParticleFieldState
st  forall a. Num a => a -> a -> a
+  DParticleFieldState -> R
dtdt DParticleFieldState
dst  forall a. Num a => a -> a -> a
* R
dt
          , position :: Position
position      = Vec -> Position -> Position
shiftPosition (DParticleFieldState -> Vec
drdt DParticleFieldState
dst Vec -> R -> Vec
^* R
dt) (ParticleFieldState -> Position
position ParticleFieldState
st)
          , velocity :: Vec
velocity      = ParticleFieldState -> Vec
velocity ParticleFieldState
st Vec -> Vec -> Vec
^+^ DParticleFieldState -> Vec
dvdt DParticleFieldState
dst Vec -> R -> Vec
^* R
dt
          , electricField :: VectorField
electricField = \Position
r -> ParticleFieldState -> VectorField
electricField ParticleFieldState
st Position
r Vec -> Vec -> Vec
^+^ DParticleFieldState -> VectorField
dEdt DParticleFieldState
dst Position
r Vec -> R -> Vec
^* R
dt
          , magneticField :: VectorField
magneticField = \Position
r -> ParticleFieldState -> VectorField
magneticField ParticleFieldState
st Position
r Vec -> Vec -> Vec
^+^ DParticleFieldState -> VectorField
dBdt DParticleFieldState
dst Position
r Vec -> R -> Vec
^* R
dt
          }

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

lorentzForce :: ParticleFieldState -> Vec
lorentzForce :: ParticleFieldState -> Vec
lorentzForce (ParticleFieldState R
_m R
q R
_t Position
r Vec
v VectorField
eF VectorField
bF)
    = R
q R -> Vec -> Vec
*^ (VectorField
eF Position
r Vec -> Vec -> Vec
^+^ Vec
v Vec -> Vec -> Vec
>< VectorField
bF Position
r)

newtonSecondPFS :: ParticleFieldState -> DParticleFieldState
newtonSecondPFS :: ParticleFieldState -> DParticleFieldState
newtonSecondPFS ParticleFieldState
st
    = let v :: Vec
v = ParticleFieldState -> Vec
velocity ParticleFieldState
st
          a :: Vec
a = ParticleFieldState -> Vec
lorentzForce ParticleFieldState
st Vec -> R -> Vec
^/ ParticleFieldState -> R
mass ParticleFieldState
st
      in DParticleFieldState { 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
a            -- dv/dt
                             , dEdt :: VectorField
dEdt = forall a b. a -> b -> a
const Vec
zeroV  -- dE/dt
                             , dBdt :: VectorField
dBdt = forall a b. a -> b -> a
const Vec
zeroV  -- dB/dt
                             }

pfsUpdate :: R  -- time step
          -> ParticleFieldState -> ParticleFieldState
pfsUpdate :: R -> ParticleFieldState -> ParticleFieldState
pfsUpdate R
dt = forall s ds. Diff s ds => R -> (s -> ds) -> s -> s
rungeKutta4 R
dt ParticleFieldState -> DParticleFieldState
newtonSecondPFS

defaultPFS :: ParticleFieldState
defaultPFS :: ParticleFieldState
defaultPFS = ParticleFieldState { mass :: R
mass          = R
0
                                , charge :: R
charge        = R
0
                                , time :: R
time          = R
0
                                , position :: Position
position      = Position
origin
                                , velocity :: Vec
velocity      = Vec
zeroV
                                , electricField :: VectorField
electricField = forall a b. a -> b -> a
const Vec
zeroV
                                , magneticField :: VectorField
magneticField = forall a b. a -> b -> a
const Vec
zeroV }

scalePos :: R -> Position -> Position
scalePos :: R -> Position -> Position
scalePos R
metersPerVis (Cart R
x R
y R
z)
    = R -> R -> R -> Position
Cart (R
xforall a. Fractional a => a -> a -> a
/R
metersPerVis) (R
yforall a. Fractional a => a -> a -> a
/R
metersPerVis) (R
zforall a. Fractional a => a -> a -> a
/R
metersPerVis)

newtonSecondPFS' :: [ParticleFieldState -> Vec]
                 -> ParticleFieldState -> DParticleFieldState
newtonSecondPFS' :: [ParticleFieldState -> Vec]
-> ParticleFieldState -> DParticleFieldState
newtonSecondPFS' [ParticleFieldState -> Vec]
fs ParticleFieldState
st = forall a. HasCallStack => a
undefined [ParticleFieldState -> Vec]
fs ParticleFieldState
st