diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,29 @@ +Copyright (c) 2022 Scott N. Walck . +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Scott N. Walck nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/LPFP.cabal b/LPFP.cabal new file mode 100644 --- /dev/null +++ b/LPFP.cabal @@ -0,0 +1,307 @@ +cabal-version: 1.12 + +name: LPFP +version: 1.0 +synopsis: Code for the book Learn Physics with Functional Programming +description: Haskell code to help the user learn mechanics of one particle, + mechanics of multiple interacting particles, and electromagnetic theory. +homepage: https://lpfp.io +author: Scott N. Walck +maintainer: walck@lvc.edu +copyright: 2023 Scott N. Walck +license: BSD3 +license-file: LICENSE +build-type: Simple +category: Physics + +library + exposed-modules: + LPFP + LPFP.SimpleVec + LPFP.Newton2 + LPFP.Mechanics1D + LPFP.Mechanics3D + LPFP.MultipleObjects + LPFP.MOExamples + LPFP.Electricity + LPFP.CoordinateSystems + LPFP.Geometry + LPFP.Integrals + LPFP.Charge + LPFP.ElectricField + LPFP.Current + LPFP.MagneticField + LPFP.Lorentz + LPFP.Maxwell + hs-source-dirs: src + build-depends: + base >=4.7 && <5, gnuplot, spatial-math, gloss, not-gloss, diagrams-lib + , diagrams-cairo, containers + default-language: Haskell2010 + +executable LPFP-hello + main-is: hello.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + default-language: Haskell2010 + +executable LPFP-MakeTrajectoryGraph + main-is: MakeTrajectoryGraph.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , gnuplot + default-language: Haskell2010 + +executable LPFP-GlossDisplay + main-is: GlossDisplay.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , gloss + default-language: Haskell2010 + +executable LPFP-GlossDisplay2 + main-is: GlossDisplay2.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , gloss + default-language: Haskell2010 + +executable LPFP-GlossAnimate + main-is: GlossAnimate.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , gloss + default-language: Haskell2010 + +executable LPFP-GlossSimulate + main-is: GlossSimulate.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , gloss + default-language: Haskell2010 + +executable LPFP-GlossSimulate2 + main-is: GlossSimulate2.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , gloss + default-language: Haskell2010 + +executable LPFP-VisDisplay + main-is: VisDisplay.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , not-gloss + default-language: Haskell2010 + +executable LPFP-VisDisplay2 + main-is: VisDisplay2.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , not-gloss + , linear + default-language: Haskell2010 + +executable LPFP-VisDisplay3 + main-is: VisDisplay3.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , not-gloss + , linear + , spatial-math + default-language: Haskell2010 + +executable LPFP-VisAnimate + main-is: VisAnimate.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , not-gloss + , spatial-math + default-language: Haskell2010 + +executable LPFP-VisSimulate + main-is: VisSimulate.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , not-gloss + default-language: Haskell2010 + +executable LPFP-GlossHalley + main-is: GlossHalley.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , gloss + default-language: Haskell2010 + +executable LPFP-GlossProjectile + main-is: GlossProjectile.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , gloss + default-language: Haskell2010 + +executable LPFP-VisProjectile + main-is: VisProjectile.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , not-gloss + default-language: Haskell2010 + +executable LPFP-Mech3Proton + main-is: Mech3Proton.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , gloss + default-language: Haskell2010 + +executable LPFP-GlossProton + main-is: GlossProton.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , gloss + default-language: Haskell2010 + +executable LPFP-VisTwoSprings + main-is: VisTwoSprings.hs + hs-source-dirs: app + build-depends: LPFP, base >=4.7 && <5, not-gloss + default-language: Haskell2010 + +executable LPFP-GlossBilliard + main-is: GlossBilliard.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , gloss + default-language: Haskell2010 + +executable LPFP-GlossWave + main-is: GlossWave.hs + hs-source-dirs: app + build-depends: LPFP, base >=4.7 && <5, gloss + default-language: Haskell2010 + +executable LPFP-GnuplotWave + main-is: GnuplotWave.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , gnuplot + default-language: Haskell2010 + +executable LPFP-GlossBilliardAnimate + main-is: GlossBilliardAnimate.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + default-language: Haskell2010 + +executable LPFP-ProtonEB + main-is: ProtonEB.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , not-gloss + default-language: Haskell2010 + +executable LPFP-Hydrogen + main-is: Hydrogen.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , not-gloss + default-language: Haskell2010 + +executable LPFP-ElectricWave + main-is: ElectricWave.hs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + LPFP + , base >=4.7 && <5 + , diagrams-cairo, diagrams-lib + default-language: Haskell2010 diff --git a/app/ElectricWave.hs b/app/ElectricWave.hs new file mode 100644 --- /dev/null +++ b/app/ElectricWave.hs @@ -0,0 +1,10 @@ +{-# OPTIONS -Wall #-} + +import LPFP.Maxwell ( makeEpng, stateUpdate, jGaussian, initialStateFDTD ) +import Diagrams.Prelude ( black, yellow ) + +main :: IO () +main = let dt = 0.02e-9 -- 0.02 ns time step + numTimeSteps = 719 + in sequence_ $ map (makeEpng (yellow,black)) $ zip [0..numTimeSteps] $ + iterate (stateUpdate dt jGaussian) (initialStateFDTD 0.108) diff --git a/app/GlossAnimate.hs b/app/GlossAnimate.hs new file mode 100644 --- /dev/null +++ b/app/GlossAnimate.hs @@ -0,0 +1,24 @@ +{-# OPTIONS -Wall #-} + +import Graphics.Gloss + +displayMode :: Display +displayMode = InWindow "My Window" (1000, 700) (10, 10) + +disk :: Float -> Picture +disk radius = ThickCircle (radius / 2) radius + +redDisk :: Picture +redDisk = Color red (disk 25) + +projectileMotion :: Float -> Picture +projectileMotion t = Translate (xDisk t) (yDisk t) redDisk + +xDisk :: Float -> Float +xDisk t = 40 * t + +yDisk :: Float -> Float +yDisk t = 80 * t - 4.9 * t**2 + +main :: IO () +main = animate displayMode black projectileMotion diff --git a/app/GlossBilliard.hs b/app/GlossBilliard.hs new file mode 100644 --- /dev/null +++ b/app/GlossBilliard.hs @@ -0,0 +1,10 @@ +{-# OPTIONS -Wall #-} + +import LPFP.Mechanics3D ( simulateGloss ) +import LPFP.MultipleObjects (eulerCromerMPS ) + +import LPFP.MOExamples ( billiardInitial, billiardPicture, billiardUpdate ) + +main :: IO () +main = simulateGloss 1 100 billiardInitial billiardPicture + (billiardUpdate eulerCromerMPS 30) diff --git a/app/GlossBilliardAnimate.hs b/app/GlossBilliardAnimate.hs new file mode 100644 --- /dev/null +++ b/app/GlossBilliardAnimate.hs @@ -0,0 +1,8 @@ +{-# OPTIONS -Wall #-} + +import LPFP.MultipleObjects ( eulerCromerMPS ) +import LPFP.MOExamples + ( animateGloss, billiardPicture, billiardStates ) + +main :: IO () +main = animateGloss 1 billiardPicture (billiardStates eulerCromerMPS 30 0.01) diff --git a/app/GlossDisplay.hs b/app/GlossDisplay.hs new file mode 100644 --- /dev/null +++ b/app/GlossDisplay.hs @@ -0,0 +1,14 @@ +{-# OPTIONS -Wall #-} + +import Graphics.Gloss + +displayMode :: Display +displayMode = InWindow "Axes" (1000, 700) (10, 10) + +axes :: Picture +axes = Pictures [Color red $ Line [(0,0),(100, 0)] + ,Color green $ Line [(0,0),( 0,100)] + ] + +main :: IO () +main = display displayMode black axes diff --git a/app/GlossDisplay2.hs b/app/GlossDisplay2.hs new file mode 100644 --- /dev/null +++ b/app/GlossDisplay2.hs @@ -0,0 +1,23 @@ +{-# OPTIONS -Wall #-} + +import Graphics.Gloss + +displayMode :: Display +displayMode = InWindow "My Window" (1000, 700) (10, 10) + +blueCircle :: Picture +blueCircle = Color blue (Circle 100) + +disk :: Float -> Picture +disk radius = ThickCircle (radius / 2) radius + +redDisk :: Picture +redDisk = Color red (disk 100) + +wholePicture :: Picture +wholePicture = Pictures [Translate (-120) 0 blueCircle + ,Translate 120 0 redDisk + ] + +main :: IO () +main = display displayMode black wholePicture diff --git a/app/GlossHalley.hs b/app/GlossHalley.hs new file mode 100644 --- /dev/null +++ b/app/GlossHalley.hs @@ -0,0 +1,26 @@ +{-# OPTIONS -Wall #-} + +import LPFP.SimpleVec + ( xComp, yComp ) +import LPFP.Mechanics3D + ( ParticleState(..), simulateGloss, disk, halleyInitial, halleyUpdate ) +import Graphics.Gloss + ( Picture(..), pictures, translate, red, yellow ) + +diskComet :: Picture +diskComet = Color red (disk 10) + +diskSun :: Picture +diskSun = Color yellow (disk 20) + +halleyPicture :: ParticleState -> Picture +halleyPicture (ParticleState _m _q _t r _v) + = pictures [diskSun, translate xPixels yPixels diskComet] + where + pixelsPerMeter = 1e-10 + xPixels = pixelsPerMeter * realToFrac (xComp r) + yPixels = pixelsPerMeter * realToFrac (yComp r) + +main :: IO () +main = simulateGloss (365.25 * 24 * 60 * 60) 400 + halleyInitial halleyPicture halleyUpdate diff --git a/app/GlossProjectile.hs b/app/GlossProjectile.hs new file mode 100644 --- /dev/null +++ b/app/GlossProjectile.hs @@ -0,0 +1,28 @@ +{-# OPTIONS -Wall #-} + +import LPFP.SimpleVec + ( yComp, zComp ) +import LPFP.Mechanics3D + ( ParticleState(..), simulateGloss, disk + , projectileInitial, projectileUpdate ) +import Graphics.Gloss + ( Picture(..), red, scale, translate ) +import System.Environment + ( getArgs ) + +projectilePicture :: ParticleState -> Picture +projectilePicture (ParticleState _m _q _t r _v) + = scale 0.2 0.2 $ translate yFloat zFloat redDisk + where + yFloat = realToFrac (yComp r) + zFloat = realToFrac (zComp r) + redDisk :: Picture + redDisk = Color red (disk 50) + +mainWithArgs :: [String] -> IO () +mainWithArgs args + = simulateGloss 3 20 + (projectileInitial args) projectilePicture projectileUpdate + +main :: IO () +main = getArgs >>= mainWithArgs diff --git a/app/GlossProton.hs b/app/GlossProton.hs new file mode 100644 --- /dev/null +++ b/app/GlossProton.hs @@ -0,0 +1,10 @@ +{-# OPTIONS -Wall #-} + +import LPFP.Mechanics3D + ( simulateGloss + , twoProtInitial, twoProtPicture, twoProtUpdate + ) + +main :: IO () +main = simulateGloss 1e-8 20 + twoProtInitial twoProtPicture twoProtUpdate diff --git a/app/GlossSimulate.hs b/app/GlossSimulate.hs new file mode 100644 --- /dev/null +++ b/app/GlossSimulate.hs @@ -0,0 +1,31 @@ +{-# OPTIONS -Wall #-} + +import Graphics.Gloss + +displayMode :: Display +displayMode = InWindow "My Window" (1000, 700) (10, 10) + +-- updates per second of real time +rate :: Int +rate = 2 + +disk :: Float -> Picture +disk radius = ThickCircle (radius / 2) radius + +redDisk :: Picture +redDisk = Color red (disk 25) + +type State = (Float,Float) + +initialState :: State +initialState = (0,0) + +displayFunc :: State -> Picture +displayFunc (x,y) = Translate x y redDisk + +updateFunc :: Float -> State -> State +updateFunc dt (x,y) = (x + 10 * dt, y - 5 * dt) + +main :: IO () +main = simulate displayMode black rate initialState displayFunc + (\_ -> updateFunc) diff --git a/app/GlossSimulate2.hs b/app/GlossSimulate2.hs new file mode 100644 --- /dev/null +++ b/app/GlossSimulate2.hs @@ -0,0 +1,35 @@ +{-# OPTIONS -Wall #-} + +import Graphics.Gloss + +displayMode :: Display +displayMode = InWindow "My Window" (1000, 700) (10, 10) + +-- updates per second of real time +rate :: Int +rate = 24 + +disk :: Float -> Picture +disk radius = ThickCircle (radius / 2) radius + +redDisk :: Picture +redDisk = Color red (disk 25) + +type Position = (Float,Float) +type Velocity = (Float,Float) +type State = (Position,Velocity) + +initialState :: State +initialState = ((0,0),(40,80)) + +displayFunc :: State -> Picture +displayFunc ((x,y),_) = Translate x y redDisk + +updateFunc :: Float -> State -> State +updateFunc dt ((x,y),(vx,vy)) + = (( x + vx * dt, y + vy * dt) + ,(vx ,vy - 9.8 * dt)) + +main :: IO () +main = simulate displayMode black rate initialState displayFunc + (\_ -> updateFunc) diff --git a/app/GlossWave.hs b/app/GlossWave.hs new file mode 100644 --- /dev/null +++ b/app/GlossWave.hs @@ -0,0 +1,19 @@ +{-# OPTIONS -Wall #-} + +import LPFP.SimpleVec ( zeroV, iHat, (*^), xComp, yComp ) +import LPFP.Mechanics3D ( ParticleState(..), simulateGloss ) +import LPFP.MultipleObjects ( MultiParticleState(..) ) +import LPFP.MOExamples +import Graphics.Gloss ( Picture(..), scale, blue ) + +stringPicture :: MultiParticleState -> Picture +stringPicture (MPS sts) + = let rs = [zeroV] ++ [posVec st | st <- sts] ++ [0.65 *^ iHat] + xy r = (realToFrac $ xComp r, realToFrac $ yComp r) + xys = map xy rs + ppm = 400 -- pixels per meter + in scale ppm (20*ppm) $ Color blue $ Line xys + +main :: IO () +main = let initialState = stringInitialOvertone 3 + in simulateGloss 0.001 40 initialState stringPicture stringUpdate diff --git a/app/GnuplotWave.hs b/app/GnuplotWave.hs new file mode 100644 --- /dev/null +++ b/app/GnuplotWave.hs @@ -0,0 +1,28 @@ +{-# OPTIONS -Wall #-} + +import LPFP.SimpleVec ( R, zeroV, iHat, (*^), xComp, yComp ) +import LPFP.Mechanics3D ( ParticleState(..) ) +import LPFP.MultipleObjects ( MultiParticleState(..) ) +import LPFP.MOExamples +import Graphics.Gnuplot.Simple + +makePNG :: (Int,MultiParticleState) -> IO () +makePNG (n,MPS sts) + = let rs = [zeroV] ++ [posVec st | st <- sts] ++ [0.65 *^ iHat] + xy r = (xComp r, yComp r) + xys :: [(R,R)] + xys = map xy rs + threeDigitString = reverse $ take 3 $ reverse ("00" ++ show n) + pngFilePath = "GnuplotWave" ++ threeDigitString ++ ".png" + in plotPath [Title "Wave" + ,XLabel "Position (m)" + ,YLabel "Displacement (m)" + ,XRange (0,0.65) + ,YRange (-0.01,0.01) + ,PNG pngFilePath + ,Key Nothing + ] xys + +main :: IO () +main = sequence_ $ map makePNG $ zip [0..999] $ + iterate (stringUpdate 25e-6) (stringInitialOvertone 3) diff --git a/app/Hydrogen.hs b/app/Hydrogen.hs new file mode 100644 --- /dev/null +++ b/app/Hydrogen.hs @@ -0,0 +1,21 @@ +{-# OPTIONS -Wall #-} + +import LPFP.SimpleVec ( vec ) +import LPFP.Electricity ( elementaryCharge ) +import LPFP.CoordinateSystems ( cart ) +import LPFP.Charge ( protonOrigin ) +import LPFP.ElectricField ( eField, epsilon0 ) +import LPFP.Lorentz ( ParticleFieldState(..), animatePFS, defaultPFS ) + +main :: IO () +main = animatePFS period 30 (4*bohrRadius) + ( defaultPFS { mass = electronMass + , charge = -elementaryCharge -- electron charge + , position = cart bohrRadius 0 0 + , velocity = vec 0 v0 0 + , electricField = eField protonOrigin } ) + where electronMass = 9.109e-31 -- kg + bohrRadius = 0.529e-10 -- meters + v0 = elementaryCharge + / sqrt (4 * pi * epsilon0 * electronMass * bohrRadius) + period = 2 * pi * bohrRadius / v0 diff --git a/app/MakeTrajectoryGraph.hs b/app/MakeTrajectoryGraph.hs new file mode 100644 --- /dev/null +++ b/app/MakeTrajectoryGraph.hs @@ -0,0 +1,17 @@ +{-# OPTIONS -Wall #-} + +import LPFP.SimpleVec ( iHat, kHat, xComp, zComp, projectilePos, (^+^), (*^) ) +import Graphics.Gnuplot.Simple ( Attribute(..), plotPath ) + +main :: IO () +main = let posInitial = 10 *^ kHat + velInitial = 20 *^ cos (pi/6) *^ iHat ^+^ 20 *^ sin (pi/6) *^ kHat + posFunc = projectilePos posInitial velInitial + pairs = [(xComp r, zComp r) | t <- [0, 0.01 ..], let r = posFunc t] + plottingPairs = takeWhile (\(_,z) -> z >= 0) pairs + in plotPath [Title "Projectile Motion" + ,XLabel "Horizontal position (m)" + ,YLabel "Height of projectile (m)" + ,PNG "projectile.png" + ,Key Nothing + ] plottingPairs diff --git a/app/Mech3Proton.hs b/app/Mech3Proton.hs new file mode 100644 --- /dev/null +++ b/app/Mech3Proton.hs @@ -0,0 +1,6 @@ +{-# OPTIONS -Wall #-} + +import LPFP.Mechanics3D (simulateVis, protonInitial, protonPicture, protonUpdate) + +main :: IO () +main = simulateVis 1 60 protonInitial protonPicture protonUpdate diff --git a/app/ProtonEB.hs b/app/ProtonEB.hs new file mode 100644 --- /dev/null +++ b/app/ProtonEB.hs @@ -0,0 +1,13 @@ +{-# OPTIONS -Wall #-} + +import LPFP.SimpleVec ( vec ) +import LPFP.Electricity ( elementaryCharge ) +import LPFP.Lorentz ( ParticleFieldState(..), animatePFS, defaultPFS ) + +main :: IO () +main = animatePFS 1e-5 30 0.05 + ( defaultPFS { mass = 1.673e-27 -- proton in kg + , charge = elementaryCharge + , velocity = vec 0 2000 0 + , electricField = \_ -> vec 0 20 0 + , magneticField = \_ -> vec 0 0 0.01 } ) diff --git a/app/VisAnimate.hs b/app/VisAnimate.hs new file mode 100644 --- /dev/null +++ b/app/VisAnimate.hs @@ -0,0 +1,13 @@ +{-# OPTIONS -Wall #-} + +import Vis +import SpatialMath + +rotatingCube :: Float -> VisObject Float +rotatingCube t = RotEulerRad (Euler 0 0 t) (Cube 1 Solid blue) + +orient :: VisObject Float -> VisObject Float +orient pict = RotEulerDeg (Euler 270 180 0) $ pict + +main :: IO () +main = animate defaultOpts (orient . rotatingCube) diff --git a/app/VisDisplay.hs b/app/VisDisplay.hs new file mode 100644 --- /dev/null +++ b/app/VisDisplay.hs @@ -0,0 +1,11 @@ +{-# OPTIONS -Wall #-} + +import Vis + +type R = Double + +blueCube :: VisObject R +blueCube = Cube 1 Solid blue + +main :: IO () +main = display defaultOpts blueCube diff --git a/app/VisDisplay2.hs b/app/VisDisplay2.hs new file mode 100644 --- /dev/null +++ b/app/VisDisplay2.hs @@ -0,0 +1,15 @@ +{-# OPTIONS -Wall #-} + +import Vis +import Linear + +type R = Double + +axes :: VisObject R +axes = VisObjects [Line Nothing [V3 0 0 0, V3 1 0 0] red + ,Line Nothing [V3 0 0 0, V3 0 1 0] green + ,Line Nothing [V3 0 0 0, V3 0 0 1] blue + ] + +main :: IO () +main = display defaultOpts axes diff --git a/app/VisDisplay3.hs b/app/VisDisplay3.hs new file mode 100644 --- /dev/null +++ b/app/VisDisplay3.hs @@ -0,0 +1,19 @@ +{-# OPTIONS -Wall #-} + +import Vis +import Linear +import SpatialMath + +type R = Double + +axes :: VisObject R +axes = VisObjects [Line Nothing [V3 0 0 0, V3 1 0 0] red + ,Line Nothing [V3 0 0 0, V3 0 1 0] green + ,Line Nothing [V3 0 0 0, V3 0 0 1] blue + ] + +orient :: VisObject R -> VisObject R +orient pict = RotEulerDeg (Euler 270 180 0) $ pict + +main :: IO () +main = display defaultOpts (orient axes) diff --git a/app/VisProjectile.hs b/app/VisProjectile.hs new file mode 100644 --- /dev/null +++ b/app/VisProjectile.hs @@ -0,0 +1,23 @@ +{-# OPTIONS -Wall #-} + +import LPFP.SimpleVec ( R, (*^) ) +import LPFP.Mechanics3D + ( ParticleState(..), simulateVis + , projectileInitial, projectileUpdate, v3FromVec ) +import Vis + ( VisObject(..), Flavour(..), red ) +import System.Environment + ( getArgs ) + +projectileVisObject :: ParticleState -> VisObject R +projectileVisObject st + = let r = posVec st + in Trans (v3FromVec (0.01 *^ r)) (Sphere 0.1 Solid red) + +mainWithArgs :: [String] -> IO () +mainWithArgs args + = simulateVis 3 20 + (projectileInitial args) projectileVisObject projectileUpdate + +main :: IO () +main = getArgs >>= mainWithArgs diff --git a/app/VisSimulate.hs b/app/VisSimulate.hs new file mode 100644 --- /dev/null +++ b/app/VisSimulate.hs @@ -0,0 +1,19 @@ +{-# OPTIONS -Wall #-} + +import Vis + +type State = (Int,[Float]) + +-- seconds / update +dt :: Double +dt = 0.5 + +displayFunc :: State -> VisObject Double +displayFunc (n,ts) = Text2d (show n ++ " " ++ show (take 4 ts)) + (100,100) Fixed9By15 orange + +updateFunc :: Float -> State -> State +updateFunc t (n,ts) = (n+1,t:ts) + +main :: IO () +main = simulate defaultOpts dt (0,[]) displayFunc updateFunc diff --git a/app/VisTwoSprings.hs b/app/VisTwoSprings.hs new file mode 100644 --- /dev/null +++ b/app/VisTwoSprings.hs @@ -0,0 +1,24 @@ +{-# OPTIONS -Wall #-} + +import LPFP.SimpleVec ( R, zeroV ) +import LPFP.Mechanics3D ( posVec, simulateVis, v3FromVec ) +import LPFP.MultipleObjects ( MultiParticleState(..) ) +import LPFP.MOExamples ( twoSpringsInitial, twoSpringsUpdate ) +import Vis ( VisObject(..), Flavour(..), red, green, blue ) + +main :: IO () +main = simulateVis 1 20 twoSpringsInitial twoSpringsVisObject twoSpringsUpdate + +twoSpringsVisObject :: MultiParticleState -> VisObject R +twoSpringsVisObject (MPS sts) + = let r0 = posVec (sts !! 0) + r1 = posVec (sts !! 1) + springsObj = Line Nothing [v3FromVec zeroV + ,v3FromVec r0 + ,v3FromVec r1] blue + objs = [Trans (v3FromVec r0) (Sphere 0.1 Solid red) + ,Trans (v3FromVec r1) (Sphere 0.1 Solid green) + ,springsObj + ] + vpm = 1 -- Vis units per meter + in Scale (vpm,vpm,vpm) $ VisObjects objs diff --git a/app/hello.hs b/app/hello.hs new file mode 100644 --- /dev/null +++ b/app/hello.hs @@ -0,0 +1,2 @@ +main :: IO () +main = putStrLn "Hello, world!" diff --git a/src/LPFP.hs b/src/LPFP.hs new file mode 100644 --- /dev/null +++ b/src/LPFP.hs @@ -0,0 +1,448 @@ +{-# OPTIONS -Wall #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : LPFP +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from the book Learn Physics with Functional Programming +-} + +module LPFP + ( + -- * (Approximations to) Real numbers + R + , Time + -- * Vectors + , Vec + , PosVec + , Velocity + , Acceleration + , vec + , (^+^) + , (^-^) + , (*^) + , (^*) + , (^/) + , (<.>) + , (><) + , magnitude + , zeroV + , negateV + , sumV + , xComp + , yComp + , zComp + , iHat + , jHat + , kHat + , positionCV + , velocityCA + , positionCA + , aParallel + , aPerp + , speedRateChange + -- * Calculus + , Derivative + , VecDerivative + , derivative + , vecDerivative + , integral + , antiDerivative + , velFromPos + , accFromVel + -- ** Differential equations + , UpdateFunction + , DifferentialEquation + , NumericalMethod + , RealVectorSpace(..) + , Diff(..) + , solver + , euler + , rungeKutta4 + -- * 3D Mechanics + -- ** Single particle state + , ParticleState(..) + , DParticleState(..) + , HasTime(..) + , defaultParticleState + , newtonSecondPS + , relativityPS + , eulerCromerPS + , statesPS + , updatePS + -- ** One-body forces + , OneBodyForce + , earthSurfaceGravity + , sunGravity + , airResistance + , windForce + , uniformLorentzForce + , fixedLinearSpring + -- ** Graphics utilities + , simulateGloss + , simulateVis + , v3FromVec + , orient + -- * Interacting particles + , Force(..) + , MultiParticleState(..) + , DMultiParticleState(..) + -- ** Two-body forces + , TwoBodyForce + , universalGravity + , linearSpring + , centralForce + , billiardForce + , newtonSecondMPS + , eulerCromerMPS + , updateMPS + , statesMPS + , Justification(..) + , Table(..) + , kineticEnergy + , systemKE + , momentum + , systemP + , linearSpringPE + , earthSurfaceGravityPE + , tenths + , sigFigs + , animateGloss + , animateVis + -- * Electricity + , elementaryCharge + , coulombForce + -- * Coordinate Systems + , Position + , Displacement + , ScalarField + , VectorField + , CoordinateSystem + , cartesian + , cylindrical + , spherical + , cart + , cyl + , sph + , cartesianCoordinates + , cylindricalCoordinates + , sphericalCoordinates + , displacement + , shiftPosition + , rHat + , thetaHat + , phiHat + , sHat + , xHat + , yHat + , zHat + , origin + , xSF + , ySF + , rSF + , rVF + , fst3 + , snd3 + , thd3 + , addScalarFields + , addVectorFields + , sf3D + , vf3D + , v3FromPos + , sfTable + , vfPNG + , vfPNGxy + , vfGrad + -- * Geometry + , Curve(..) + , unitCircle + , straightLine + , Surface(..) + , unitSphere + , centeredSphere + , sphere + , northernHemisphere + , disk + , shiftSurface + , Volume(..) + , unitBall + , centeredBall + , northernHalfBall + , centeredCylinder + -- * Electromagnetic Theory + -- ** Charge + , Charge + , ChargeDistribution(..) + , totalCharge + , electricDipoleMoment + -- ** Electric Field + , epsilon0 + , cSI + , mu0 + , eField + , ScalarLineIntegral + , ScalarSurfaceIntegral + , ScalarVolumeIntegral + , VectorLineIntegral + , VectorSurfaceIntegral + , VectorVolumeIntegral + , CurveApprox + , SurfaceApprox + , VolumeApprox + , scalarLineIntegral + , scalarSurfaceIntegral + , scalarVolumeIntegral + , vectorLineIntegral + , vectorSurfaceIntegral + , vectorVolumeIntegral + , dottedLineIntegral + , dottedSurfaceIntegral + , curveSample + , surfaceSample + , volumeSample + , Field + -- ** Current + , Current + , CurrentDistribution(..) + , crossedLineIntegral + , totalCurrent + , magneticDipoleMoment + -- ** Magnetic Field + , bField + -- ** Lorentz Force Law + , lorentzForce + , newtonSecondPFS + , defaultPFS + , pfsVisObject + , animatePFS + -- ** Maxwell Equations + , directionalDerivative + , curl + , FieldState + ) + where + +import LPFP.SimpleVec + ( R + , Vec + , Time + , PosVec + , Velocity + , Acceleration + , Derivative + , VecDerivative + , vec + , (^+^) + , (^-^) + , (*^) + , (^*) + , (^/) + , (<.>) + , (><) + , magnitude + , zeroV + , negateV + , sumV + , xComp + , yComp + , zComp + , iHat + , jHat + , kHat + , positionCV + , velocityCA + , positionCA + , derivative + , vecDerivative + , velFromPos + , accFromVel + , aParallel + , aPerp + , speedRateChange + ) +import LPFP.Newton2 + ( integral + , antiDerivative + ) +import LPFP.Mechanics1D + ( UpdateFunction + , DifferentialEquation + , NumericalMethod + , RealVectorSpace(..) + , Diff(..) + , solver + , euler + , rungeKutta4 + ) +import LPFP.Mechanics3D + ( ParticleState(..) + , DParticleState(..) + , HasTime(..) + , OneBodyForce + , defaultParticleState + , newtonSecondPS + , relativityPS + , earthSurfaceGravity + , sunGravity + , airResistance + , windForce + , uniformLorentzForce + , eulerCromerPS + , statesPS + , updatePS + , simulateGloss + , simulateVis + , v3FromVec + , orient + ) +import LPFP.MultipleObjects + ( TwoBodyForce + , Force(..) + , MultiParticleState(..) + , DMultiParticleState(..) + , universalGravity + , linearSpring + , fixedLinearSpring + , centralForce + , billiardForce + , newtonSecondMPS + , eulerCromerMPS + , updateMPS + , statesMPS + ) +import LPFP.MOExamples + ( Justification(..) + , Table(..) + , kineticEnergy + , systemKE + , momentum + , systemP + , linearSpringPE + , earthSurfaceGravityPE + , tenths + , sigFigs + , animateGloss + , animateVis + ) +import LPFP.Electricity + ( Charge + , elementaryCharge + , coulombForce + ) +import LPFP.CoordinateSystems + ( Position + , Displacement + , ScalarField + , VectorField + , CoordinateSystem + , cartesian + , cylindrical + , spherical + , cart + , cyl + , sph + , cartesianCoordinates + , cylindricalCoordinates + , sphericalCoordinates + , displacement + , shiftPosition + , rHat + , thetaHat + , phiHat + , sHat + , xHat + , yHat + , zHat + , origin + , xSF + , ySF + , rSF + , rVF + , fst3 + , snd3 + , thd3 + , addScalarFields + , addVectorFields + , sf3D + , vf3D + , v3FromPos + , sfTable + , vfPNG + , vfPNGxy + , vfGrad + ) +import LPFP.Geometry + ( Curve(..) + , unitCircle + , straightLine + , Surface(..) + , unitSphere + , centeredSphere + , sphere + , northernHemisphere + , disk + , shiftSurface + , Volume(..) + , unitBall + , centeredBall + , northernHalfBall + , centeredCylinder + ) +import LPFP.Charge + ( ChargeDistribution(..) + , totalCharge + , electricDipoleMoment + ) +import LPFP.ElectricField + ( epsilon0 + , cSI + , mu0 + , eField + , ScalarLineIntegral + , ScalarSurfaceIntegral + , ScalarVolumeIntegral + , VectorLineIntegral + , VectorSurfaceIntegral + , VectorVolumeIntegral + , CurveApprox + , SurfaceApprox + , VolumeApprox + , scalarLineIntegral + , scalarSurfaceIntegral + , scalarVolumeIntegral + , vectorLineIntegral + , vectorSurfaceIntegral + , vectorVolumeIntegral + , dottedLineIntegral + , dottedSurfaceIntegral + , curveSample + , surfaceSample + , volumeSample + , Field + ) +import LPFP.Current + ( Current + , CurrentDistribution(..) + , crossedLineIntegral + , totalCurrent + , magneticDipoleMoment + ) +import LPFP.MagneticField + ( bField + ) +import LPFP.Lorentz + ( lorentzForce + , newtonSecondPFS + , defaultPFS + , pfsVisObject + , animatePFS + ) +import LPFP.Maxwell + ( directionalDerivative + , curl + , FieldState + ) diff --git a/src/LPFP/Charge.hs b/src/LPFP/Charge.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/Charge.hs @@ -0,0 +1,123 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFP.Charge +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 24 of the book Learn Physics with Functional Programming +-} + +module LPFP.Charge where + +import LPFP.SimpleVec ( R, Vec, vec, sumV, (*^), (^/), (<.>), magnitude, negateV ) +import LPFP.Electricity ( elementaryCharge ) +import LPFP.CoordinateSystems ( Position, ScalarField, origin, cart, sph + , rVF, displacement, shiftPosition ) +import LPFP.Geometry ( Curve(..), Surface(..), Volume(..) + , straightLine, shiftSurface, disk ) +import LPFP.Integrals + ( scalarLineIntegral, scalarSurfaceIntegral, scalarVolumeIntegral + , vectorLineIntegral, vectorSurfaceIntegral, vectorVolumeIntegral + , curveSample, surfaceSample, volumeSample ) + +type Charge = R + +data ChargeDistribution + = PointCharge Charge Position + | LineCharge ScalarField Curve + | SurfaceCharge ScalarField Surface + | VolumeCharge ScalarField Volume + | MultipleCharges [ChargeDistribution] + +protonOrigin :: ChargeDistribution +protonOrigin = PointCharge elementaryCharge origin + +chargedLine :: Charge -> R -> ChargeDistribution +chargedLine q len + = LineCharge (const $ q / len) $ + Curve (\z -> cart 0 0 z) (-len/2) (len/2) + +chargedBall :: Charge -> R -> ChargeDistribution +chargedBall q radius + = VolumeCharge (const $ q / (4/3*pi*radius**3)) $ + Volume (\(r,theta,phi) -> sph r theta phi) + 0 radius (const 0) (const pi) (\_ _ -> 0) (\_ _ -> 2*pi) + +diskCap :: R -> R -> R -> ChargeDistribution +diskCap radius plateSep sigma + = MultipleCharges + [SurfaceCharge (const sigma) $ + shiftSurface (vec 0 0 (plateSep/2)) (disk radius) + ,SurfaceCharge (const $ -sigma) $ + shiftSurface (vec 0 0 (-plateSep/2)) (disk radius) + ] + +totalCharge :: ChargeDistribution -> Charge +totalCharge (PointCharge q _) + = q +totalCharge (LineCharge lambda c) + = scalarLineIntegral (curveSample 1000) lambda c +totalCharge (SurfaceCharge sigma s) + = scalarSurfaceIntegral (surfaceSample 200) sigma s +totalCharge (VolumeCharge rho v) + = scalarVolumeIntegral (volumeSample 50) rho v +totalCharge (MultipleCharges ds ) + = sum [totalCharge d | d <- ds] + +simpleDipole :: Vec -- electric dipole moment + -> R -- charge separation + -> ChargeDistribution +simpleDipole p sep + = let q = magnitude p / sep + disp = (sep/2) *^ (p ^/ magnitude p) + in MultipleCharges + [PointCharge q (shiftPosition disp origin) + ,PointCharge (-q) (shiftPosition (negateV disp) origin) + ] + +electricDipoleMoment :: ChargeDistribution -> Vec +electricDipoleMoment (PointCharge q r) + = q *^ displacement origin r +electricDipoleMoment (LineCharge lambda c) + = vectorLineIntegral (curveSample 1000) (\r -> lambda r *^ rVF r) c +electricDipoleMoment (SurfaceCharge sigma s) + = vectorSurfaceIntegral (surfaceSample 200) (\r -> sigma r *^ rVF r) s +electricDipoleMoment (VolumeCharge rho v) + = vectorVolumeIntegral (volumeSample 50) (\r -> rho r *^ rVF r) v +electricDipoleMoment (MultipleCharges ds ) + = sumV [electricDipoleMoment d | d <- ds] + +lineDipole :: Vec -- dipole moment + -> R -- charge separation + -> ChargeDistribution +lineDipole p sep + = let disp = (sep/2) *^ (p ^/ magnitude p) + curve = straightLine (shiftPosition (negateV disp) origin) + (shiftPosition disp origin) + coeff = 12 / sep**3 + lambda r = coeff * (displacement origin r <.> p) + in LineCharge lambda curve + +chargedDisk :: Charge -> R -> ChargeDistribution +chargedDisk q radius = undefined q radius + +circularLineCharge :: Charge -> R -> ChargeDistribution +circularLineCharge q radius = undefined q radius + +chargedSquarePlate :: Charge -> R -> ChargeDistribution +chargedSquarePlate q side = undefined q side + +chargedSphericalShell :: Charge -> R -> ChargeDistribution +chargedSphericalShell q radius = undefined q radius + +chargedCube :: Charge -> R -> ChargeDistribution +chargedCube q side = undefined q side + +squareCap :: R -> R -> R -> ChargeDistribution +squareCap side plateSep sigma = undefined side plateSep sigma + +hydrogen :: ChargeDistribution +hydrogen = undefined diff --git a/src/LPFP/CoordinateSystems.hs b/src/LPFP/CoordinateSystems.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/CoordinateSystems.hs @@ -0,0 +1,301 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFP.CoordinateSystems +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 22 of the book Learn Physics with Functional Programming +-} + +module LPFP.CoordinateSystems where + +import LPFP.SimpleVec + ( R, Vec, (^/), vec, xComp, yComp, zComp, iHat, jHat, kHat + , magnitude, sumV, zeroV ) +import LPFP.Mechanics3D ( orient, v3FromVec ) +import LPFP.MOExamples ( Table(..), Justification(..) ) +import qualified Vis as V +import SpatialMath ( V3(..) ) +import Diagrams.Prelude + ( Diagram, V2(..), PolyType(..), PolyOrientation(..), PolygonOpts(..) + , (#), (@@), dims, p2, r2, arrowAt, position, fc, black, white + , blend, none, lw, rotate, deg, rad, scale, polygon, sinA ) +import Diagrams.Backend.Cairo ( B, renderCairo ) + +data Position = Cart R R R + deriving (Show) + +type CoordinateSystem = (R,R,R) -> Position + +cartesian :: CoordinateSystem +cartesian (x,y,z) + = Cart x y z + +cylindrical :: CoordinateSystem +cylindrical (s,phi,z) + = Cart (s * cos phi) (s * sin phi) z + +spherical :: CoordinateSystem +spherical (r,theta,phi) + = Cart (r * sin theta * cos phi) + (r * sin theta * sin phi) + (r * cos theta) + +cart :: R -- x coordinate + -> R -- y coordinate + -> R -- z coordinate + -> Position +cart = Cart + +cyl :: R -- s coordinate + -> R -- phi coordinate + -> R -- z coordinate + -> Position +cyl s phi z = cylindrical (s,phi,z) + +sph :: R -- r coordinate + -> R -- theta coordinate + -> R -- phi coordinate + -> Position +sph r theta phi = spherical (r,theta,phi) + +origin :: Position +origin = cart 0 0 0 + +cartesianCoordinates :: Position -> (R,R,R) +cartesianCoordinates (Cart x y z) = (x,y,z) + +cylindricalCoordinates :: Position -> (R,R,R) +cylindricalCoordinates (Cart x y z) = (s,phi,z) + where + s = sqrt(x**2 + y**2) + phi = atan2 y x + +sphericalCoordinates :: Position -> (R,R,R) +sphericalCoordinates (Cart x y z) = (r,theta,phi) + where + r = sqrt(x**2 + y**2 + z**2) + theta = atan2 s z + s = sqrt(x**2 + y**2) + phi = atan2 y x + +type Displacement = Vec + +displacement :: Position -- source position + -> Position -- target position + -> Displacement +displacement (Cart x' y' z') (Cart x y z) + = vec (x-x') (y-y') (z-z') + +shiftPosition :: Displacement -> Position -> Position +shiftPosition v (Cart x y z) + = Cart (x + xComp v) (y + yComp v) (z + zComp v) + +type ScalarField = Position -> R + +xSF :: ScalarField +xSF p = x + where + (x,_,_) = cartesianCoordinates p + +rSF :: ScalarField +rSF p = r + where + (r,_,_) = sphericalCoordinates p + +fst3 :: (a,b,c) -> a +fst3 (u,_,_) = u + +snd3 :: (a,b,c) -> b +snd3 (_,u,_) = u + +thd3 :: (a,b,c) -> c +thd3 (_,_,u) = u + +ySF :: ScalarField +ySF = snd3 . cartesianCoordinates + +type VectorField = Position -> Vec + +sHat :: VectorField +sHat r = vec ( cos phi) (sin phi) 0 + where + (_,phi,_) = cylindricalCoordinates r + +phiHat :: VectorField +phiHat r = vec (-sin phi) (cos phi) 0 + where + (_,phi,_) = cylindricalCoordinates r + +rHat :: VectorField +rHat rv = let d = displacement origin rv + in if d == zeroV + then zeroV + else d ^/ magnitude d + +thetaHat :: VectorField +thetaHat r = vec ( cos theta * cos phi) + ( cos theta * sin phi) + (-sin theta ) + where + (_,theta,phi) = sphericalCoordinates r + +xHat :: VectorField +xHat = const iHat + +yHat :: VectorField +yHat = const jHat + +zHat :: VectorField +zHat = const kHat + +rVF :: VectorField +rVF = displacement origin + +addScalarFields :: [ScalarField] -> ScalarField +addScalarFields flds r = sum [fld r | fld <- flds] + +addVectorFields :: [VectorField] -> VectorField +addVectorFields flds r = sumV [fld r | fld <- flds] + +sf3D :: [Position] -- positions to use + -> ScalarField -- to display + -> IO () +sf3D ps sf + = V.display whiteBackground $ orient $ + V.VisObjects [V.Text3d (show (round $ sf p :: Int)) + (v3FromPos p) V.Fixed9By15 V.black + | p <- ps] + +v3FromPos :: Position -> V3 R +v3FromPos p = V3 x y z + where + (x,y,z) = cartesianCoordinates p + +whiteBackground :: V.Options +whiteBackground = V.defaultOpts {V.optBackgroundColor = Just V.white} + +whiteBackground' :: V.Options +whiteBackground' + = V.defaultOpts {V.optBackgroundColor = Just V.white, + V.optInitialCamera = Just V.Camera0 {V.rho0 = 40.0, + V.theta0 = 45.0, + V.phi0 = 20.0}} + +ySF3D :: IO () +ySF3D = sf3D [cart x y z | x <- [-6,-2..6] + , y <- [-6,-2..6] + , z <- [-6,-2..6]] ySF + +sfTable :: ((R,R) -> Position) + -> [R] -- horizontal + -> [R] -- vertical + -> ScalarField + -> Table Int +sfTable toPos ss ts sf + = Table RJ [[round $ sf $ toPos (s,t) | s <- ss] | t <- reverse ts] + +vf3D :: R -- scale factor, vector field units per meter + -> [Position] -- positions to show the field + -> VectorField -- vector field to display + -> IO () +vf3D unitsPerMeter ps vf + = V.display whiteBackground $ orient $ + V.VisObjects [V.Trans (v3FromPos p) $ + visVec V.black (vf p ^/ unitsPerMeter) + | p <- ps] + +visVec :: V.Color -> Vec -> V.VisObject R +visVec color v = let vmag = magnitude v + in V.Arrow (vmag,20*vmag) (v3FromVec v) color + +phiHat3D :: IO () +phiHat3D = vf3D 1 [cyl r ph z | r <- [1,2,3] + , ph <- [0,pi/4..2*pi] + , z <- [-2..2]] phiHat + +vfPNG :: ((R,R) -> Position) + -> (Vec -> (R,R)) + -> FilePath -- file name + -> R -- scale factor in units per meter + -> [(R,R)] -- positions to use + -> VectorField + -> IO () +vfPNG toPos fromVec fileName unitsPerMeter pts vf + = let vf2d = r2 . fromVec . (^/ unitsPerMeter) . vf . toPos + pic = mconcat [arrowAt (p2 pt) (vf2d pt) | pt <- pts] + in renderCairo fileName (dims (V2 1024 1024)) pic + +vfPNGxy :: FilePath -- file name + -> R -- scale factor + -> [(R,R)] -- positions to use + -> VectorField + -> IO () +vfPNGxy = vfPNG (\(x,y) -> cart x y 0) (\v -> (xComp v, yComp v)) + +phiHatPNG :: IO () +phiHatPNG + = vfPNGxy "phiHatPNG.png" 1 + [(r * cos ph, r * sin ph) | r <- [1,2] + , ph <- [0,pi/4..2*pi]] phiHat + +rVFpng :: IO () +rVFpng + = vfPNGxy "rVFpng.png" 2 + [(r * cos ph, r * sin ph) | r <- [1,2] + , ph <- [0,pi/4..2*pi]] rVF + +vfGrad :: (R -> R) + -> ((R,R) -> Position) + -> (Vec -> (R,R)) + -> FilePath + -> Int -- n for n x n + -> VectorField + -> IO () +vfGrad curve toPos fromVec fileName n vf + = let step = 2 / fromIntegral n + xs = [-1+step/2, -1+3*step/2 .. 1-step/2] + pts = [(x, y) | x <- xs, y <- xs] + array = [(pt,magRad $ fromVec $ vf $ toPos pt) | pt <- pts] + maxMag = maximum (map (fst . snd) array) + scaledArrow m th = scale step $ arrowMagRad (curve (m/maxMag)) th + pic = position [(p2 pt, scaledArrow m th) | (pt,(m,th)) <- array] + in renderCairo fileName (dims (V2 1024 1024)) pic + +magRad :: (R,R) -> (R,R) +magRad (x,y) = (sqrt (x*x + y*y), atan2 y x) + +-- magnitude from 0 to 1 +arrowMagRad :: R -- magnitude + -> R -- angle in radians, counterclockwise from x axis + -> Diagram B +arrowMagRad mag th + = let r = sinA (15 @@ deg) / sinA (60 @@ deg) + myType = PolyPolar [120 @@ deg, 0 @@ deg, 45 @@ deg, 30 @@ deg, + 45 @@ deg, 0 @@ deg, 120 @@ deg] + [1,1,r,1,1,r,1,1] + myOpts = PolygonOpts myType NoOrient (p2 (0,0)) + in scale 0.5 $ polygon myOpts # lw none # fc (blend mag black white) # + rotate (th @@ rad) + +rVFGrad :: IO () +rVFGrad = vfGrad id + (\(x,y) -> cart x y 0) + (\v -> (xComp v,yComp v)) + "rVFGrad.png" 20 + rVF + +thetaSF :: ScalarField +thetaSF = undefined + +thetaHat3D :: IO () +thetaHat3D = undefined + +thetaHatGrad :: IO () +thetaHatGrad = vfGrad id undefined undefined "thetaHatGrad.png" 20 thetaHat + +phiHatGrad :: IO () +phiHatGrad = undefined diff --git a/src/LPFP/Current.hs b/src/LPFP/Current.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/Current.hs @@ -0,0 +1,102 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFP.Current +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 26 of the book Learn Physics with Functional Programming +-} + +module LPFP.Current where + +import LPFP.SimpleVec + ( R, Vec, sumV, (><), (*^) ) +import LPFP.CoordinateSystems + ( VectorField, rVF, cyl, phiHat ) +import LPFP.Geometry + ( Curve(..), Surface(..), Volume(..) ) +import LPFP.ElectricField + ( CurveApprox, curveSample, surfaceSample, volumeSample + , vectorSurfaceIntegral, vectorVolumeIntegral ) + +type Current = R + +data CurrentDistribution + = LineCurrent Current Curve + | SurfaceCurrent VectorField Surface + | VolumeCurrent VectorField Volume + | MultipleCurrents [CurrentDistribution] + +circularCurrentLoop :: R -- radius + -> R -- current + -> CurrentDistribution +circularCurrentLoop radius i + = LineCurrent i (Curve (\phi -> cyl radius phi 0) 0 (2*pi)) + +wireSolenoid :: R -- radius + -> R -- length + -> R -- turns/length + -> R -- current + -> CurrentDistribution +wireSolenoid radius len n i + = LineCurrent i (Curve (\phi -> cyl radius phi (phi/(2*pi*n))) + (-pi*n*len) (pi*n*len)) + +sheetSolenoid :: R -- radius + -> R -- length + -> R -- turns/length + -> R -- current + -> CurrentDistribution +sheetSolenoid radius len n i + = SurfaceCurrent (\r -> (n*i) *^ phiHat r) + (Surface (\(phi,z) -> cyl radius phi z) + 0 (2*pi) (const $ -len/2) (const $ len/2)) + +wireToroid :: R -- small radius + -> R -- big radius + -> R -- number of turns + -> R -- current + -> CurrentDistribution +wireToroid smallR bigR n i + = let alpha phi = n * phi + curve phi = cyl (bigR + smallR * cos (alpha phi)) phi + (smallR * sin (alpha phi)) + in LineCurrent i (Curve curve 0 (2*pi)) + +crossedLineIntegral :: CurveApprox -> VectorField -> Curve -> Vec +crossedLineIntegral approx vF c + = sumV [vF r' >< dl' | (r',dl') <- approx c] + +magneticDipoleMoment :: CurrentDistribution -> Vec +magneticDipoleMoment (LineCurrent i c) + = crossedLineIntegral (curveSample 1000) (\r -> 0.5 *^ i *^ rVF r) c +magneticDipoleMoment (SurfaceCurrent k s) + = vectorSurfaceIntegral (surfaceSample 200) (\r -> 0.5 *^ (rVF r >< k r)) s +magneticDipoleMoment (VolumeCurrent j v) + = vectorVolumeIntegral (volumeSample 50) (\r -> 0.5 *^ (rVF r >< j r)) v +magneticDipoleMoment (MultipleCurrents ds ) + = sumV [magneticDipoleMoment d | d <- ds] + +helmholtzCoil :: R -- radius + -> R -- current + -> CurrentDistribution +helmholtzCoil radius i = undefined radius i + +longStraightWire :: R -- wire length + -> R -- current + -> CurrentDistribution +longStraightWire len i = undefined len i + +torus :: R -> R -> Surface +torus smallR bigR + = Surface (\(phi,alpha) -> cyl (bigR + smallR * cos alpha) phi + (smallR * sin alpha)) + 0 (2*pi) (const 0) (const $ 2*pi) + +totalCurrent :: VectorField -- volume current density + -> Surface + -> Current -- total current through surface +totalCurrent j s = undefined j s diff --git a/src/LPFP/ElectricField.hs b/src/LPFP/ElectricField.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/ElectricField.hs @@ -0,0 +1,316 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFP.ElectricField +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 25 of the book Learn Physics with Functional Programming +-} + +module LPFP.ElectricField where + +import LPFP.SimpleVec + ( R, Vec, (^+^), (^-^), (*^), (^*), (^/), (<.>), (><) + , sumV, magnitude, vec, xComp, yComp, zComp, kHat ) +import LPFP.CoordinateSystems + ( Position, ScalarField, VectorField + , displacement, shiftPosition, addVectorFields + , cart, sph, vf3D, vfPNGxy, vfGrad, origin, rVF ) +import LPFP.Geometry ( Curve(..), Surface(..), Volume(..) ) +import LPFP.Charge + ( Charge, ChargeDistribution(..) + , diskCap, protonOrigin, simpleDipole, lineDipole ) + +epsilon0 :: R +epsilon0 = 1/(mu0 * cSI**2) + +cSI :: R +cSI = 299792458 -- m/s + +mu0 :: R +mu0 = 4e-7 * pi -- N/A^2 + +eFieldFromPointCharge + :: Charge -- in Coulombs + -> Position -- of point charge (in m) + -> VectorField -- electric field (in V/m) +eFieldFromPointCharge q1 r1 r + = let k = 1 / (4 * pi * epsilon0) + d = displacement r1 r + in (k * q1) *^ d ^/ magnitude d ** 3 + +eField :: ChargeDistribution -> VectorField +eField (PointCharge q r) = eFieldFromPointCharge q r +eField (LineCharge lam c) = eFieldFromLineCharge lam c +eField (SurfaceCharge sig s) = eFieldFromSurfaceCharge sig s +eField (VolumeCharge rho v) = eFieldFromVolumeCharge rho v +eField (MultipleCharges cds) = addVectorFields $ map eField cds + +eFieldPicProton2D :: IO () +eFieldPicProton2D + = vfPNGxy "eFieldPicProton2D.png" 3e-9 pts (eField protonOrigin) + where + pts = [(r * cos th, r * sin th) | r <- [1,1.5,2] + , th <- [0,pi/4 .. 2*pi]] + +eFieldPicProtonGrad :: IO () +eFieldPicProtonGrad + = vfGrad (**0.2) (\(x,y) -> cart x y 0) (\v -> (xComp v, yComp v)) + "eFieldPicProtonGrad.png" 20 (eField protonOrigin) + +eFieldPicProton3D :: IO () +eFieldPicProton3D = vf3D 4e-9 + [sph r th ph | r <- [1,1.5,2] + , th <- [0,pi/4..pi] + , ph <- [0,pi/4..2*pi]] (eField protonOrigin) + +simpleDipoleSodiumChloride :: ChargeDistribution +simpleDipoleSodiumChloride = simpleDipole (vec 0 0 2.99e-29) 2.36e-10 + +eFieldSodiumChloride :: VectorField +eFieldSodiumChloride = eField simpleDipoleSodiumChloride + +eFieldPicSimpleDipole :: IO () +eFieldPicSimpleDipole + = vfGrad (**0.2) (\(y,z) -> cart 0 (3e-10*y) (3e-10*z)) + (\v -> (yComp v, zComp v)) "eFieldPicSimpleDipole.png" 20 + eFieldSodiumChloride + +eFieldIdealDipole :: Vec -- electric dipole moment + -> VectorField -- electric field +eFieldIdealDipole p r + = let k = 1 / (4 * pi * epsilon0) -- SI units + rMag = magnitude (rVF r) + rUnit = rVF r ^/ rMag + in k *^ (1 / rMag**3) *^ (3 *^ (p <.> rUnit) *^ rUnit ^-^ p) + +eFieldPicIdealDipole :: IO () +eFieldPicIdealDipole + = vfGrad (**0.2) (\(y,z) -> cart 0 (3e-10*y) (3e-10*z)) + (\v -> (yComp v, zComp v)) "eFieldPicIdealDipole.png" 20 + (eFieldIdealDipole kHat) + +type VectorLineIntegral = VectorField -> Curve -> Vec + +type CurveApprox = Curve -> [(Position,Vec)] + +vectorLineIntegral :: CurveApprox -> VectorField -> Curve -> Vec +vectorLineIntegral approx vF c + = sumV [vF r' ^* magnitude dl' | (r',dl') <- approx c] + +eFieldFromLineCharge + :: ScalarField -- linear charge density lambda + -> Curve -- geometry of the line charge + -> VectorField -- electric field (in V/m) +eFieldFromLineCharge lambda c r + = let k = 1 / (4 * pi * epsilon0) + integrand r' = lambda r' *^ d ^/ magnitude d ** 3 + where d = displacement r' r + in k *^ vectorLineIntegral (curveSample 1000) integrand c + +lineDipoleSodiumChloride :: ChargeDistribution +lineDipoleSodiumChloride = lineDipole (vec 0 0 2.99e-29) 2.36e-10 + +eFieldLineDipole :: VectorField +eFieldLineDipole = eField lineDipoleSodiumChloride + +type VectorSurfaceIntegral = VectorField -> Surface -> Vec + +type SurfaceApprox = Surface -> [(Position,Vec)] + +vectorSurfaceIntegral :: SurfaceApprox -> VectorField -> Surface -> Vec +vectorSurfaceIntegral approx vF s + = sumV [vF r' ^* magnitude da' | (r',da') <- approx s] + +eFieldFromSurfaceCharge + :: ScalarField -- surface charge density sigma + -> Surface -- geometry of the surface charge + -> VectorField -- electric field (in V/m) +eFieldFromSurfaceCharge sigma s r + = let k = 1 / (4 * pi * epsilon0) + integrand r' = sigma r' *^ d ^/ magnitude d ** 3 + where d = displacement r' r + in k *^ vectorSurfaceIntegral (surfaceSample 200) integrand s + +eFieldDiskCap :: VectorField +eFieldDiskCap = eField $ diskCap 0.05 0.04 2e-8 + +eFieldPicDiskCap :: IO () +eFieldPicDiskCap = vfGrad (**0.2) (\(x,z) -> cart (0.1*x) 0 (0.1*z)) + (\v -> (xComp v, zComp v)) "eFieldPicDiskCap.png" 20 + eFieldDiskCap + +type VectorVolumeIntegral = VectorField -> Volume -> Vec + +type VolumeApprox = Volume -> [(Position,R)] + +vectorVolumeIntegral :: VolumeApprox -> VectorField -> Volume -> Vec +vectorVolumeIntegral approx vF vol + = sumV [vF r' ^* dv' | (r',dv') <- approx vol] + +eFieldFromVolumeCharge + :: ScalarField -- volume charge density rho + -> Volume -- geometry of the volume charge + -> VectorField -- electric field (in V/m) +eFieldFromVolumeCharge rho v r + = let k = 1 / (4 * pi * epsilon0) + integrand r' = rho r' *^ d ^/ magnitude d ** 3 + where d = displacement r' r + in k *^ vectorVolumeIntegral (volumeSample 50) integrand v + +type ScalarLineIntegral = ScalarField -> Curve -> R + +scalarLineIntegral :: CurveApprox -> ScalarField -> Curve -> R +scalarLineIntegral approx f c + = sum [f r' * magnitude dl' | (r',dl') <- approx c] + +type ScalarSurfaceIntegral = ScalarField -> Surface -> R + +scalarSurfaceIntegral :: SurfaceApprox -> ScalarField -> Surface -> R +scalarSurfaceIntegral approx f s + = sum [f r' * magnitude da' | (r',da') <- approx s] + +type ScalarVolumeIntegral = ScalarField -> Volume -> R + +scalarVolumeIntegral :: VolumeApprox -> ScalarField -> Volume -> R +scalarVolumeIntegral approx f vol + = sum [f r' * dv' | (r',dv') <- approx vol] + +curveSample :: Int -> Curve -> [(Position,Vec)] +curveSample n c + = let segCent :: Segment -> Position + segCent (p1,p2) = shiftPosition ((rVF p1 ^+^ rVF p2) ^/ 2) origin + segDisp :: Segment -> Vec + segDisp = uncurry displacement + in [(segCent seg, segDisp seg) | seg <- segments n c] + +type Segment = (Position,Position) + +segments :: Int -> Curve -> [Segment] +segments n (Curve g a b) + = let ps = map g $ linSpaced n a b + in zip ps (tail ps) + +linSpaced :: Int -> R -> R -> [R] +linSpaced n x0 x1 = take (n+1) [x0, x0+dx .. x1] + where dx = (x1 - x0) / fromIntegral n + +surfaceSample :: Int -> Surface -> [(Position,Vec)] +surfaceSample n s = [(triCenter tri, triArea tri) | tri <- triangles n s] + +data Triangle = Tri Position Position Position + +triCenter :: Triangle -> Position +triCenter (Tri p1 p2 p3) + = shiftPosition ((rVF p1 ^+^ rVF p2 ^+^ rVF p3) ^/ 3) origin + +triArea :: Triangle -> Vec -- vector area +triArea (Tri p1 p2 p3) = 0.5 *^ (displacement p1 p2 >< displacement p2 p3) + +triangles :: Int -> Surface -> [Triangle] +triangles n (Surface g sl su tl tu) + = let sts = [[(s,t) | t <- linSpaced n (tl s) (tu s)] + | s <- linSpaced n sl su] + stSquares = [( sts !! j !! k + , sts !! (j+1) !! k + , sts !! (j+1) !! (k+1) + , sts !! j !! (k+1)) + | j <- [0..n-1], k <- [0..n-1]] + twoTriangles (pp1,pp2,pp3,pp4) + = [Tri (g pp1) (g pp2) (g pp3),Tri (g pp1) (g pp3) (g pp4)] + in concatMap twoTriangles stSquares + +volumeSample :: Int -> Volume -> [(Position,R)] +volumeSample n v = [(tetCenter tet, tetVolume tet) | tet <- tetrahedrons n v] + +data Tet = Tet Position Position Position Position + +tetCenter :: Tet -> Position +tetCenter (Tet p1 p2 p3 p4) + = shiftPosition ((rVF p1 ^+^ rVF p2 ^+^ rVF p3 ^+^ rVF p4) ^/ 4) origin + +tetVolume :: Tet -> R +tetVolume (Tet p1 p2 p3 p4) + = abs $ (d1 <.> (d2 >< d3)) / 6 + where + d1 = displacement p1 p4 + d2 = displacement p2 p4 + d3 = displacement p3 p4 + +data ParamCube + = PC { v000 :: (R,R,R) + , v001 :: (R,R,R) + , v010 :: (R,R,R) + , v011 :: (R,R,R) + , v100 :: (R,R,R) + , v101 :: (R,R,R) + , v110 :: (R,R,R) + , v111 :: (R,R,R) + } + +tetrahedrons :: Int -> Volume -> [Tet] +tetrahedrons n (Volume g sl su tl tu ul uu) + = let stus = [[[(s,t,u) | u <- linSpaced n (ul s t) (uu s t)] + | t <- linSpaced n (tl s) (tu s)] + | s <- linSpaced n sl su] + stCubes = [PC (stus !! j !! k !! l ) + (stus !! j !! k !! (l+1)) + (stus !! j !! (k+1) !! l ) + (stus !! j !! (k+1) !! (l+1)) + (stus !! (j+1) !! k !! l ) + (stus !! (j+1) !! k !! (l+1)) + (stus !! (j+1) !! (k+1) !! l ) + (stus !! (j+1) !! (k+1) !! (l+1)) + | j <- [0..n-1], k <- [0..n-1], l <- [0..n-1]] + tets (PC c000 c001 c010 c011 c100 c101 c110 c111) + = [Tet (g c000) (g c100) (g c010) (g c001) + ,Tet (g c011) (g c111) (g c001) (g c010) + ,Tet (g c110) (g c010) (g c100) (g c111) + ,Tet (g c101) (g c001) (g c111) (g c100) + ,Tet (g c111) (g c100) (g c010) (g c001) + ] + in concatMap tets stCubes + +type Field a = Position -> a + +class AbstractVector a where + zeroVector :: a + add :: a -> a -> a + scale :: R -> a -> a + +sumG :: AbstractVector a => [a] -> a +sumG = foldr add zeroVector + +generalLineIntegral + :: AbstractVector a => CurveApprox -> Field a -> Curve -> a +generalLineIntegral approx f c + = sumG [scale (magnitude dl') (f r') | (r',dl') <- approx c] + +dottedSurfaceIntegral :: SurfaceApprox -> VectorField -> Surface -> R +dottedSurfaceIntegral approx vF s + = sum [vF r' <.> da' | (r',da') <- approx s] + +electricFluxFromField :: VectorField -> Surface -> R +electricFluxFromField = undefined + +electricFluxFromCharge :: ChargeDistribution -> Surface -> R +electricFluxFromCharge dist = undefined dist + +eFieldFromSurfaceChargeP :: SurfaceApprox -> ScalarField -> Surface + -> VectorField +eFieldFromSurfaceChargeP approx sigma s r + = sumV [eFieldFromPointCharge (sigma r' * magnitude da') r' r + | (r',da') <- approx s] + +surfaceArea :: Surface -> R +surfaceArea = undefined + +dottedLineIntegral :: CurveApprox -> VectorField -> Curve -> R +dottedLineIntegral approx f c = sum [f r' <.> dl' | (r',dl') <- approx c] + +electricPotentialFromField :: VectorField -- electric field + -> ScalarField -- electric potential +electricPotentialFromField ef r = undefined ef r diff --git a/src/LPFP/Electricity.hs b/src/LPFP/Electricity.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/Electricity.hs @@ -0,0 +1,99 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFP.Electricity +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 21 of the book Learn Physics with Functional Programming +-} + +module LPFP.Electricity where + +import LPFP.SimpleVec + ( Vec(..), R, (*^), iHat ) +import LPFP.Mechanics3D + ( ParticleState(..), defaultParticleState ) +import LPFP.MultipleObjects + ( TwoBodyForce, MultiParticleState(..), Force(..), statesMPS + , eulerCromerMPS, centralForce ) +import Graphics.Gnuplot.Simple + ( Attribute(..), plotPaths ) + +type Charge = R + +elementaryCharge :: Charge +elementaryCharge = 1.602176634e-19 -- in Coulombs + +coulombMagnitude :: Charge -> Charge -> R -> R +coulombMagnitude q1 q2 r + = let k = 9e9 -- in N m^2 / C^2 + in k * abs (q1 * q2) / r**2 + +coulombForce :: TwoBodyForce +coulombForce st1 st2 + = let k = 9e9 -- N m^2 / C^2 + q1 = charge st1 + q2 = charge st2 + in centralForce (\r -> k * q1 * q2 / r**2) st1 st2 + +twoProtonStates :: R -- time step + -> MultiParticleState -- initial 2-particle state + -> [MultiParticleState] -- infinite list of states +twoProtonStates dt + = statesMPS (eulerCromerMPS dt) [InternalForce 1 0 coulombForce] + +-- protons are released from rest +initialTwoProtonState :: R -- initial separation + -> MultiParticleState +initialTwoProtonState d + = let protonMass = 1.673e-27 -- in kg + in MPS [defaultParticleState { mass = protonMass + , charge = elementaryCharge + , posVec = (-d/2) *^ iHat + } + ,defaultParticleState { mass = protonMass + , charge = elementaryCharge + , posVec = ( d/2) *^ iHat + } + ] + +oneProtonVelocity :: R -- dt + -> R -- starting separation + -> [(R,R)] -- (time,velocity) pairs +oneProtonVelocity dt d + = let state0 = initialTwoProtonState d + in [(time st2, xComp $ velocity st2) + | MPS [_,st2] <- twoProtonStates dt state0] + +tvPairs :: [(R,R)] +tvPairs = takeWhile (\(t,_) -> t <= 2e-2) $ + oneProtonVelocity 1e-5 1e-2 + +velocityPlot :: IO () +velocityPlot + = plotPaths [Title "Two protons released from 1 cm" + ,XLabel "Time (s)" + ,YLabel "Proton velocity (m/s)" + ,PNG "protons.png" + ,Key Nothing + ] $ [tvPairs + ,[(t,1379*t) | t <- [0,1e-5..4e-3]] + ,[(t,3.71) | t <- [0,1e-3..2e-2]]] + +oneProtonPosition :: R -- dt + -> R -- starting separation + -> [(R,R)] -- (time,position) pairs +oneProtonPosition dt d + = undefined dt d + +positionPlot :: IO () +positionPlot = plotPaths [Title "Two protons released from 1 cm" + ,XLabel "Time (s)" + ,YLabel "Proton position (m)" + ,PNG "ProtonPosition.png" + ,Key Nothing + ] $ [undefined $ oneProtonPosition 1e-5 1e-2 + ,undefined :: [(R,R)]] diff --git a/src/LPFP/Geometry.hs b/src/LPFP/Geometry.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/Geometry.hs @@ -0,0 +1,134 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFP.Geometry +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 23 of the book Learn Physics with Functional Programming +-} + +module LPFP.Geometry where + +import LPFP.SimpleVec ( R, Vec, (*^) ) +import LPFP.CoordinateSystems ( Position, cylindrical, spherical, cart, cyl, sph + , shiftPosition, displacement ) + +data Curve = Curve { curveFunc :: R -> Position + , startingCurveParam :: R -- t_a + , endingCurveParam :: R -- t_b + } + +circle2 :: Curve +circle2 = Curve (\t -> cart (2 * cos t) (2 * sin t) 0) 0 (2*pi) + +circle2' :: Curve +circle2' = Curve (\phi -> cyl 2 phi 0) 0 (2*pi) + +unitCircle :: Curve +unitCircle = Curve (\t -> cyl 1 t 0) 0 (2 * pi) + +straightLine :: Position -- starting position + -> Position -- ending position + -> Curve -- straight-line curve +straightLine r1 r2 = let d = displacement r1 r2 + f t = shiftPosition (t *^ d) r1 + in Curve f 0 1 + +data Surface = Surface { surfaceFunc :: (R,R) -> Position + , lowerLimit :: R -- s_l + , upperLimit :: R -- s_u + , lowerCurve :: R -> R -- t_l(s) + , upperCurve :: R -> R -- t_u(s) + } + +unitSphere :: Surface +unitSphere = Surface (\(th,phi) -> cart (sin th * cos phi) + (sin th * sin phi) + (cos th)) + 0 pi (const 0) (const $ 2*pi) + +unitSphere' :: Surface +unitSphere' = Surface (\(th,phi) -> sph 1 th phi) + 0 pi (const 0) (const $ 2*pi) + +parabolaSurface :: Surface +parabolaSurface = Surface (\(x,y) -> cart x y 0) + (-2) 2 (\x -> x*x) (const 4) + +shiftSurface :: Vec -> Surface -> Surface +shiftSurface d (Surface g sl su tl tu) + = Surface (shiftPosition d . g) sl su tl tu + +centeredSphere :: R -> Surface +centeredSphere r = Surface (\(th,phi) -> sph r th phi) + 0 pi (const 0) (const $ 2*pi) + +sphere :: R -> Position -> Surface +sphere radius center + = shiftSurface (displacement (cart 0 0 0) center) + (centeredSphere radius) + +northernHemisphere :: Surface +northernHemisphere = Surface (\(th,phi) -> sph 1 th phi) + 0 (pi/2) (const 0) (const $ 2*pi) + +disk :: R -> Surface +disk radius = Surface (\(s,phi) -> cyl s phi 0) + 0 radius (const 0) (const (2*pi)) + +unitCone :: R -> Surface +unitCone theta = Surface (\(r,phi) -> sph r theta phi) + 0 1 (const 0) (const (2*pi)) + +data Volume = Volume { volumeFunc :: (R,R,R) -> Position + , loLimit :: R -- s_l + , upLimit :: R -- s_u + , loCurve :: R -> R -- t_l(s) + , upCurve :: R -> R -- t_u(s) + , loSurf :: R -> R -> R -- u_l(s,t) + , upSurf :: R -> R -> R -- u_u(s,t) + } + +unitBall :: Volume +unitBall = Volume spherical 0 1 (const 0) (const pi) + (\_ _ -> 0) (\_ _ -> 2*pi) + +centeredCylinder :: R -- radius + -> R -- height + -> Volume -- cylinder +centeredCylinder radius height + = Volume cylindrical 0 radius (const 0) (const (2*pi)) + (\_ _ -> 0) (\_ _ -> height) + +circle :: Position -- center position + -> R -- radius + -> Curve +circle r radius = undefined r radius + +square :: Curve +square = Curve squareFunc 0 4 + +squareFunc :: R -> Position +squareFunc t + | t < 1 = cart undefined (-1) 0 + | 1 <= t && t < 2 = cart 1 undefined 0 + | 2 <= t && t < 3 = cart undefined 1 0 + | otherwise = cart (-1) undefined 0 + +northernHalfBall :: Volume +northernHalfBall = undefined + +centeredBall :: R -> Volume +centeredBall = undefined + +shiftVolume :: Vec -> Volume -> Volume +shiftVolume = undefined + +quarterDiskBoundary :: R -> Curve +quarterDiskBoundary = undefined + +quarterCylinder :: R -> R -> Volume +quarterCylinder = undefined diff --git a/src/LPFP/Integrals.hs b/src/LPFP/Integrals.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/Integrals.hs @@ -0,0 +1,143 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFP.Integrals +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code needed in chapter 24 of the book Learn Physics with Functional Programming +-} + +module LPFP.Integrals where + +import LPFP.SimpleVec + ( R, Vec, (^+^), (*^), (^*), (^/), (<.>), (><), sumV, magnitude ) +import LPFP.CoordinateSystems ( Position, ScalarField, VectorField + , displacement, shiftPosition, origin, rVF ) +import LPFP.Geometry ( Curve(..), Surface(..), Volume(..) ) + +type CurveApprox = Curve -> [(Position,Vec)] + +type SurfaceApprox = Surface -> [(Position,Vec)] + +type VolumeApprox = Volume -> [(Position,R)] + +scalarLineIntegral :: CurveApprox -> ScalarField -> Curve -> R +scalarLineIntegral approx f c + = sum [f r' * magnitude dl' | (r',dl') <- approx c] + +scalarSurfaceIntegral :: SurfaceApprox -> ScalarField -> Surface -> R +scalarSurfaceIntegral approx f s + = sum [f r' * magnitude da' | (r',da') <- approx s] + +scalarVolumeIntegral :: VolumeApprox -> ScalarField -> Volume -> R +scalarVolumeIntegral approx f vol + = sum [f r' * dv' | (r',dv') <- approx vol] + +vectorLineIntegral :: CurveApprox -> VectorField -> Curve -> Vec +vectorLineIntegral approx vF c + = sumV [vF r' ^* magnitude dl' | (r',dl') <- approx c] + +vectorSurfaceIntegral :: SurfaceApprox -> VectorField -> Surface -> Vec +vectorSurfaceIntegral approx vF s + = sumV [vF r' ^* magnitude da' | (r',da') <- approx s] + +vectorVolumeIntegral :: VolumeApprox -> VectorField -> Volume -> Vec +vectorVolumeIntegral approx vF vol + = sumV [vF r' ^* dv' | (r',dv') <- approx vol] + +curveSample :: Int -> Curve -> [(Position,Vec)] +curveSample n c + = let segCent :: Segment -> Position + segCent (p1,p2) = shiftPosition ((rVF p1 ^+^ rVF p2) ^/ 2) origin + segDisp :: Segment -> Vec + segDisp = uncurry displacement + in [(segCent seg, segDisp seg) | seg <- segments n c] + +type Segment = (Position,Position) +segments :: Int -> Curve -> [Segment] +segments n (Curve g a b) + = let ps = map g $ linSpaced n a b + in zip ps (tail ps) + +linSpaced :: Int -> R -> R -> [R] +linSpaced n x0 x1 = take (n+1) [x0, x0+dx .. x1] + where dx = (x1 - x0) / fromIntegral n + +surfaceSample :: Int -> Surface -> [(Position,Vec)] +surfaceSample n s = [(triCenter tri, triArea tri) | tri <- triangles n s] + +data Triangle = Tri Position Position Position + +triCenter :: Triangle -> Position +triCenter (Tri p1 p2 p3) + = shiftPosition ((rVF p1 ^+^ rVF p2 ^+^ rVF p3) ^/ 3) origin + +triArea :: Triangle -> Vec -- vector area +triArea (Tri p1 p2 p3) = 0.5 *^ (displacement p1 p2 >< displacement p2 p3) + +triangles :: Int -> Surface -> [Triangle] +triangles n (Surface g sl su tl tu) + = let sts = [[(s,t) | t <- linSpaced n (tl s) (tu s)] + | s <- linSpaced n sl su] + stSquares = [( sts !! j !! k + , sts !! (j+1) !! k + , sts !! (j+1) !! (k+1) + , sts !! j !! (k+1)) + | j <- [0..n-1], k <- [0..n-1]] + twoTriangles (pp1,pp2,pp3,pp4) + = [Tri (g pp1) (g pp2) (g pp3),Tri (g pp1) (g pp3) (g pp4)] + in concatMap twoTriangles stSquares + +volumeSample :: Int -> Volume -> [(Position,R)] +volumeSample n v = [(tetCenter tet, tetVolume tet) | tet <- tetrahedrons n v] + +data Tet = Tet Position Position Position Position + +tetCenter :: Tet -> Position +tetCenter (Tet p1 p2 p3 p4) + = shiftPosition ((rVF p1 ^+^ rVF p2 ^+^ rVF p3 ^+^ rVF p4) ^/ 4) origin + +tetVolume :: Tet -> R +tetVolume (Tet p1 p2 p3 p4) + = abs $ (d1 <.> (d2 >< d3)) / 6 + where + d1 = displacement p1 p4 + d2 = displacement p2 p4 + d3 = displacement p3 p4 + +data ParamCube + = PC { v000 :: (R,R,R) + , v001 :: (R,R,R) + , v010 :: (R,R,R) + , v011 :: (R,R,R) + , v100 :: (R,R,R) + , v101 :: (R,R,R) + , v110 :: (R,R,R) + , v111 :: (R,R,R) + } + +tetrahedrons :: Int -> Volume -> [Tet] +tetrahedrons n (Volume g sl su tl tu ul uu) + = let stus = [[[(s,t,u) | u <- linSpaced n (ul s t) (uu s t)] + | t <- linSpaced n (tl s) (tu s)] + | s <- linSpaced n sl su] + stCubes = [PC (stus !! j !! k !! l ) + (stus !! j !! k !! (l+1)) + (stus !! j !! (k+1) !! l ) + (stus !! j !! (k+1) !! (l+1)) + (stus !! (j+1) !! k !! l ) + (stus !! (j+1) !! k !! (l+1)) + (stus !! (j+1) !! (k+1) !! l ) + (stus !! (j+1) !! (k+1) !! (l+1)) + | j <- [0..n-1], k <- [0..n-1], l <- [0..n-1]] + tets (PC c000 c001 c010 c011 c100 c101 c110 c111) + = [Tet (g c000) (g c100) (g c010) (g c001) + ,Tet (g c011) (g c111) (g c001) (g c010) + ,Tet (g c110) (g c010) (g c100) (g c111) + ,Tet (g c101) (g c001) (g c111) (g c100) + ,Tet (g c111) (g c100) (g c010) (g c001) + ] + in concatMap tets stCubes diff --git a/src/LPFP/Lorentz.hs b/src/LPFP/Lorentz.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/Lorentz.hs @@ -0,0 +1,147 @@ +{-# OPTIONS -Wall #-} +{-# LANGUAGE MultiParamTypeClasses #-} + +{- | +Module : LPFP.Lorentz +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 28 of the book Learn Physics with Functional Programming +-} + +module LPFP.Lorentz where + +import LPFP.SimpleVec ( R, Vec, (^+^), (*^), (^*), (^/), (><), zeroV, magnitude ) +import LPFP.Mechanics1D ( RealVectorSpace(..), Diff(..), rungeKutta4 ) +import LPFP.Mechanics3D ( HasTime(..), simulateVis ) +import LPFP.CoordinateSystems ( Position(..), VectorField, cart, v3FromPos, origin + , shiftPosition, addVectorFields, visVec ) +import qualified Vis as V + +data ParticleFieldState = ParticleFieldState { mass :: R + , charge :: R + , time :: R + , position :: Position + , velocity :: Vec + , electricField :: VectorField + , magneticField :: VectorField } + +data DParticleFieldState = DParticleFieldState { dmdt :: R + , dqdt :: R + , dtdt :: R + , drdt :: Vec + , dvdt :: Vec + , dEdt :: VectorField + , dBdt :: VectorField } + +instance RealVectorSpace DParticleFieldState where + dst1 +++ dst2 + = DParticleFieldState { dmdt = dmdt dst1 + dmdt dst2 + , dqdt = dqdt dst1 + dqdt dst2 + , dtdt = dtdt dst1 + dtdt dst2 + , drdt = drdt dst1 ^+^ drdt dst2 + , dvdt = dvdt dst1 ^+^ dvdt dst2 + , dEdt = addVectorFields [dEdt dst1, dEdt dst2] + , dBdt = addVectorFields [dBdt dst1, dBdt dst2] + } + scale w dst + = DParticleFieldState { dmdt = w * dmdt dst + , dqdt = w * dqdt dst + , dtdt = w * dtdt dst + , drdt = w *^ drdt dst + , dvdt = w *^ dvdt dst + , dEdt = (w *^) . (dEdt dst) + , dBdt = (w *^) . (dBdt dst) + } + +instance Diff ParticleFieldState DParticleFieldState where + shift dt dst st + = ParticleFieldState + { mass = mass st + dmdt dst * dt + , charge = charge st + dqdt dst * dt + , time = time st + dtdt dst * dt + , position = shiftPosition (drdt dst ^* dt) (position st) + , velocity = velocity st ^+^ dvdt dst ^* dt + , electricField = \r -> electricField st r ^+^ dEdt dst r ^* dt + , magneticField = \r -> magneticField st r ^+^ dBdt dst r ^* dt + } + +instance HasTime ParticleFieldState where + timeOf = time + +lorentzForce :: ParticleFieldState -> Vec +lorentzForce (ParticleFieldState _m q _t r v eF bF) + = q *^ (eF r ^+^ v >< bF r) + +newtonSecondPFS :: ParticleFieldState -> DParticleFieldState +newtonSecondPFS st + = let v = velocity st + a = lorentzForce st ^/ mass st + in DParticleFieldState { dmdt = 0 -- dm/dt + , dqdt = 0 -- dq/dt + , dtdt = 1 -- dt/dt + , drdt = v -- dr/dt + , dvdt = a -- dv/dt + , dEdt = const zeroV -- dE/dt + , dBdt = const zeroV -- dB/dt + } + +pfsUpdate :: R -- time step + -> ParticleFieldState -> ParticleFieldState +pfsUpdate dt = rungeKutta4 dt newtonSecondPFS + +defaultPFS :: ParticleFieldState +defaultPFS = ParticleFieldState { mass = 0 + , charge = 0 + , time = 0 + , position = origin + , velocity = zeroV + , electricField = const zeroV + , magneticField = const zeroV } + +pfsVisObject :: R -- cube width + -> ParticleFieldState -> V.VisObject R +pfsVisObject width st + = let r = position st + xs = [-width/2, width/2] + es :: [(Position,Vec)] + es = [(cart x y z, electricField st (cart x y z)) + | x <- xs, y <- xs, z <- xs] + maxE = maximum $ map (magnitude . snd) es + bs :: [(Position,Vec)] + bs = [(cart x y z, magneticField st (cart x y z)) + | x <- xs, y <- xs, z <- xs] + maxB = maximum $ map (magnitude . snd) bs + metersPerVis = width/2 + in V.VisObjects [ vectorsVisObject metersPerVis (2*maxE) es V.blue + , vectorsVisObject metersPerVis (2*maxB) bs V.red + , V.Trans (v3FromPos (scalePos metersPerVis r)) + (V.Sphere 0.1 V.Solid V.green) + ] + +vectorsVisObject :: R -- scale factor, meters per Vis unit + -> R -- scale factor, vector field units per Vis unit + -> [(Position,Vec)] -- positions to show the field + -> V.Color + -> V.VisObject R +vectorsVisObject metersPerVis unitsPerVis pvs color + = V.VisObjects [V.Trans (v3FromPos (scalePos metersPerVis r)) $ + visVec color (v ^/ unitsPerVis) | (r,v) <- pvs] + +scalePos :: R -> Position -> Position +scalePos metersPerVis (Cart x y z) + = Cart (x/metersPerVis) (y/metersPerVis) (z/metersPerVis) + +animatePFS :: R -- time scale factor + -> Int -- animation rate + -> R -- display width + -> ParticleFieldState -- initial state + -> IO () +animatePFS tsf ar width st + = simulateVis tsf ar st (pfsVisObject width) pfsUpdate + +newtonSecondPFS' :: [ParticleFieldState -> Vec] + -> ParticleFieldState -> DParticleFieldState +newtonSecondPFS' fs st = undefined fs st diff --git a/src/LPFP/MOExamples.hs b/src/LPFP/MOExamples.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/MOExamples.hs @@ -0,0 +1,353 @@ +{-# 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) diff --git a/src/LPFP/MagneticField.hs b/src/LPFP/MagneticField.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/MagneticField.hs @@ -0,0 +1,117 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFP.MagneticField +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 27 of the book Learn Physics with Functional Programming +-} + +module LPFP.MagneticField where + +import LPFP.SimpleVec ( Vec(..), R + , (^-^), (*^), (^/), (<.>), (><) + , magnitude, kHat, zComp ) +import LPFP.CoordinateSystems + ( VectorField + , rVF, displacement, addVectorFields, cart, vfGrad ) +import LPFP.Geometry ( Curve(..), Surface(..), Volume(..) ) +import LPFP.ElectricField + ( curveSample, surfaceSample, volumeSample + , vectorSurfaceIntegral, vectorVolumeIntegral, mu0 ) +import LPFP.Current + ( Current, CurrentDistribution(..) + , wireSolenoid, wireToroid, crossedLineIntegral, circularCurrentLoop ) + +bFieldFromLineCurrent + :: Current -- current (in Amps) + -> Curve + -> VectorField -- magnetic field (in Tesla) +bFieldFromLineCurrent i c r + = let coeff = -mu0 * i / (4 * pi) -- SI units + integrand r' = d ^/ magnitude d ** 3 + where d = displacement r' r + in coeff *^ crossedLineIntegral (curveSample 1000) integrand c + +bField :: CurrentDistribution -> VectorField +bField (LineCurrent i c) = bFieldFromLineCurrent i c +bField (SurfaceCurrent kC s) = bFieldFromSurfaceCurrent kC s +bField (VolumeCurrent j v) = bFieldFromVolumeCurrent j v +bField (MultipleCurrents cds) = addVectorFields $ map bField cds + +circleB :: VectorField -- magnetic field +circleB = bField $ circularCurrentLoop 0.25 10 + +bFieldPicLoop :: IO () +bFieldPicLoop + = vfGrad (**0.2) (\(y,z) -> cart 0 y z) (\v -> (yComp v, zComp v)) + "bFieldPicLoop.png" 20 circleB + +bFieldIdealDipole :: Vec -- magnetic dipole moment + -> VectorField -- magnetic field +bFieldIdealDipole m r + = let coeff = mu0 / (4 * pi) -- SI units + rMag = magnitude (rVF r) + rUnit = rVF r ^/ rMag + in coeff *^ (1 / rMag**3) *^ (3 *^ (m <.> rUnit) *^ rUnit ^-^ m) + +bFieldPicIdealDipole :: IO () +bFieldPicIdealDipole + = vfGrad (**0.2) (\(y,z) -> cart 0 y z) (\v -> (yComp v, zComp v)) + "bFieldPicIdealDipole.png" 20 (bFieldIdealDipole kHat) + +bFieldPicSolenoid10 :: IO () +bFieldPicSolenoid10 = vfGrad (**0.2) (\(y,z) -> cart 0 (0.02*y) (0.02*z)) + (\v -> (yComp v, zComp v)) "bFieldPicSolenoid10.png" 20 + (bField $ wireSolenoid 0.01 0.1 100 10) + +bFieldPicSolenoid100 :: IO () +bFieldPicSolenoid100 = vfGrad (**0.2) (\(y,z) -> cart 0 (0.02*y) (0.02*z)) + (\v -> (yComp v, zComp v)) "bFieldPicSolenoid100.png" 20 + (bField $ wireSolenoid 0.01 0.1 1000 10) + +bFieldWireToroid :: VectorField +bFieldWireToroid = bField (wireToroid 0.3 1 50 10) + +bFieldPicWireToroid :: IO () +bFieldPicWireToroid + = vfGrad (**0.2) (\(x,y) -> cart (1.5*x) (1.5*y) 0) + (\v -> (xComp v, yComp v)) "bFieldPicWireToroid.png" 20 bFieldWireToroid + +bFieldFromSurfaceCurrent + :: VectorField -- surface current density + -> Surface -- surface across which current flows + -> VectorField -- magnetic field (in T) +bFieldFromSurfaceCurrent kCurrent s r + = let coeff = mu0 / (4 * pi) -- SI units + integrand r' = (kCurrent r' >< d) ^/ magnitude d ** 3 + where d = displacement r' r + in coeff *^ vectorSurfaceIntegral (surfaceSample 200) integrand s + +bFieldFromVolumeCurrent + :: VectorField -- volume current density + -> Volume -- volume throughout which current flows + -> VectorField -- magnetic field (in T) +bFieldFromVolumeCurrent j vol r + = let coeff = mu0 / (4 * pi) -- SI units + integrand r' = (j r' >< d) ^/ magnitude d ** 3 + where d = displacement r' r + in coeff *^ vectorVolumeIntegral (volumeSample 50) integrand vol + +magneticFluxFromField :: VectorField -> Surface -> R +magneticFluxFromField = undefined + +magneticFluxFromCurrent :: CurrentDistribution -> Surface -> R +magneticFluxFromCurrent = undefined + +visLoop :: IO () +visLoop = undefined + +bFieldPicSolenoid1000 :: IO () +bFieldPicSolenoid1000 + = vfGrad (**0.2) (\(y,z) -> cart 0 (0.02*y) (0.02*z)) + (\v -> (yComp v, zComp v)) "bFieldPicSolenoid1000.png" 20 + (bField $ wireSolenoid 0.01 0.1 10000 10) diff --git a/src/LPFP/Maxwell.hs b/src/LPFP/Maxwell.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/Maxwell.hs @@ -0,0 +1,222 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFP.Maxwell +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 29 of the book Learn Physics with Functional Programming +-} + +module LPFP.Maxwell where + +import LPFP.SimpleVec + ( R, Vec(..), (^/), (^+^), (^-^), (*^) + , vec, negateV, magnitude, xComp, yComp, zComp, iHat, jHat, kHat ) +import LPFP.CoordinateSystems + ( ScalarField, VectorField + , cart, shiftPosition, rVF, magRad ) +import LPFP.ElectricField ( cSI, mu0 ) +import qualified Data.Map.Strict as M +import qualified Diagrams.Prelude as D +import Diagrams.Prelude + ( Diagram, Colour + , PolyType(..), PolyOrientation(..), PolygonOpts(..), V2(..) + , (#), rotate, deg, rad, polygon, sinA, dims, p2 + , fc, none, lw, blend ) +import Diagrams.Backend.Cairo ( B, renderCairo ) + +directionalDerivative :: Vec -> ScalarField -> ScalarField +directionalDerivative d f r + = (f (shiftPosition (d ^/ 2) r) - f (shiftPosition (negateV d ^/ 2) r)) + / magnitude d + +curl :: R -> VectorField -> VectorField +curl a vf r + = let vx = xComp . vf + vy = yComp . vf + vz = zComp . vf + derivX = directionalDerivative (a *^ iHat) + derivY = directionalDerivative (a *^ jHat) + derivZ = directionalDerivative (a *^ kHat) + in (derivY vz r - derivZ vy r) *^ iHat + ^+^ (derivZ vx r - derivX vz r) *^ jHat + ^+^ (derivX vy r - derivY vx r) *^ kHat + +type FieldState = (R -- time t + ,VectorField -- electric field E + ,VectorField -- magnetic field B + ) + +maxwellUpdate :: R -- dx + -> R -- dt + -> (R -> VectorField) -- J + -> FieldState -> FieldState +maxwellUpdate dx dt j (t,eF,bF) + = let t' = t + dt + eF' r = eF r ^+^ cSI**2 *^ dt *^ (curl dx bF r ^-^ mu0 *^ j t r) + bF' r = bF r ^-^ dt *^ curl dx eF r + in (t',eF',bF') + +maxwellEvolve :: R -- dx + -> R -- dt + -> (R -> VectorField) -- J + -> FieldState -> [FieldState] +maxwellEvolve dx dt j st0 = iterate (maxwellUpdate dx dt j) st0 + +exLocs, eyLocs, ezLocs, bxLocs, byLocs, bzLocs :: [(Int,Int,Int)] +exLocs = [(nx,ny,nz) | nx <- odds , ny <- evens, nz <- evens] +eyLocs = [(nx,ny,nz) | nx <- evens, ny <- odds , nz <- evens] +ezLocs = [(nx,ny,nz) | nx <- evens, ny <- evens, nz <- odds ] +bxLocs = [(nx,ny,nz) | nx <- evens, ny <- odds , nz <- odds ] +byLocs = [(nx,ny,nz) | nx <- odds , ny <- evens, nz <- odds ] +bzLocs = [(nx,ny,nz) | nx <- odds , ny <- odds , nz <- evens] + +spaceStepsCE :: Int +spaceStepsCE = 40 + +hiEven :: Int +hiEven = 2 * spaceStepsCE + +evens :: [Int] +evens = [-hiEven, -hiEven + 2 .. hiEven] + +odds :: [Int] +odds = [-hiEven + 1, -hiEven + 3 .. hiEven - 1] + +data StateFDTD = StateFDTD {timeFDTD :: R + ,stepX :: R + ,stepY :: R + ,stepZ :: R + ,eField :: M.Map (Int,Int,Int) R + ,bField :: M.Map (Int,Int,Int) R + } deriving Show + +initialStateFDTD :: R -> StateFDTD +initialStateFDTD spatialStep + = StateFDTD {timeFDTD = 0 + ,stepX = spatialStep + ,stepY = spatialStep + ,stepZ = spatialStep + ,eField = M.fromList [(loc,0) | loc <- exLocs++eyLocs++ezLocs] + ,bField = M.fromList [(loc,0) | loc <- bxLocs++byLocs++bzLocs] + } + +lookupAZ :: Ord k => k -> M.Map k R -> R +lookupAZ key m = case M.lookup key m of + Nothing -> 0 + Just x -> x + +partialX,partialY,partialZ :: R -> M.Map (Int,Int,Int) R -> (Int,Int,Int) -> R +partialX dx m (i,j,k) = (lookupAZ (i+1,j,k) m - lookupAZ (i-1,j,k) m) / dx +partialY dy m (i,j,k) = (lookupAZ (i,j+1,k) m - lookupAZ (i,j-1,k) m) / dy +partialZ dz m (i,j,k) = (lookupAZ (i,j,k+1) m - lookupAZ (i,j,k-1) m) / dz + +curlEx,curlEy,curlEz,curlBx,curlBy,curlBz :: StateFDTD -> (Int,Int,Int) -> R +curlBx (StateFDTD _ _ dy dz _ b) loc = partialY dy b loc - partialZ dz b loc +curlBy (StateFDTD _ dx _ dz _ b) loc = partialZ dz b loc - partialX dx b loc +curlBz (StateFDTD _ dx dy _ _ b) loc = partialX dx b loc - partialY dy b loc +curlEx (StateFDTD _ _ dy dz e _) loc = partialY dy e loc - partialZ dz e loc +curlEy (StateFDTD _ dx _ dz e _) loc = partialZ dz e loc - partialX dx e loc +curlEz (StateFDTD _ dx dy _ e _) loc = partialX dx e loc - partialY dy e loc + +stateUpdate :: R -- dt + -> (R -> VectorField) -- current density J + -> StateFDTD -> StateFDTD +stateUpdate dt j st0@(StateFDTD t _dx _dy _dz _e _b) + = let st1 = updateE dt (j t) st0 + st2 = updateB dt st1 + in st2 + +updateE :: R -- time step dt + -> VectorField -- current density J + -> StateFDTD -> StateFDTD +updateE dt jVF st + = st { timeFDTD = timeFDTD st + dt / 2 + , eField = M.mapWithKey (updateEOneLoc dt jVF st) (eField st) } + +updateB :: R -> StateFDTD -> StateFDTD +updateB dt st + = st { timeFDTD = timeFDTD st + dt / 2 + , bField = M.mapWithKey (updateBOneLoc dt st) (bField st) } + +updateEOneLoc :: R -> VectorField -> StateFDTD -> (Int,Int,Int) -> R -> R +updateEOneLoc dt jVF st (nx,ny,nz) ec + = let r = cart (fromIntegral nx * stepX st / 2) + (fromIntegral ny * stepY st / 2) + (fromIntegral nz * stepZ st / 2) + Vec jx jy jz = jVF r + in case (odd nx, odd ny, odd nz) of + (True , False, False) + -> ec + cSI**2 * (curlBx st (nx,ny,nz) - mu0 * jx) * dt -- Ex + (False, True , False) + -> ec + cSI**2 * (curlBy st (nx,ny,nz) - mu0 * jy) * dt -- Ey + (False, False, True ) + -> ec + cSI**2 * (curlBz st (nx,ny,nz) - mu0 * jz) * dt -- Ez + _ -> error "updateEOneLoc passed bad indices" + +updateBOneLoc :: R -> StateFDTD -> (Int,Int,Int) -> R -> R +updateBOneLoc dt st (nx,ny,nz) bc + = case (odd nx, odd ny, odd nz) of + (False, True , True ) -> bc - curlEx st (nx,ny,nz) * dt -- Bx + (True , False, True ) -> bc - curlEy st (nx,ny,nz) * dt -- By + (True , True , False) -> bc - curlEz st (nx,ny,nz) * dt -- Bz + _ -> error "updateBOneLoc passed bad indices" + +jGaussian :: R -> VectorField +jGaussian t r + = let wavelength = 1.08 -- meters + frequency = cSI / wavelength -- Hz + j0 = 77.5 -- A/m^2 + l = 0.108 -- meters + rMag = magnitude (rVF r) -- meters + in j0 *^ exp (-rMag**2 / l**2) *^ cos (2*pi*frequency*t) *^ kHat + +makeEpng :: (Colour R, Colour R) -> (Int,StateFDTD) -> IO () +makeEpng (scol,zcol) (n,StateFDTD _ _ _ _ em _) + = let threeDigitString = reverse $ take 3 $ reverse ("00" ++ show n) + pngFilePath = "MaxVF" ++ threeDigitString ++ ".png" + strongE = 176 -- V/m + vs = [((fromIntegral nx, fromIntegral nz),(xComp ev, zComp ev)) + | nx <- evens, nz <- evens, abs nx <= 50, abs nz <= 50 + , let ev = getAverage (nx,0,nz) em ^/ strongE] + in gradientVectorPNG pngFilePath (scol,zcol) vs + +getAverage :: (Int,Int,Int) -- (even,even,even) or (odd,odd,odd) + -> M.Map (Int,Int,Int) R + -> Vec +getAverage (i,j,k) m + = let vXl = lookupAZ (i-1,j ,k ) m + vYl = lookupAZ (i ,j-1,k ) m + vZl = lookupAZ (i ,j ,k-1) m + vXr = lookupAZ (i+1,j ,k ) m + vYr = lookupAZ (i ,j+1,k ) m + vZr = lookupAZ (i ,j ,k+1) m + in vec ((vXl+vXr)/2) ((vYl+vYr)/2) ((vZl+vZr)/2) + +gradientVectorPNG :: FilePath + -> (Colour R, Colour R) + -> [((R,R),(R,R))] + -> IO () +gradientVectorPNG fileName (scol,zcol) vs + = let maxX = maximum $ map fst $ map fst $ vs + normalize (x,y) = (x/maxX,y/maxX) + array = [(normalize (x,y), magRad v) | ((x,y),v) <- vs] + arrowMagRadColors :: R -- magnitude + -> R -- angle in radians, ccw from x axis + -> Diagram B + arrowMagRadColors mag th + = let r = sinA (15 D.@@ deg) / sinA (60 D.@@ deg) + myType = PolyPolar [120 D.@@ deg, 0 D.@@ deg, 45 D.@@ deg + , 30 D.@@ deg, 45 D.@@ deg, 0 D.@@ deg + ,120 D.@@ deg] + [1,1,r,1,1,r,1,1] + myOpts = PolygonOpts myType NoOrient (p2 (0,0)) + in D.scale 0.5 $ polygon myOpts # lw none # + fc (blend mag scol zcol) # rotate (th D.@@ rad) + step = 2 / (sqrt $ fromIntegral $ length vs) + scaledArrow m th = D.scale step $ arrowMagRadColors m th + pic = D.position [(p2 pt, scaledArrow m th) | (pt,(m,th)) <- array] + in renderCairo fileName (dims (V2 1024 1024)) pic diff --git a/src/LPFP/Mechanics1D.hs b/src/LPFP/Mechanics1D.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/Mechanics1D.hs @@ -0,0 +1,273 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-} + +{- | +Module : LPFP.Mechanics1D +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 15 of the book Learn Physics with Functional Programming +-} + +module LPFP.Mechanics1D where + +import Graphics.Gnuplot.Simple + +import LPFP.Newton2 ( fAir ) + +import LPFP.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 m fs (t,x0,v0) + = let fNet = sum [f (t,x0,v0) | f <- fs] + acc = fNet / m + in (1,v0,acc) + +euler1D :: R -- time step dt + -> (State1D -> (R,R,R)) -- differential equation + -> State1D -> State1D -- state-update function +euler1D dt deriv (t0,x0,v0) + = let (_, _, dvdt) = deriv (t0,x0,v0) + t1 = t0 + dt + x1 = x0 + v0 * dt + v1 = v0 + dvdt * dt + in (t1,x1,v1) + +updateTXV :: R -- time interval dt + -> Mass + -> [State1D -> Force] -- list of force funcs + -> State1D -> State1D -- state-update function +updateTXV dt m fs = euler1D dt (newtonSecond1D m fs) + +statesTXV :: R -- time step + -> Mass + -> State1D -- initial state + -> [State1D -> Force] -- list of force funcs + -> [State1D] -- infinite list of states +statesTXV dt m txv0 fs = iterate (updateTXV dt m fs) txv0 + +-- assume that dt is the same between adjacent pairs +velocity1D :: [State1D] -- infinite list + -> Time -> Velocity -- velocity function +velocity1D sts t + = let (t0,_,_) = sts !! 0 + (t1,_,_) = sts !! 1 + dt = t1 - t0 + numSteps = abs $ round (t / dt) + (_,_,v0) = sts !! numSteps + in v0 + +velocityFtxv :: R -- time step + -> Mass + -> State1D -- initial state + -> [State1D -> Force] -- list of force funcs + -> Time -> Velocity -- velocity function +velocityFtxv dt m txv0 fs = velocity1D (statesTXV dt m txv0 fs) + +-- assume that dt is the same between adjacent pairs +position1D :: [State1D] -- infinite list + -> Time -> Position -- position function +position1D sts t + = let (t0,_,_) = sts !! 0 + (t1,_,_) = sts !! 1 + dt = t1 - t0 + numSteps = abs $ round (t / dt) + (_,x0,_) = sts !! numSteps + in x0 + +positionFtxv :: R -- time step + -> Mass + -> State1D -- initial state + -> [State1D -> Force] -- list of force funcs + -> Time -> Position -- position function +positionFtxv dt m txv0 fs = position1D (statesTXV dt m txv0 fs) + +springForce :: R -> State1D -> Force +springForce k (_,x0,_) = -k * x0 + +dampedHOForces :: [State1D -> Force] +dampedHOForces = [springForce 0.8 + ,\(_,_,v0) -> fAir 2 1.225 (pi * 0.02**2) v0 + ,\_ -> -0.0027 * 9.80665 + ] + +dampedHOStates :: [State1D] +dampedHOStates = statesTXV 0.001 0.0027 (0.0,0.1,0.0) dampedHOForces + +dampedHOGraph :: IO () +dampedHOGraph + = plotPath [Title "Ping Pong Ball on a Slinky" + ,XLabel "Time (s)" + ,YLabel "Position (m)" + ,PNG "dho.png" + ,Key Nothing + ] [(t,x) | (t,x,_) <- take 3000 dampedHOStates] + +pingpongPosition :: Time -> Velocity +pingpongPosition = positionFtxv 0.001 0.0027 (0,0.1,0) dampedHOForces + +dampedHOGraph2 :: IO () +dampedHOGraph2 + = plotFunc [Title "Ping Pong Ball on a Slinky" + ,XLabel "Time (s)" + ,YLabel "Position (m)" + ,Key Nothing + ] [0,0.01..3] pingpongPosition + +pingpongVelocity :: Time -> Velocity +pingpongVelocity = velocityFtxv 0.001 0.0027 (0,0.1,0) dampedHOForces + +dampedHOGraph3 :: IO () +dampedHOGraph3 + = plotFunc [Title "Ping Pong Ball on a Slinky" + ,XLabel "Time (s)" + ,YLabel "Velocity (m/s)" + ,PNG "dho2.png" + ,Key Nothing + ] [0,0.01..3] pingpongVelocity + +eulerCromer1D :: R -- time step dt + -> (State1D -> (R,R,R)) -- differential equation + -> State1D -> State1D -- state-update function +eulerCromer1D dt deriv (t0,x0,v0) + = let (_, _, dvdt) = deriv (t0,x0,v0) + t1 = t0 + dt + x1 = x0 + v1 * dt + v1 = v0 + dvdt * dt + in (t1,x1,v1) + +updateTXVEC :: R -- time interval dt + -> Mass + -> [State1D -> Force] -- list of force funcs + -> State1D -> State1D -- state-update function +updateTXVEC dt m fs = eulerCromer1D dt (newtonSecond1D m 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 method = iterate . 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 + (dtdt0, dxdt0, dvdt0) +++ (dtdt1, dxdt1, dvdt1) + = (dtdt0 + dtdt1, dxdt0 + dxdt1, dvdt0 + dvdt1) + scale w (dtdt0, dxdt0, dvdt0) = (w * dtdt0, w * dxdt0, w * 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 dt (dtdt,dxdt,dvdt) (t,x,v) + = (t + dtdt * dt, x + dxdt * dt, v + dvdt * 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 dt deriv st0 = shift dt (deriv st0) 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 dt deriv st0 + = let m0 = deriv st0 + m1 = deriv (shift (dt/2) m0 st0) + m2 = deriv (shift (dt/2) m1 st0) + m3 = deriv (shift dt m2 st0) + in shift (dt/6) (m0 +++ m1 +++ m1 +++ m2 +++ m2 +++ m3) st0 + +exponential :: DifferentialEquation (R,R,R) (R,R,R) +exponential (_,x0,v0) = (1,v0,x0) + +update2 :: (R,R,R) -- starting state + -> (R,R,R) -- ending state +update2 = undefined + +earthGravity :: Mass -> State1D -> Force +earthGravity m _ = let g = 9.80665 + in -m * g + +type MState = (Time,Mass,Position,Velocity) + +earthGravity2 :: MState -> Force +earthGravity2 (_,m,_,_) = let g = 9.80665 + in -m * g + +positionFtxv2 :: R -- time step + -> MState -- initial state + -> [MState -> Force] -- list of force funcs + -> Time -> Position -- position function +positionFtxv2 = undefined + +statesTXV2 :: R -- time step + -> MState -- initial state + -> [MState -> Force] -- list of force funcs + -> [MState] -- infinite list of states +statesTXV2 = undefined + +updateTXV2 :: R -- dt for stepping + -> [MState -> Force] -- list of force funcs + -> MState -- current state + -> MState -- new state +updateTXV2 = undefined + +instance RealVectorSpace (R,R) where + (dtdt0, dvdt0) +++ (dtdt1, dvdt1) = (dtdt0 + dtdt1, dvdt0 + dvdt1) + scale w (dtdt0, dvdt0) = (w * dtdt0, w * dvdt0) + +instance Diff (Time,Velocity) (R,R) where + shift dt (dtdt,dvdt) (t,v) + = (t + dtdt * dt, v + dvdt * dt) + +updateTV' :: R -- dt for stepping + -> Mass + -> [(Time,Velocity) -> Force] -- list of force funcs + -> (Time,Velocity) -- current state + -> (Time,Velocity) -- new state +updateTV' = undefined + +forces :: R -> [State1D -> R] +forces mu = [\(_t,x,_v) -> undefined x + ,\(_t,x, v) -> undefined mu x v] + +vdp :: R -> [(R,R)] +vdp mu = map (\(_,x,v) -> (x,v)) $ take 10000 $ + solver (rungeKutta4 0.01) (newtonSecond1D 1 $ forces mu) (0,2,0) + +vdpPhasePlanePlot :: IO () +vdpPhasePlanePlot = plotPaths [Title "Van der Pol oscillator" + ,XLabel "x" + ,YLabel "v" + ,PNG "VanderPol.png" + ,Key Nothing] (undefined :: [[(R,R)]]) diff --git a/src/LPFP/Mechanics3D.hs b/src/LPFP/Mechanics3D.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/Mechanics3D.hs @@ -0,0 +1,537 @@ +{-# OPTIONS -Wall #-} +{-# LANGUAGE MultiParamTypeClasses #-} + +{- | +Module : LPFP.Mechanics3D +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapters 16, 17, and 18 of the book Learn Physics with Functional Programming +-} + +module LPFP.Mechanics3D where + +import LPFP.SimpleVec + ( R, Vec, PosVec, (^+^), (^-^), (*^), (^*), (^/), (<.>), (><) + , vec, sumV, magnitude, zeroV, xComp, yComp, zComp, iHat, jHat, kHat) +import LPFP.Mechanics1D + ( RealVectorSpace(..), Diff(..), NumericalMethod + , Time, TimeStep, rungeKutta4, solver ) +import SpatialMath + ( V3(..), Euler(..) ) +import Graphics.Gnuplot.Simple + ( Attribute(..), Aspect(..), plotFunc, plotPaths ) +import qualified Graphics.Gloss as G +import qualified Vis as V + +-- | Data type for the state of a single particle in three-dimensional space. +data ParticleState = ParticleState { mass :: R + , charge :: R + , time :: R + , posVec :: Vec + , velocity :: Vec } + deriving Show + +-- | A default particle state. +defaultParticleState :: ParticleState +defaultParticleState = ParticleState { mass = 1 + , charge = 0 + , time = 0 + , posVec = zeroV + , velocity = zeroV } + +rockState :: ParticleState +rockState + = defaultParticleState { mass = 2 -- kg + , velocity = 3 *^ iHat ^+^ 4 *^ 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 { dmdt :: R + , dqdt :: R + , dtdt :: R + , drdt :: Vec + , dvdt :: Vec } + deriving Show + +-- | Given a list of forces, return a differential equation +-- based on Newton's second law. +newtonSecondPS :: [OneBodyForce] + -> ParticleState -> DParticleState -- ^ a differential equation +newtonSecondPS fs st + = let fNet = sumV [f st | f <- fs] + m = mass st + v = velocity st + acc = fNet ^/ m + in DParticleState { dmdt = 0 -- dm/dt + , dqdt = 0 -- dq/dt + , dtdt = 1 -- dt/dt + , drdt = v -- dr/dt + , dvdt = acc -- dv/dt + } + +-- | The force of gravity near Earth's surface. +-- The z direction is toward the sky. +-- Assumes SI units. +earthSurfaceGravity :: OneBodyForce +earthSurfaceGravity st + = let g = 9.80665 -- m/s^2 + in (-mass st * g) *^ 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 m _q _t r _v) + = let bigG = 6.67408e-11 -- N m^2/kg^2 + sunMass = 1.98848e30 -- kg + in (-bigG * sunMass * m) *^ r ^/ magnitude 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 drag rho area (ParticleState _m _q _t _r v) + = (-0.5 * drag * rho * area * magnitude v) *^ 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 vWind drag rho area (ParticleState _m _q _t _r v) + = let vRel = v ^-^ vWind + in (-0.5 * drag * rho * area * magnitude vRel) *^ vRel + +-- | The force of uniform electric and magnetic fields on an object. +uniformLorentzForce :: Vec -- ^ E + -> Vec -- ^ B + -> OneBodyForce +uniformLorentzForce vE vB (ParticleState _m q _t _r v) + = q *^ (vE ^+^ v >< vB) + +-- | Euler-Cromer method for the 'ParticleState' data type. +eulerCromerPS :: TimeStep -- dt for stepping + -> NumericalMethod ParticleState DParticleState +eulerCromerPS dt deriv st + = let t = time st + r = posVec st + v = velocity st + dst = deriv st + acc = dvdt dst + v' = v ^+^ acc ^* dt + in st { time = t + dt + , posVec = r ^+^ v' ^* dt + , velocity = v ^+^ acc ^* dt + } + +instance RealVectorSpace DParticleState where + dst1 +++ dst2 + = DParticleState { dmdt = dmdt dst1 + dmdt dst2 + , dqdt = dqdt dst1 + dqdt dst2 + , dtdt = dtdt dst1 + dtdt dst2 + , drdt = drdt dst1 ^+^ drdt dst2 + , dvdt = dvdt dst1 ^+^ dvdt dst2 + } + scale w dst + = DParticleState { dmdt = w * dmdt dst + , dqdt = w * dqdt dst + , dtdt = w * dtdt dst + , drdt = w *^ drdt dst + , dvdt = w *^ dvdt dst + } + +instance Diff ParticleState DParticleState where + shift dt dps (ParticleState m q t r v) + = ParticleState (m + dmdt dps * dt) + (q + dqdt dps * dt) + (t + dtdt dps * dt) + (r ^+^ drdt dps ^* dt) + (v ^+^ dvdt dps ^* 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 method = iterate . method . 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 method = method . 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 method fs st t + = let states = statesPS method fs st + dt = time (states !! 1) - time (states !! 0) + numSteps = abs $ round (t / dt) + st1 = solver method (newtonSecondPS fs) st !! numSteps + in posVec st1 + +-- | Given a time-scale factor, +-- an animation rate, +-- an initial state, +-- a display function, +-- and an update function, +-- use gloss to produce an animation. +simulateGloss :: R -- ^ time-scale factor + -> Int -- ^ animation rate + -> s -- ^ initial state + -> (s -> G.Picture) -- ^ display function + -> (TimeStep -> s -> s) -- ^ update function + -> IO () +simulateGloss tsFactor rate initialState picFunc updateFunc + = G.simulate (G.InWindow "" (1000, 750) (10, 10)) G.black rate + initialState picFunc + (\_ -> updateFunc . (* tsFactor) . realToFrac) + +-- | Given a time-scale factor, +-- an animation rate, +-- an initial state, +-- a display function, +-- and an update function, +-- use Vis (not-gloss) to produce an animation. +simulateVis :: HasTime s => R -- ^ time-scale factor + -> Int -- ^ animation rate + -> s -- ^ initial state + -> (s -> V.VisObject R) + -> (TimeStep -> s -> s) + -> IO () +simulateVis tsFactor rate initialState picFunc updateFunc + = let visUpdateFunc ta st + = let dtp = tsFactor * realToFrac ta - timeOf st + in updateFunc dtp st + in V.simulate V.defaultOpts (1/fromIntegral rate) + initialState (orient . picFunc) visUpdateFunc + +v3FromVec :: Vec -> V3 R +v3FromVec v = V3 x y z + where + x = xComp v + y = yComp v + z = zComp v + +orient :: V.VisObject R -> V.VisObject R +orient pict = V.RotEulerDeg (Euler 270 180 0) $ pict + +class HasTime s where + timeOf :: s -> Time + +instance HasTime ParticleState where + timeOf = time + +constantForce :: Vec -> OneBodyForce +constantForce f = undefined f + +moonSurfaceGravity :: OneBodyForce +moonSurfaceGravity = undefined + +earthGravity :: OneBodyForce +earthGravity = undefined + +tvyPair :: ParticleState -> (R,R) +tvyPair st = undefined st + +tvyPairs :: [ParticleState] -> [(R,R)] +tvyPairs sts = undefined sts + +tle1yr :: ParticleState -> Bool +tle1yr st = undefined st + +stateFunc :: [ParticleState] + -> Time -> ParticleState +stateFunc sts t + = let t0 = undefined sts + t1 = undefined sts + dt = undefined t0 t1 + numSteps = undefined t dt + in undefined sts numSteps + +airResAtAltitude :: R -- ^ drag coefficient + -> R -- ^ air density at sea level + -> R -- ^ cross-sectional area of object + -> OneBodyForce +airResAtAltitude drag rho0 area (ParticleState _m _q _t r v) + = undefined drag rho0 area r v + +projectileRangeComparison :: R -> R -> (R,R,R) +projectileRangeComparison v0 thetaDeg + = let vx0 = v0 * cos (thetaDeg / 180 * pi) + vz0 = v0 * sin (thetaDeg / 180 * pi) + drag = 1 + ballRadius = 0.05 -- meters + area = pi * ballRadius**2 + airDensity = 1.225 -- kg/m^3 @ sea level + leadDensity = 11342 -- kg/m^3 + m = leadDensity * 4 * pi * ballRadius**3 / 3 + stateInitial = undefined m vx0 vz0 + aboveSeaLevel :: ParticleState -> Bool + aboveSeaLevel st = zComp (posVec st) >= 0 + range :: [ParticleState] -> R + range = xComp . posVec . last . takeWhile aboveSeaLevel + method = rungeKutta4 0.01 + forcesNoAir + = [earthSurfaceGravity] + forcesConstAir + = [earthSurfaceGravity, airResistance drag airDensity area] + forcesVarAir + = [earthSurfaceGravity, airResAtAltitude drag airDensity area] + rangeNoAir = range $ statesPS method forcesNoAir stateInitial + rangeConstAir = range $ statesPS method forcesConstAir stateInitial + rangeVarAir = range $ statesPS method forcesVarAir stateInitial + in undefined rangeNoAir rangeConstAir rangeVarAir + +halleyUpdate :: TimeStep + -> ParticleState -> ParticleState +halleyUpdate dt + = updatePS (eulerCromerPS dt) [sunGravity] + +halleyInitial :: ParticleState +halleyInitial = ParticleState { mass = 2.2e14 -- kg + , charge = 0 + , time = 0 + , posVec = 8.766e10 *^ iHat -- m + , velocity = 54569 *^ jHat } -- m/s + +disk :: Float -> G.Picture +disk radius = G.ThickCircle (radius/2) radius + +baseballForces :: [OneBodyForce] +baseballForces + = let area = pi * (0.074 / 2) ** 2 + in [earthSurfaceGravity + ,airResistance 0.3 1.225 area] + +baseballTrajectory :: R -- time step + -> R -- initial speed + -> R -- launch angle in degrees + -> [(R,R)] -- (y,z) pairs +baseballTrajectory dt v0 thetaDeg + = let thetaRad = thetaDeg * pi / 180 + vy0 = v0 * cos thetaRad + vz0 = v0 * sin thetaRad + initialState + = ParticleState { mass = 0.145 + , charge = 0 + , time = 0 + , posVec = zeroV + , velocity = vec 0 vy0 vz0 } + in trajectory $ zGE0 $ + statesPS (eulerCromerPS dt) baseballForces initialState + +zGE0 :: [ParticleState] -> [ParticleState] +zGE0 = takeWhile (\(ParticleState _ _ _ r _) -> zComp r >= 0) + +trajectory :: [ParticleState] -> [(R,R)] +trajectory sts = [(yComp r,zComp r) | (ParticleState _ _ _ r _) <- sts] + +baseballRange :: R -- time step + -> R -- initial speed + -> R -- launch angle in degrees + -> R -- range +baseballRange dt v0 thetaDeg + = let (y,_) = last $ baseballTrajectory dt v0 thetaDeg + in y + +baseballRangeGraph :: IO () +baseballRangeGraph + = plotFunc [Title "Range for baseball hit at 45 m/s" + ,XLabel "Angle above horizontal (degrees)" + ,YLabel "Horizontal range (m)" + ,PNG "baseballrange.png" + ,Key Nothing + ] [10,11..80] $ baseballRange 0.01 45 + +bestAngle :: (R,R) +bestAngle + = maximum [(baseballRange 0.01 45 thetaDeg,thetaDeg) | + thetaDeg <- [30,31..60]] + +projectileUpdate :: TimeStep + -> ParticleState -- old state + -> ParticleState -- new state +projectileUpdate dt + = updatePS (eulerCromerPS dt) baseballForces + +projectileInitial :: [String] -> ParticleState +projectileInitial [] = error "Please supply initial speed and angle." +projectileInitial [_] = error "Please supply initial speed and angle." +projectileInitial (_:_:_:_) + = error "First argument is speed. Second is angle in degrees." +projectileInitial (arg1:arg2:_) + = let v0 = read arg1 :: R -- initial speed, m/s + angleDeg = read arg2 :: R -- initial angle, degrees + theta = angleDeg * pi / 180 -- in radians + in defaultParticleState + { mass = 0.145 -- kg + , posVec = zeroV + , velocity = vec 0 (v0 * cos theta) (v0 * sin theta) + } + +protonUpdate :: TimeStep -> ParticleState -> ParticleState +protonUpdate dt + = updatePS (rungeKutta4 dt) [uniformLorentzForce zeroV (3e-8 *^ kHat)] + +protonInitial :: ParticleState +protonInitial + = defaultParticleState { mass = 1.672621898e-27 -- kg + , charge = 1.602176621e-19 -- C + , posVec = zeroV + , velocity = 1.5*^jHat ^+^ 0.3*^kHat -- m/s + } + +protonPicture :: ParticleState -> V.VisObject R +protonPicture st + = let r0 = v3FromVec (posVec st) + in V.Trans r0 (V.Sphere 0.1 V.Solid V.red) + +apR :: R +apR = 0.04 -- meters + +wallForce :: OneBodyForce +wallForce ps + = let m = mass ps + r = posVec ps + x = xComp r + y = yComp r + z = zComp r + v = velocity ps + timeStep = 5e-4 / 60 + in if y >= 1 && y < 1.1 && sqrt (x**2 + z**2) > apR + then (-m) *^ (v ^/ timeStep) + else zeroV + +zOut :: V.VisObject R -> V.VisObject R +zOut = V.RotEulerDeg (Euler 90 0 90) + +energy :: ParticleState -> R +energy ps = undefined ps + +firstOrbit :: ParticleState -> Bool +firstOrbit st + = let year = 365.25 * 24 * 60 * 60 + in time st < 50 * year || yComp (posVec st) <= 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 fs st + = let fNet = sumV [f st | f <- fs] + c = 299792458 -- m / s + m = mass st + v = velocity st + u = v ^/ c + acc = sqrt (1 - u <.> u) *^ (fNet ^-^ (fNet <.> u) *^ u) ^/ m + in DParticleState { dmdt = 0 -- dm/dt + , dqdt = 0 -- dq/dt + , dtdt = 1 -- dt/dt + , drdt = v -- dr/dt + , dvdt = acc -- dv/vt + } + +constantForcePlot :: IO () +constantForcePlot + = let year = 365.25 * 24 * 60 * 60 -- seconds + c = 299792458 -- m/s + method = rungeKutta4 1000 + forces = [const (10 *^ iHat)] + initialState = defaultParticleState { mass = 1 } + newtonStates = solver method (newtonSecondPS forces) initialState + relativityStates = solver method (relativityPS forces) initialState + newtonTVs = [(time st / year, xComp (velocity st) / c) + | st <- takeWhile tle1yr newtonStates] + relativityTVs = [(time st / year, xComp (velocity st) / c) + | st <- takeWhile tle1yr relativityStates] + in plotPaths [Key Nothing + ,Title "Response to a constant force" + ,XLabel "Time (years)" + ,YLabel "Velocity (multiples of c)" + ,PNG "constantForceComp.png" + ,customLabel (0.1,1) "mass = 1 kg" + ,customLabel (0.1,0.9) "force = 10 N" + ,customLabel (0.5,0.7) "Newtonian" + ,customLabel (0.8,0.6) "relativistic" + ] [newtonTVs,relativityTVs] + +customLabel :: (R,R) -> String -> Attribute +customLabel (x,y) label + = Custom "label" + ["\"" ++ label ++ "\"" ++ " at " ++ show x ++ "," ++ show y] + +circularPlot :: IO () +circularPlot + = let c = 299792458 -- m/s + method = rungeKutta4 1e-9 + forces = [uniformLorentzForce zeroV kHat] -- 1 T + initialState = defaultParticleState + { mass = 1.672621898e-27 -- kg + , charge = 1.602176621e-19 -- C + , velocity = 0.8 *^ c *^ jHat + } + newtonStates = solver method (newtonSecondPS forces) initialState + relativityStates = solver method (relativityPS forces) initialState + newtonXYs = [(xComp (posVec st), yComp (posVec st)) + | st <- take 100 newtonStates] + relativityXYs = [(xComp (posVec st), yComp (posVec st)) + | st <- take 120 relativityStates] + in plotPaths [Key Nothing + ,Aspect (Ratio 1) + ,Title "Proton in a 1-T magnetic field" + ,XLabel "x (m)" + ,YLabel "y (m)" + ,PNG "circularComp.png" + ,customLabel (0.5,4.5) "v = 0.8 c" + ,customLabel (2.5,0.0) "Newtonian" + ,customLabel (3.0,3.5) "relativistic" + ] [newtonXYs,relativityXYs] + +twoProtUpdate :: TimeStep + -> (ParticleState,ParticleState) + -> (ParticleState,ParticleState) +twoProtUpdate dt (stN,stR) + = let forces = [uniformLorentzForce zeroV kHat] + in (rungeKutta4 dt (newtonSecondPS forces) stN + ,rungeKutta4 dt (relativityPS forces) stR) + +twoProtInitial :: (ParticleState,ParticleState) +twoProtInitial + = let c = 299792458 -- m/s + pInit = protonInitial { velocity = 0.8 *^ c *^ jHat } + in (pInit,pInit) + +twoProtPicture :: (ParticleState,ParticleState) -> G.Picture +twoProtPicture (stN,stR) + = G.scale 50 50 $ G.pictures [G.translate xN yN protonNewtonian + ,G.translate xR yR protonRelativistic] + where + xN = realToFrac $ xComp $ posVec stN + yN = realToFrac $ yComp $ posVec stN + xR = realToFrac $ xComp $ posVec stR + yR = realToFrac $ yComp $ posVec stR + protonNewtonian = G.Color G.blue (disk 0.1) + protonRelativistic = G.Color G.red (disk 0.1) + +relativityPS' :: R -- c + -> [OneBodyForce] + -> ParticleState -> DParticleState +relativityPS' c fs st = undefined c fs st diff --git a/src/LPFP/MultipleObjects.hs b/src/LPFP/MultipleObjects.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/MultipleObjects.hs @@ -0,0 +1,180 @@ +{-# OPTIONS -Wall #-} +{-# LANGUAGE MultiParamTypeClasses #-} + +{- | +Module : LPFP.MultipleObjects +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 19 of the book Learn Physics with Functional Programming +-} + +module LPFP.MultipleObjects where + +import LPFP.SimpleVec + ( Vec, R, (^+^), (^-^), (*^), (^*), (^/), zeroV, magnitude ) +import LPFP.Mechanics1D + ( RealVectorSpace(..), Diff(..), NumericalMethod, Mass, TimeStep, euler ) +import LPFP.Mechanics3D + ( OneBodyForce, ParticleState(..), DParticleState(..), HasTime(..) + , defaultParticleState, newtonSecondPS ) + +type TwoBodyForce + = ParticleState -- force is produced BY particle with this state + -> ParticleState -- force acts ON particle with this state + -> ForceVector + +type ForceVector = Vec + +oneFromTwo :: ParticleState -- state of particle PRODUCING the force + -> TwoBodyForce + -> OneBodyForce +oneFromTwo stBy f = f stBy + +gravityMagnitude :: Mass -> Mass -> R -> R +gravityMagnitude m1 m2 r = let gg = 6.67408e-11 -- N m^2 / kg^2 + in gg * m1 * m2 / r**2 + +universalGravity :: TwoBodyForce +universalGravity st1 st2 + = let gg = 6.67408e-11 -- N m^2 / kg^2 + m1 = mass st1 + m2 = mass st2 + r1 = posVec st1 + r2 = posVec st2 + r21 = r2 ^-^ r1 + in (-gg) *^ m1 *^ m2 *^ r21 ^/ magnitude r21 ** 3 + +constantRepulsiveForceWrong :: ForceVector -> TwoBodyForce +constantRepulsiveForceWrong force = \_ _ -> force + +constantRepulsiveForce :: R -> TwoBodyForce +constantRepulsiveForce force st1 st2 + = let r1 = posVec st1 + r2 = posVec st2 + r21 = r2 ^-^ r1 + in force *^ r21 ^/ magnitude r21 + +linearSpring :: R -- spring constant + -> R -- equilibrium length + -> TwoBodyForce +linearSpring k re st1 st2 + = let r1 = posVec st1 + r2 = posVec st2 + r21 = r2 ^-^ r1 + r21mag = magnitude r21 + in (-k) *^ (r21mag - re) *^ r21 ^/ r21mag + +-- | Force provided by a spring that is fixed at one end. +fixedLinearSpring :: R -> R -> Vec -> OneBodyForce +fixedLinearSpring k re r1 + = oneFromTwo (defaultParticleState { posVec = r1 }) (linearSpring k re) + +centralForce :: (R -> R) -> TwoBodyForce +centralForce f st1 st2 + = let r1 = posVec st1 + r2 = posVec st2 + r21 = r2 ^-^ r1 + r21mag = magnitude r21 + in f r21mag *^ r21 ^/ r21mag + +linearSpringCentral :: R -- spring constant + -> R -- equilibrium length + -> TwoBodyForce +linearSpringCentral k re = centralForce (\r -> -k * (r - re)) + +billiardForce :: R -- spring constant + -> R -- threshold center separation + -> TwoBodyForce +billiardForce k re + = centralForce $ \r -> if r >= re + then 0 + else (-k * (r - re)) + +data Force = ExternalForce Int OneBodyForce + | InternalForce Int Int TwoBodyForce + +data MultiParticleState + = MPS { particleStates :: [ParticleState] } deriving Show + +instance HasTime MultiParticleState where + timeOf (MPS sts) = time (sts !! 0) + +data DMultiParticleState = DMPS [DParticleState] deriving Show + +newtonSecondMPS :: [Force] + -> MultiParticleState -> DMultiParticleState -- a diff eqn + +newtonSecondMPS fs mpst@(MPS sts) + = let deriv (n,st) = newtonSecondPS (forcesOn n mpst fs) st + in DMPS $ map deriv (zip [0..] sts) + +forcesOn :: Int -> MultiParticleState -> [Force] -> [OneBodyForce] +forcesOn n mpst = map (forceOn n mpst) + +forceOn :: Int -> MultiParticleState -> Force -> OneBodyForce +forceOn n _ (ExternalForce n0 fOneBody) + | n == n0 = fOneBody + | otherwise = const zeroV +forceOn n (MPS sts) (InternalForce n0 n1 fTwoBody) + | n == n0 = oneFromTwo (sts !! n1) fTwoBody -- n1 acts on n0 + | n == n1 = oneFromTwo (sts !! n0) fTwoBody -- n0 acts on n1 + | otherwise = const zeroV + +instance RealVectorSpace DMultiParticleState where + DMPS dsts1 +++ DMPS dsts2 = DMPS $ zipWith (+++) dsts1 dsts2 + scale w (DMPS dsts) = DMPS $ map (scale w) dsts + +instance Diff MultiParticleState DMultiParticleState where + shift dt (DMPS dsts) (MPS sts) = MPS $ zipWith (shift dt) dsts sts + +eulerCromerMPS :: TimeStep -- dt for stepping + -> NumericalMethod MultiParticleState DMultiParticleState +eulerCromerMPS dt deriv mpst0 + = let mpst1 = euler dt deriv mpst0 + sts0 = particleStates mpst0 + sts1 = particleStates mpst1 + -- now update positions + in MPS $ [ st1 { posVec = posVec st0 ^+^ velocity st1 ^* dt } + | (st0,st1) <- zip sts0 sts1 ] + +updateMPS :: NumericalMethod MultiParticleState DMultiParticleState + -> [Force] + -> MultiParticleState -> MultiParticleState +updateMPS method = method . newtonSecondMPS + +statesMPS :: NumericalMethod MultiParticleState DMultiParticleState + -> [Force] + -> MultiParticleState -> [MultiParticleState] +statesMPS method = iterate . method . newtonSecondMPS + +speed :: ParticleState -> R +speed st = undefined st + +universalGravity' :: TwoBodyForce +universalGravity' (ParticleState m1 _ _ r1 _) (ParticleState m2 _ _ r2 _) + = undefined m1 r1 m2 r2 + +universalGravityCentral :: TwoBodyForce +universalGravityCentral = undefined + +lennardJones :: R -- dissociation energy + -> R -- equilibrium length + -> TwoBodyForce +lennardJones de re = centralForce $ \r -> undefined de re r + +systemKE :: MultiParticleState -> R +systemKE mpst = undefined mpst + +forcesOn' :: Int -> MultiParticleState -> [Force] -> [OneBodyForce] +forcesOn' n mpst fs = externalForcesOn n fs ++ internalForcesOn n mpst fs + +externalForcesOn :: Int -> [Force] -> [OneBodyForce] +externalForcesOn n fs = undefined n fs + +internalForcesOn :: Int -> MultiParticleState -> [Force] -> [OneBodyForce] +internalForcesOn n (MPS sts) fs + = [oneFromTwo (sts !! n1) f | InternalForce n0 n1 f <- fs, n == n0] ++ + [oneFromTwo (sts !! n0) f | InternalForce n0 n1 f <- fs, n == n1] diff --git a/src/LPFP/Newton2.hs b/src/LPFP/Newton2.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/Newton2.hs @@ -0,0 +1,230 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFP.Newton2 +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 14 of the book Learn Physics with Functional Programming +-} + +module LPFP.Newton2 where + +import Graphics.Gnuplot.Simple + +velocityCF :: Mass + -> Velocity -- initial velocity + -> [Force] -- list of forces + -> Time -> Velocity -- velocity function + +type R = Double + +type Mass = R +type Time = R +type Position = R +type Velocity = R +type Force = R + +velocityCF m v0 fs + = let fNet = sum fs -- net force + a0 = fNet / m -- Newton's second law + v t = v0 + a0 * t -- constant acceleration eqn + in v + +positionCF :: Mass + -> Position -- initial position + -> Velocity -- initial velocity + -> [Force] -- list of forces + -> Time -> Position -- position function +positionCF m x0 v0 fs + = let fNet = sum fs + a0 = fNet / m + x t = x0 + v0 * t + a0*t**2 / 2 + in x + +carGraph :: IO () +carGraph + = plotFunc [Title "Car on an air track" + ,XLabel "Time (s)" + ,YLabel "Velocity of Car (m/s)" + ,PNG "CarVelocity.png" + ,Key Nothing + ] [0..4 :: Time] (velocityCF 0.1 0.6 [0.04, -0.08]) + +velocityFt :: R -- dt for integral + -> Mass + -> Velocity -- initial velocity + -> [Time -> Force] -- list of force functions + -> Time -> Velocity -- velocity function +velocityFt dt m v0 fs + = let fNet t = sum [f t | f <- fs] + a t = fNet t / m + in antiDerivative dt v0 a + +-- | Given a step size, a y-intercept, and a function, return a function +-- with the given y-intercept whose +-- derivative is the given function. +antiDerivative :: R -> R -> (R -> R) -> (R -> R) +antiDerivative dt v0 a t = v0 + integral dt a 0 t + +-- | Given a step size, a function, a lower limit, and an upper limit, return +-- the definite integral of the function. +integral :: R -> (R -> R) -> R -> R -> R +integral dt f a b + = sum [f t * dt | t <- [a+dt/2, a+3*dt/2 .. b - dt/2]] + +positionFt :: R -- dt for integral + -> Mass + -> Position -- initial position + -> Velocity -- initial velocity + -> [Time -> Force] -- list of force functions + -> Time -> Position -- position function +positionFt dt m x0 v0 fs + = antiDerivative dt x0 (velocityFt dt m v0 fs) + +pedalCoast :: Time -> Force +pedalCoast t + = let tCycle = 20 + nComplete :: Int + nComplete = truncate (t / tCycle) + remainder = t - fromIntegral nComplete * tCycle + in if remainder < 10 + then 10 + else 0 + +childGraph :: IO () +childGraph + = plotFunc [Title "Child pedaling then coasting" + ,XLabel "Time (s)" + ,YLabel "Position of Bike (m)" + ,PNG "ChildPosition.png" + ,Key Nothing + ] [0..40 :: R] (positionFt 0.1 20 0 0 [pedalCoast]) + +fAir :: R -- drag coefficient + -> R -- air density + -> R -- cross-sectional area of object + -> Velocity + -> Force +fAir drag rho area v = -drag * rho * area * abs v * v / 2 + +newtonSecondV :: Mass + -> [Velocity -> Force] -- list of force functions + -> Velocity -- current velocity + -> R -- derivative of velocity +newtonSecondV m fs v0 = sum [f v0 | f <- fs] / m + +updateVelocity :: R -- time interval dt + -> Mass + -> [Velocity -> Force] -- list of force functions + -> Velocity -- current velocity + -> Velocity -- new velocity +updateVelocity dt m fs v0 + = v0 + (newtonSecondV m fs v0) * dt + +velocityFv :: R -- time step + -> Mass + -> Velocity -- initial velocity v(0) + -> [Velocity -> Force] -- list of force functions + -> Time -> Velocity -- velocity function +velocityFv dt m v0 fs t + = let numSteps = abs $ round (t / dt) + in iterate (updateVelocity dt m fs) v0 !! numSteps + +bikeVelocity :: Time -> Velocity +bikeVelocity = velocityFv 1 70 0 [const 100,fAir 2 1.225 0.6] + +bikeGraph :: IO () +bikeGraph = plotFunc [Title "Bike velocity" + ,XLabel "Time (s)" + ,YLabel "Velocity of Bike (m/s)" + ,PNG "BikeVelocity1.png" + ,Key Nothing + ] [0,0.5..60] bikeVelocity + +newtonSecondTV :: Mass + -> [(Time,Velocity) -> Force] -- force funcs + -> (Time,Velocity) -- current state + -> (R,R) -- deriv of state +newtonSecondTV m fs (t,v0) + = let fNet = sum [f (t,v0) | f <- fs] + acc = fNet / m + in (1,acc) + +updateTV :: R -- time interval dt + -> Mass + -> [(Time,Velocity) -> Force] -- list of force funcs + -> (Time,Velocity) -- current state + -> (Time,Velocity) -- new state +updateTV dt m fs (t,v0) + = let (dtdt, dvdt) = newtonSecondTV m fs (t,v0) + in (t + dtdt * dt + ,v0 + dvdt * dt) + +statesTV :: R -- time step + -> Mass + -> (Time,Velocity) -- initial state + -> [(Time,Velocity) -> Force] -- list of force funcs + -> [(Time,Velocity)] -- infinite list of states +statesTV dt m tv0 fs + = iterate (updateTV dt m fs) tv0 + +velocityFtv :: R -- time step + -> Mass + -> (Time,Velocity) -- initial state + -> [(Time,Velocity) -> Force] -- list of force funcs + -> Time -> Velocity -- velocity function +velocityFtv dt m tv0 fs t + = let numSteps = abs $ round (t / dt) + in snd $ statesTV dt m tv0 fs !! numSteps + +pedalCoastAir :: [(Time,Velocity)] +pedalCoastAir = statesTV 0.1 20 (0,0) + [\(t,_) -> pedalCoast t + ,\(_,v) -> fAir 2 1.225 0.5 v] + +pedalCoastAirGraph :: IO () +pedalCoastAirGraph + = plotPath [Title "Pedaling and coasting with air" + ,XLabel "Time (s)" + ,YLabel "Velocity of Bike (m/s)" + ,PNG "pedalCoastAirGraph.png" + ,Key Nothing + ] (takeWhile (\(t,_) -> t <= 100) + pedalCoastAir) + +pedalCoastAir2 :: Time -> Velocity +pedalCoastAir2 = velocityFtv 0.1 20 (0,0) + [\( t,_v) -> pedalCoast t + ,\(_t, v) -> fAir 1 1.225 0.5 v] + +velocityCF' :: Mass + -> Velocity -- initial velocity + -> [Force] -- list of forces + -> Time -> Velocity -- velocity function +velocityCF' m v0 fs t = undefined m v0 fs t + +sumF :: [R -> R] -> R -> R +sumF = undefined + +positionFv :: R -- time step + -> Mass + -> Position -- initial position x(0) + -> Velocity -- initial velocity v(0) + -> [Velocity -> Force] -- list of force functions + -> Time -> Position -- position function +positionFv = undefined + +positionFtv :: R -- time step + -> Mass + -> Position -- initial position x(0) + -> Velocity -- initial velocity v(0) + -> [(Time,Velocity) -> Force] -- force functions + -> Time -> Position -- position function +positionFtv = undefined + +updateExample :: (Time,Velocity) -- starting state + -> (Time,Velocity) -- ending state +updateExample = undefined diff --git a/src/LPFP/SimpleVec.hs b/src/LPFP/SimpleVec.hs new file mode 100644 --- /dev/null +++ b/src/LPFP/SimpleVec.hs @@ -0,0 +1,256 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFP.SimpleVec +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 10 of the book Learn Physics with Functional Programming +-} + +module LPFP.SimpleVec where + +infixl 6 ^+^ +infixl 6 ^-^ +infixr 7 *^ +infixl 7 ^* +infixr 7 ^/ +infixr 7 <.> +infixl 7 >< + +-- | A vector derivative takes a vector-valued function of a real variable (usually time) as input, +-- and produces a vector-valued function of a real variable as output. +type VecDerivative = (R -> Vec) -> R -> Vec + +-- | Given a step size, calculate the vector derivative of a vector-valued function of a real variable +-- (usually time). +vecDerivative :: R -> VecDerivative +vecDerivative dt v t = (v (t + dt/2) ^-^ v (t - dt/2)) ^/ dt + +v1 :: R -> Vec +v1 t = 2 *^ t**2 *^ iHat ^+^ 3 *^ t**3 *^ jHat ^+^ t**4 *^ kHat + +xCompFunc :: (R -> Vec) -> R -> R +xCompFunc v t = xComp (v t) + +-- | A derivative takes a real-valued function of a real variable (often time) as input, +-- and produces a real-valued function of a real variable as output. +type Derivative = (R -> R) -> R -> R + +-- | Given a step size, calculate the derivative of a real-valued function of a real variable +-- (often time). +derivative :: R -> Derivative +derivative dt x t = (x (t + dt/2) - x (t - dt/2)) / dt + +-- | Time is a real number. +type Time = R +-- | The position of a particle can be represented as a vector. +type PosVec = Vec +-- | Velocity is a vector. +type Velocity = Vec +-- | Acceleration is a vector. +type Acceleration = Vec + +-- | Given a time step and a position function, return a velocity function. +velFromPos :: R -- ^ dt + -> (Time -> PosVec ) -- ^ position function + -> (Time -> Velocity) -- ^ velocity function +velFromPos = vecDerivative + +-- | Given a time step and a velocity function, return an acceleration function. +accFromVel :: R -- dt + -> (Time -> Velocity) -- velocity function + -> (Time -> Acceleration) -- acceleration function +accFromVel = vecDerivative + +-- | Given initial position and a constant velocity, return a position function. +positionCV :: PosVec -> Velocity -> Time -> PosVec +positionCV r0 v0 t = v0 ^* t ^+^ r0 + +-- | Given initial velocity and a constant acceleration, return a velocity function. +velocityCA :: Velocity -> Acceleration -> Time -> Velocity +velocityCA v0 a0 t = a0 ^* t ^+^ v0 + +-- | Given initial position, initial velocity, and a constant acceleration, return a position function. +positionCA :: PosVec -> Velocity -> Acceleration + -> Time -> PosVec +positionCA r0 v0 a0 t = 0.5 *^ t**2 *^ a0 ^+^ v0 ^* t ^+^ r0 + +-- | Given a nonzero velocity and an acceleration, return the component of acceleration +-- parallel to the velocity. +aParallel :: Vec -> Vec -> Vec +aParallel v a = let vHat = v ^/ magnitude v + in (vHat <.> a) *^ vHat + +-- | Given a nonzero velocity and an acceleration, return the component of acceleration +-- perpendicular to the velocity. +aPerp :: Vec -> Vec -> Vec +aPerp v a = a ^-^ aParallel v a + +-- | Given velocity and acceleration, return the rate at which speed is changing. +speedRateChange :: Vec -> Vec -> R +speedRateChange v a = (v <.> a) / magnitude v + +radiusOfCurvature :: Vec -> Vec -> R +radiusOfCurvature v a = (v <.> v) / magnitude (aPerp v a) + +projectilePos :: PosVec -> Velocity -> Time -> PosVec +projectilePos r0 v0 = positionCA r0 v0 (9.81 *^ negateV kHat) + +-- | An approximation to a real number. +type R = Double + +data Mass = Mass R + deriving (Eq,Show) + +data Grade = Grade String Int + deriving (Eq,Show) + +grades :: [Grade] +grades = [Grade "Albert Einstein" 89 + ,Grade "Isaac Newton" 95 + ,Grade "Alan Turing" 91 + ] + +data GradeRecord = GradeRecord { name :: String + , grade :: Int + } deriving (Eq,Show) + +gradeRecords1 :: [GradeRecord] +gradeRecords1 = [GradeRecord "Albert Einstein" 89 + ,GradeRecord "Isaac Newton" 95 + ,GradeRecord "Alan Turing" 91 + ] + +gradeRecords2 :: [GradeRecord] +gradeRecords2 = [GradeRecord {name = "Albert Einstein", grade = 89} + ,GradeRecord {name = "Isaac Newton" , grade = 95} + ,GradeRecord {name = "Alan Turing" , grade = 91} + ] + +data MyBool = MyFalse | MyTrue + deriving (Eq,Show) + +data MyMaybe a = MyNothing + | MyJust a + deriving (Eq,Show) + +-- | A type for three-dimensional vectors. +data Vec = Vec { xComp :: R -- ^ x component of a vector + , yComp :: R -- ^ y component of a vector + , zComp :: R -- ^ z component of a vector + } deriving (Eq) + +instance Show Vec where + show (Vec x y z) = "vec " ++ showDouble x ++ " " + ++ showDouble y ++ " " + ++ showDouble z + +showDouble :: R -> String +showDouble x + | x < 0 = "(" ++ show x ++ ")" + | otherwise = show x + +-- | Form a vector by giving its x, y, and z components. +vec :: R -- ^ x component + -> R -- ^ y component + -> R -- ^ z component + -> Vec +vec = Vec + +-- | A unit vector in the x direction. +iHat :: Vec +iHat = vec 1 0 0 + +-- | A unit vector in the y direction. +jHat :: Vec +jHat = vec 0 1 0 + +-- | A unit vector in the z direction. +kHat :: Vec +kHat = vec 0 0 1 + +-- | The zero vector. +zeroV :: Vec +zeroV = vec 0 0 0 + +-- | Negate a vector. +negateV :: Vec -> Vec +negateV (Vec ax ay az) = Vec (-ax) (-ay) (-az) + +-- | Vector addition. +(^+^) :: Vec -> Vec -> Vec +Vec ax ay az ^+^ Vec bx by bz = Vec (ax+bx) (ay+by) (az+bz) + +-- | Vector subtraction. +(^-^) :: Vec -> Vec -> Vec +Vec ax ay az ^-^ Vec bx by bz = Vec (ax-bx) (ay-by) (az-bz) + +-- | Add a list of vectors. +sumV :: [Vec] -> Vec +sumV = foldr (^+^) zeroV + +-- | Scalar multiplication of a number and a vector. +(*^) :: R -> Vec -> Vec +c *^ Vec ax ay az = Vec (c*ax) (c*ay) (c*az) + +-- | Scalar multiplication of a vector and a number. +(^*) :: Vec -> R -> Vec +Vec ax ay az ^* c = Vec (c*ax) (c*ay) (c*az) + +-- | Dot product of two vectors. +(<.>) :: Vec -> Vec -> R +Vec ax ay az <.> Vec bx by bz = ax*bx + ay*by + az*bz + +-- | Cross product of two vectors. +(><) :: Vec -> Vec -> Vec +Vec ax ay az >< Vec bx by bz + = Vec (ay*bz - az*by) (az*bx - ax*bz) (ax*by - ay*bx) + +-- | Division of a vector by a number. +(^/) :: Vec -> R -> Vec +Vec ax ay az ^/ c = Vec (ax/c) (ay/c) (az/c) + +-- | Magnitude of a vector. +magnitude :: Vec -> R +magnitude v = sqrt(v <.> v) + +-- | Definite integral of a vector-valued function of a real number. +vecIntegral :: R -- ^ step size dt + -> (R -> Vec) -- ^ vector-valued function + -> R -- ^ lower limit + -> R -- ^ upper limit + -> Vec -- ^ result +vecIntegral = undefined + +maxHeight :: PosVec -> Velocity -> R +maxHeight = undefined + +speedCA :: Velocity -> Acceleration -> Time -> R +speedCA = undefined + +xyProj :: Vec -> Vec +xyProj = undefined + +magAngles :: Vec -> (R,R,R) +magAngles = undefined + +gEarth :: Vec +gEarth = undefined + +vBall :: R -> Vec +vBall t = undefined t + +speedRateChangeBall :: R -> R +speedRateChangeBall t = undefined t + +rNCM :: (R, R -> R) -> R -> Vec +rNCM (radius, theta) t = undefined radius theta t + +aPerpFromPosition :: R -> (R -> Vec) -> R -> Vec +aPerpFromPosition epsilon r t + = let v = vecDerivative epsilon r + a = vecDerivative epsilon v + in aPerp (v t) (a t)