{-# OPTIONS -Wall #-} {- | Module : LPFP.MOExamples Copyright : (c) Scott N. Walck 2023 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck Stability : stable Code from chapter 20 of the book Learn Physics with Functional Programming -} module LPFP.MOExamples where import LPFP.SimpleVec ( R, Vec, (^+^), (^-^), (*^), vec, zeroV, magnitude , sumV, iHat, jHat, kHat, xComp, yComp, zComp ) import LPFP.Mechanics1D ( TimeStep, NumericalMethod, euler, rungeKutta4 ) import LPFP.Mechanics3D ( ParticleState(..), HasTime(..), defaultParticleState , earthSurfaceGravity, customLabel, orient, disk ) import LPFP.MultipleObjects ( MultiParticleState(..), DMultiParticleState, Force(..), TwoBodyForce , newtonSecondMPS, updateMPS, statesMPS, eulerCromerMPS , linearSpring, fixedLinearSpring, billiardForce ) import Graphics.Gnuplot.Simple import qualified Graphics.Gloss as G import qualified Vis as V twoSpringsForces :: [Force] twoSpringsForces = [ExternalForce 0 (fixedLinearSpring 100 0.5 zeroV) ,InternalForce 0 1 (linearSpring 100 0.5) ,ExternalForce 0 earthSurfaceGravity ,ExternalForce 1 earthSurfaceGravity ] twoSpringsInitial :: MultiParticleState twoSpringsInitial = MPS [defaultParticleState { mass = 2 , posVec = 0.4 *^ jHat ^-^ 0.3 *^ kHat } ,defaultParticleState { mass = 3 , posVec = 0.4 *^ jHat ^-^ 0.8 *^ kHat } ] twoSpringsUpdate :: TimeStep -> MultiParticleState -- old state -> MultiParticleState -- new state twoSpringsUpdate dt = updateMPS (eulerCromerMPS dt) twoSpringsForces kineticEnergy :: ParticleState -> R kineticEnergy st = let m = mass st v = magnitude (velocity st) in (1/2) * m * v**2 systemKE :: MultiParticleState -> R systemKE (MPS sts) = sum [kineticEnergy st | st <- 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 k re st1 st2 = let r1 = posVec st1 r2 = posVec st2 r21 = r2 ^-^ r1 r21mag = magnitude r21 in k * (r21mag - re)**2 / 2 -- z direction is toward the sky -- assumes SI units earthSurfaceGravityPE :: ParticleState -> R earthSurfaceGravityPE st = let g = 9.80665 -- m/s^2 m = mass st z = zComp (posVec st) in m * g * z twoSpringsPE :: MultiParticleState -> R twoSpringsPE (MPS sts) = linearSpringPE 100 0.5 defaultParticleState (sts !! 0) + linearSpringPE 100 0.5 (sts !! 0) (sts !! 1) + earthSurfaceGravityPE (sts !! 0) + earthSurfaceGravityPE (sts !! 1) twoSpringsME :: MultiParticleState -> R twoSpringsME mpst = systemKE mpst + twoSpringsPE mpst billiardForces :: R -> [Force] billiardForces k = [InternalForce 0 1 (billiardForce k (2*ballRadius))] ballRadius :: R ballRadius = 0.03 -- 6cm diameter = 0.03m radius billiardDiffEq :: R -> MultiParticleState -> DMultiParticleState billiardDiffEq k = newtonSecondMPS $ billiardForces k billiardUpdate :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState) -> R -- k -> TimeStep -- dt -> MultiParticleState -> MultiParticleState billiardUpdate nMethod k dt = updateMPS (nMethod dt) (billiardForces k) billiardEvolver :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState) -> R -- k -> TimeStep -- dt -> MultiParticleState -> [MultiParticleState] billiardEvolver nMethod k dt = statesMPS (nMethod dt) (billiardForces k) billiardInitial :: MultiParticleState billiardInitial = let ballMass = 0.160 -- 160g in MPS [defaultParticleState { mass = ballMass , posVec = zeroV , velocity = 0.2 *^ iHat } ,defaultParticleState { mass = ballMass , posVec = iHat ^+^ 0.02 *^ jHat , velocity = zeroV } ] billiardStates :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState) -> R -- k -> TimeStep -- dt -> [MultiParticleState] billiardStates nMethod k dt = statesMPS (nMethod dt) (billiardForces k) billiardInitial billiardStatesFinite :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState) -> R -- k -> TimeStep -- dt -> [MultiParticleState] billiardStatesFinite nMethod k dt = takeWhile (\st -> timeOf st <= 10) (billiardStates nMethod k dt) momentum :: ParticleState -> Vec momentum st = let m = mass st v = velocity st in m *^ v systemP :: MultiParticleState -> Vec systemP (MPS sts) = sumV [momentum st | st <- sts] percentChangePMag :: [MultiParticleState] -> R percentChangePMag mpsts = let p0 = systemP (head mpsts) p1 = systemP (last mpsts) in 100 * magnitude (p1 ^-^ p0) / magnitude p0 sigFigs :: Int -> R -> Float sigFigs n x = let expon :: Int expon = floor (logBase 10 x) - n + 1 toInt :: R -> Int toInt = round in (10^^expon *) $ fromIntegral $ toInt (10^^(-expon) * x) data Justification = LJ | RJ deriving Show data Table a = Table Justification [[a]] instance Show a => Show (Table a) where show (Table j xss) = let pairWithLength x = let str = show x in (str, length str) pairss = map (map pairWithLength) xss maxLength = maximum (map maximum (map (map snd) pairss)) showPair (str,len) = case j of LJ -> str ++ replicate (maxLength + 1 - len) ' ' RJ -> replicate (maxLength + 1 - len) ' ' ++ str showLine pairs = concatMap showPair pairs ++ "\n" in init $ concatMap showLine pairss pTable :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState) -> [R] -- ks -> [TimeStep] -- dts -> Table Float pTable nMethod ks dts = Table LJ [[sigFigs 2 $ percentChangePMag (billiardStatesFinite nMethod k dt) | dt <- dts] | k <- ks] pTableEu :: [R] -- ks -> [TimeStep] -- dts -> Table Float pTableEu = pTable euler systemKEWithTime :: IO () systemKEWithTime = let timeKEPairsEC = [(timeOf mpst, systemKE mpst) | mpst <- billiardStatesFinite eulerCromerMPS 30 0.03] timeKEPairsRK4 = [(timeOf mpst, systemKE mpst) | mpst <- billiardStatesFinite rungeKutta4 30 0.03] in plotPaths [Key Nothing ,Title "System Kinetic Energy versus Time" ,XLabel "Time (s)" ,YLabel "System Kinetic Energy (J)" ,XRange (4,6) ,PNG "SystemKE.png" ,customLabel (4.1,0.0026) "dt = 0.03 s" ,customLabel (4.1,0.0025) "k = 30 N/m" ,customLabel (5.4,0.00329) "Euler-Cromer" ,customLabel (5.4,0.00309) "Runge-Kutta 4" ] [timeKEPairsEC,timeKEPairsRK4] percentChangeKE :: [MultiParticleState] -> R percentChangeKE mpsts = let ke0 = systemKE (head mpsts) ke1 = systemKE (last mpsts) in 100 * (ke1 - ke0) / ke0 tenths :: R -> Float tenths = let toInt :: R -> Int toInt = round in (/ 10) . fromIntegral . toInt . (* 10) keTable :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState) -> [R] -- ks -> [TimeStep] -- dts -> Table Float keTable nMethod ks dts = Table RJ [[tenths $ percentChangeKE (billiardStatesFinite nMethod k dt) | dt <- dts] | k <- ks] contactSteps :: [MultiParticleState] -> Int contactSteps = length . takeWhile inContact . dropWhile (not . inContact) inContact :: MultiParticleState -> Bool inContact (MPS sts) = let r = magnitude $ posVec (sts !! 0) ^-^ posVec (sts !! 1) in r < 2 * ballRadius contactTable :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState) -> [R] -- ks -> [TimeStep] -- dts -> Table Int contactTable nMethod ks dts = Table RJ [[contactSteps (billiardStatesFinite nMethod k dt) | dt <- dts] | k <- ks] closest :: [MultiParticleState] -> R closest = minimum . map separation separation :: MultiParticleState -> R separation (MPS sts) = magnitude $ posVec (sts !! 0) ^-^ posVec (sts !! 1) closestTable :: (TimeStep -> NumericalMethod MultiParticleState DMultiParticleState) -> [R] -- ks -> [TimeStep] -- dts -> Table Float closestTable nMethod ks dts = Table RJ [[tenths $ (100*) $ closest (billiardStatesFinite nMethod k dt) | dt <- dts] | k <- ks] billiardPicture :: MultiParticleState -> G.Picture billiardPicture (MPS sts) = G.scale ppm ppm $ G.pictures [place st | st <- sts] where ppm = 300 -- pixels per meter place st = G.translate (xSt st) (ySt st) blueBall xSt = realToFrac . xComp . posVec ySt = realToFrac . yComp . posVec blueBall = G.Color G.blue (disk $ realToFrac ballRadius) -- 64 masses (0 to 63) -- There are 63 internal springs, 2 external springs forcesString :: [Force] forcesString = [ExternalForce 0 (fixedLinearSpring 5384 0 (vec 0 0 0)) ,ExternalForce 63 (fixedLinearSpring 5384 0 (vec 0.65 0 0))] ++ [InternalForce n (n+1) (linearSpring 5384 0) | n <- [0..62]] stringUpdate :: TimeStep -> MultiParticleState -- old state -> MultiParticleState -- new state stringUpdate dt = updateMPS (rungeKutta4 dt) forcesString stringInitialOvertone :: Int -> MultiParticleState stringInitialOvertone n = MPS [defaultParticleState { mass = 0.8293e-3 * 0.65 / 64 , posVec = x *^ iHat ^+^ y *^ jHat , velocity = zeroV } | x <- [0.01, 0.02 .. 0.64], let y = 0.005 * sin (fromIntegral n * pi * x / 0.65)] stringInitialPluck :: MultiParticleState stringInitialPluck = MPS [defaultParticleState { mass = 0.8293e-3 * 0.65 / 64 , posVec = x *^ iHat ^+^ y *^ jHat , velocity = zeroV } | x <- [0.01, 0.02 .. 0.64], let y = pluckEq x] where pluckEq :: R -> R pluckEq x | x <= 0.51 = 0.005 / (0.51 - 0.00) * (x - 0.00) | otherwise = 0.005 / (0.51 - 0.65) * (x - 0.65) mpsPos :: MultiParticleState -> IO () mpsPos = undefined mpsVel :: MultiParticleState -> IO () mpsVel = undefined dissipation :: R -- damping constant -> R -- threshold center separation -> TwoBodyForce dissipation b re st1 st2 = let r1 = posVec st1 r2 = posVec st2 v1 = velocity st1 v2 = velocity st2 r21 = r2 ^-^ r1 v21 = v2 ^-^ v1 in if magnitude r21 >= re then zeroV else (-b) *^ v21 animateGloss :: HasTime s => R -- time-scale factor -> (s -> G.Picture) -> [s] -> IO () animateGloss tsFactor displayFunc mpsts = let dtp = timeOf (mpsts !! 1) - timeOf (mpsts !! 0) n tp = round (tp / dtp) picFromAnimTime :: Float -> G.Picture picFromAnimTime ta = displayFunc (mpsts !! n (tsFactor * realToFrac ta)) displayMode = G.InWindow "My Window" (1000, 700) (10, 10) in G.animate displayMode G.black picFromAnimTime animateVis :: HasTime s => R -- time-scale factor -> (s -> V.VisObject R) -> [s] -> IO () animateVis tsFactor displayFunc mpsts = let dtp = timeOf (mpsts !! 1) - timeOf (mpsts !! 0) n tp = round (tp / dtp) picFromAnimTime :: Float -> V.VisObject R picFromAnimTime ta = displayFunc (mpsts !! n (tsFactor * realToFrac ta)) in V.animate V.defaultOpts (orient . picFromAnimTime)