LPFP (empty) → 1.0

raw patch · 44 files changed

+4766/−0 lines, 44 filesdep +LPFPdep +basedep +containers

Dependencies added: LPFP, base, containers, diagrams-cairo, diagrams-lib, gloss, gnuplot, linear, not-gloss, spatial-math

Files

• LICENSE +29/−0
• LPFP.cabal +307/−0
• app/ElectricWave.hs +10/−0
• app/GlossAnimate.hs +24/−0
• app/GlossBilliard.hs +10/−0
• app/GlossBilliardAnimate.hs +8/−0
• app/GlossDisplay.hs +14/−0
• app/GlossDisplay2.hs +23/−0
• app/GlossHalley.hs +26/−0
• app/GlossProjectile.hs +28/−0
• app/GlossProton.hs +10/−0
• app/GlossSimulate.hs +31/−0
• app/GlossSimulate2.hs +35/−0
• app/GlossWave.hs +19/−0
• app/GnuplotWave.hs +28/−0
• app/Hydrogen.hs +21/−0
• app/MakeTrajectoryGraph.hs +17/−0
• app/Mech3Proton.hs +6/−0
• app/ProtonEB.hs +13/−0
• app/VisAnimate.hs +13/−0
• app/VisDisplay.hs +11/−0
• app/VisDisplay2.hs +15/−0
• app/VisDisplay3.hs +19/−0
• app/VisProjectile.hs +23/−0
• app/VisSimulate.hs +19/−0
• app/VisTwoSprings.hs +24/−0
• app/hello.hs +2/−0
• src/LPFP.hs +448/−0
• src/LPFP/Charge.hs +123/−0
• src/LPFP/CoordinateSystems.hs +301/−0
• src/LPFP/Current.hs +102/−0
• src/LPFP/ElectricField.hs +316/−0
• src/LPFP/Electricity.hs +99/−0
• src/LPFP/Geometry.hs +134/−0
• src/LPFP/Integrals.hs +143/−0
• src/LPFP/Lorentz.hs +147/−0
• src/LPFP/MOExamples.hs +353/−0
• src/LPFP/MagneticField.hs +117/−0
• src/LPFP/Maxwell.hs +222/−0
• src/LPFP/Mechanics1D.hs +273/−0
• src/LPFP/Mechanics3D.hs +537/−0
• src/LPFP/MultipleObjects.hs +180/−0
• src/LPFP/Newton2.hs +230/−0
• src/LPFP/SimpleVec.hs +256/−0

@@ -0,0 +1,29 @@+Copyright (c) 2022 Scott N. Walck <walck@lvc.edu>.+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.

@@ -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

@@ -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)

@@ -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

@@ -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)

@@ -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)

@@ -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

@@ -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

@@ -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

@@ -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

@@ -0,0 +1,10 @@+{-# OPTIONS -Wall #-}++import LPFP.Mechanics3D+    ( simulateGloss+    , twoProtInitial, twoProtPicture, twoProtUpdate+    )++main :: IO ()+main = simulateGloss 1e-8 20+       twoProtInitial twoProtPicture twoProtUpdate

@@ -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)

@@ -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)

@@ -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

@@ -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)

@@ -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

@@ -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

@@ -0,0 +1,6 @@+{-# OPTIONS -Wall #-}++import LPFP.Mechanics3D (simulateVis, protonInitial, protonPicture, protonUpdate)++main :: IO ()+main = simulateVis 1 60 protonInitial protonPicture protonUpdate

@@ -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 } )

@@ -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)

@@ -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

@@ -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

@@ -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)

@@ -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

@@ -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

@@ -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

@@ -0,0 +1,2 @@+main :: IO ()+main = putStrLn "Hello, world!"

