{-# OPTIONS -Wall #-}

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

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

module LPFPCore.Electricity where

import LPFPCore.SimpleVec
    ( Vec(..), R, (*^), iHat )
import LPFPCore.Mechanics3D
    ( ParticleState(..), defaultParticleState )
import LPFPCore.MultipleObjects
    ( TwoBodyForce, MultiParticleState(..), Force(..), statesMPS
    , eulerCromerMPS, centralForce )

type Charge = R

elementaryCharge :: Charge
elementaryCharge :: R
elementaryCharge = R
1.602176634e-19  -- in Coulombs

coulombMagnitude :: Charge -> Charge -> R -> R
coulombMagnitude :: R -> R -> R -> R
coulombMagnitude R
q1 R
q2 R
r
    = let k :: R
k = R
9e9  -- in N m^2 / C^2
      in R
k forall a. Num a => a -> a -> a
* forall a. Num a => a -> a
abs (R
q1 forall a. Num a => a -> a -> a
* R
q2) forall a. Fractional a => a -> a -> a
/ R
rforall a. Floating a => a -> a -> a
**R
2

coulombForce :: TwoBodyForce
coulombForce :: TwoBodyForce
coulombForce ParticleState
st1 ParticleState
st2
    = let k :: R
k = R
9e9  -- N m^2 / C^2
          q1 :: R
q1 = ParticleState -> R
charge ParticleState
st1
          q2 :: R
q2 = ParticleState -> R
charge ParticleState
st2
      in (R -> R) -> TwoBodyForce
centralForce (\R
r -> R
k forall a. Num a => a -> a -> a
* R
q1 forall a. Num a => a -> a -> a
* R
q2 forall a. Fractional a => a -> a -> a
/ R
rforall a. Floating a => a -> a -> a
**R
2) ParticleState
st1 ParticleState
st2

twoProtonStates :: R                     -- time step
                -> MultiParticleState    -- initial 2-particle state
                -> [MultiParticleState]  -- infinite list of states
twoProtonStates :: R -> MultiParticleState -> [MultiParticleState]
twoProtonStates R
dt
    = NumericalMethod MultiParticleState DMultiParticleState
-> [Force] -> MultiParticleState -> [MultiParticleState]
statesMPS (R -> NumericalMethod MultiParticleState DMultiParticleState
eulerCromerMPS R
dt) [Int -> Int -> TwoBodyForce -> Force
InternalForce Int
1 Int
0 TwoBodyForce
coulombForce]

-- protons are released from rest
initialTwoProtonState :: R  -- initial separation
                      -> MultiParticleState
initialTwoProtonState :: R -> MultiParticleState
initialTwoProtonState R
d
    = let protonMass :: R
protonMass = R
1.673e-27  -- in kg
      in [ParticleState] -> MultiParticleState
MPS [ParticleState
defaultParticleState { mass :: R
mass   = R
protonMass
                                   , charge :: R
charge = R
elementaryCharge
                                   , posVec :: Vec
posVec = (-R
dforall a. Fractional a => a -> a -> a
/R
2) R -> Vec -> Vec
*^ Vec
iHat
                                   }
             ,ParticleState
defaultParticleState { mass :: R
mass   = R
protonMass
                                   , charge :: R
charge = R
elementaryCharge
                                   , posVec :: Vec
posVec = ( R
dforall a. Fractional a => a -> a -> a
/R
2) R -> Vec -> Vec
*^ Vec
iHat
                                   }
             ]

oneProtonVelocity :: R        -- dt
                  -> R        -- starting separation
                  -> [(R,R)]  -- (time,velocity) pairs
oneProtonVelocity :: R -> R -> [(R, R)]
oneProtonVelocity R
dt R
d
    = let state0 :: MultiParticleState
state0 = R -> MultiParticleState
initialTwoProtonState R
d
      in [(ParticleState -> R
time ParticleState
st2, Vec -> R
xComp forall a b. (a -> b) -> a -> b
$ ParticleState -> Vec
velocity ParticleState
st2)
              | MPS [ParticleState
_,ParticleState
st2] <- R -> MultiParticleState -> [MultiParticleState]
twoProtonStates R
dt MultiParticleState
state0]

tvPairs :: [(R,R)]
tvPairs :: [(R, R)]
tvPairs = forall a. (a -> Bool) -> [a] -> [a]
takeWhile (\(R
t,R
_) -> R
t forall a. Ord a => a -> a -> Bool
<= R
2e-2) forall a b. (a -> b) -> a -> b
$
          R -> R -> [(R, R)]
oneProtonVelocity R
1e-5 R
1e-2

oneProtonPosition :: R        -- dt
                  -> R        -- starting separation
                  -> [(R,R)]  -- (time,position) pairs
oneProtonPosition :: R -> R -> [(R, R)]
oneProtonPosition R
dt R
d
    = forall a. HasCallStack => a
undefined R
dt R
d