@@ -0,0 +1,448 @@+{-# OPTIONS -Wall #-}+{-# LANGUAGE Trustworthy #-}++{- | +Module      :  LPFP+Copyright   :  (c) Scott N. Walck 2023+License     :  BSD3 (see LICENSE)+Maintainer  :  Scott N. Walck <walck@lvc.edu>+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+    )

@@ -0,0 +1,123 @@+{-# OPTIONS -Wall #-}++{- | +Module      :  LPFP.Charge+Copyright   :  (c) Scott N. Walck 2023+License     :  BSD3 (see LICENSE)+Maintainer  :  Scott N. Walck <walck@lvc.edu>+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

@@ -0,0 +1,301 @@+{-# OPTIONS -Wall #-}++{- | +Module      :  LPFP.CoordinateSystems+Copyright   :  (c) Scott N. Walck 2023+License     :  BSD3 (see LICENSE)+Maintainer  :  Scott N. Walck <walck@lvc.edu>+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

@@ -0,0 +1,102 @@+{-# OPTIONS -Wall #-}++{- | +Module      :  LPFP.Current+Copyright   :  (c) Scott N. Walck 2023+License     :  BSD3 (see LICENSE)+Maintainer  :  Scott N. Walck <walck@lvc.edu>+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

@@ -0,0 +1,316 @@+{-# OPTIONS -Wall #-}++{- | +Module      :  LPFP.ElectricField+Copyright   :  (c) Scott N. Walck 2023+License     :  BSD3 (see LICENSE)+Maintainer  :  Scott N. Walck <walck@lvc.edu>+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

@@ -0,0 +1,99 @@+{-# OPTIONS -Wall #-}++{- | +Module      :  LPFP.Electricity+Copyright   :  (c) Scott N. Walck 2023+License     :  BSD3 (see LICENSE)+Maintainer  :  Scott N. Walck <walck@lvc.edu>+Stability   :  stable++Code from chapter 21 of the book Learn Physics with Functional Programming+-}++module 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)]]

@@ -0,0 +1,134 @@+{-# OPTIONS -Wall #-}++{- | +Module      :  LPFP.Geometry+Copyright   :  (c) Scott N. Walck 2023+License     :  BSD3 (see LICENSE)+Maintainer  :  Scott N. Walck <walck@lvc.edu>+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

@@ -0,0 +1,143 @@+{-# OPTIONS -Wall #-}++{- | +Module      :  LPFP.Integrals+Copyright   :  (c) Scott N. Walck 2023+License     :  BSD3 (see LICENSE)+Maintainer  :  Scott N. Walck <walck@lvc.edu>+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

@@ -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 <walck@lvc.edu>+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

@@ -0,0 +1,353 @@+{-# OPTIONS -Wall #-}++{- | +Module      :  LPFP.MOExamples+Copyright   :  (c) Scott N. Walck 2023+License     :  BSD3 (see LICENSE)+Maintainer  :  Scott N. Walck <walck@lvc.edu>+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)

@@ -0,0 +1,117 @@+{-# OPTIONS -Wall #-}++{- | +Module      :  LPFP.MagneticField+Copyright   :  (c) Scott N. Walck 2023+License     :  BSD3 (see LICENSE)+Maintainer  :  Scott N. Walck <walck@lvc.edu>+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)

@@ -0,0 +1,222 @@+{-# OPTIONS -Wall #-}++{- | +Module      :  LPFP.Maxwell+Copyright   :  (c) Scott N. Walck 2023+License     :  BSD3 (see LICENSE)+Maintainer  :  Scott N. Walck <walck@lvc.edu>+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

@@ -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 <walck@lvc.edu>+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)]])

@@ -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 <walck@lvc.edu>+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

@@ -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 <walck@lvc.edu>+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]

@@ -0,0 +1,230 @@+{-# OPTIONS -Wall #-}++{- | +Module      :  LPFP.Newton2+Copyright   :  (c) Scott N. Walck 2023+License     :  BSD3 (see LICENSE)+Maintainer  :  Scott N. Walck <walck@lvc.edu>+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

@@ -0,0 +1,256 @@+{-# OPTIONS -Wall #-}++{- | +Module      :  LPFP.SimpleVec+Copyright   :  (c) Scott N. Walck 2023+License     :  BSD3 (see LICENSE)+Maintainer  :  Scott N. Walck <walck@lvc.edu>+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)