LPFP-core (empty) → 1.1.1
raw patch · 19 files changed
+3457/−0 lines, 19 filesdep +basedep +containers
Dependencies added: base, containers
Files
- LICENSE +29/−0
- LPFP-core.cabal +40/−0
- src/LPFPCore.hs +419/−0
- src/LPFPCore/Charge.hs +123/−0
- src/LPFPCore/CoordinateSystems.hs +174/−0
- src/LPFPCore/Current.hs +102/−0
- src/LPFPCore/ElectricField.hs +281/−0
- src/LPFPCore/Electricity.hs +77/−0
- src/LPFPCore/Geometry.hs +134/−0
- src/LPFPCore/Integrals.hs +143/−0
- src/LPFPCore/Lorentz.hs +109/−0
- src/LPFPCore/MOExamples.hs +297/−0
- src/LPFPCore/MagneticField.hs +86/−0
- src/LPFPCore/Maxwell.hs +180/−0
- src/LPFPCore/Mechanics1D.hs +238/−0
- src/LPFPCore/Mechanics3D.hs +397/−0
- src/LPFPCore/MultipleObjects.hs +180/−0
- src/LPFPCore/Newton2.hs +192/−0
- src/LPFPCore/SimpleVec.hs +256/−0
+ LICENSE view
@@ -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.
+ LPFP-core.cabal view
@@ -0,0 +1,40 @@+name: LPFP-core+version: 1.1.1+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.+ LPFP-core elides all of the graphics dependencies of LPFP,+ so it has a much better chance of building without problems.+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+cabal-version: 1.12++library+ exposed-modules:+ LPFPCore+ LPFPCore.SimpleVec+ LPFPCore.Newton2+ LPFPCore.Mechanics1D+ LPFPCore.Mechanics3D+ LPFPCore.MultipleObjects+ LPFPCore.MOExamples+ LPFPCore.Electricity+ LPFPCore.CoordinateSystems+ LPFPCore.Geometry+ LPFPCore.Integrals+ LPFPCore.Charge+ LPFPCore.ElectricField+ LPFPCore.Current+ LPFPCore.MagneticField+ LPFPCore.Lorentz+ LPFPCore.Maxwell+ hs-source-dirs: src+ build-depends: base >= 4.7 && < 5,+ containers >= 0.6 && < 0.7+ default-language: Haskell2010
+ src/LPFPCore.hs view
@@ -0,0 +1,419 @@+{-# OPTIONS -Wall #-}+{-# LANGUAGE Trustworthy #-}++{- | +Module : LPFPCore+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 LPFPCore+ (+ -- * (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+ -- * 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+ -- * 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+ , sfTable+ -- * 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+ -- ** Maxwell Equations+ , directionalDerivative+ , curl+ , FieldState+ )+ where++import LPFPCore.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 LPFPCore.Newton2+ ( integral+ , antiDerivative+ )+import LPFPCore.Mechanics1D+ ( UpdateFunction+ , DifferentialEquation+ , NumericalMethod+ , RealVectorSpace(..)+ , Diff(..)+ , solver+ , euler+ , rungeKutta4+ )+import LPFPCore.Mechanics3D+ ( ParticleState(..)+ , DParticleState(..)+ , HasTime(..)+ , OneBodyForce+ , defaultParticleState+ , newtonSecondPS+ , relativityPS+ , earthSurfaceGravity+ , sunGravity+ , airResistance+ , windForce+ , uniformLorentzForce+ , eulerCromerPS+ , statesPS+ , updatePS+ )+import LPFPCore.MultipleObjects+ ( TwoBodyForce+ , Force(..)+ , MultiParticleState(..)+ , DMultiParticleState(..)+ , universalGravity+ , linearSpring+ , fixedLinearSpring+ , centralForce+ , billiardForce+ , newtonSecondMPS+ , eulerCromerMPS+ , updateMPS+ , statesMPS+ )+import LPFPCore.MOExamples+ ( Justification(..)+ , Table(..)+ , kineticEnergy+ , systemKE+ , momentum+ , systemP+ , linearSpringPE+ , earthSurfaceGravityPE+ , tenths+ , sigFigs+ )+import LPFPCore.Electricity+ ( Charge+ , elementaryCharge+ , coulombForce+ )+import LPFPCore.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+ , sfTable+ )+import LPFPCore.Geometry+ ( Curve(..)+ , unitCircle+ , straightLine+ , Surface(..)+ , unitSphere+ , centeredSphere+ , sphere+ , northernHemisphere+ , disk+ , shiftSurface+ , Volume(..)+ , unitBall+ , centeredBall+ , northernHalfBall+ , centeredCylinder+ )+import LPFPCore.Charge+ ( ChargeDistribution(..)+ , totalCharge+ , electricDipoleMoment+ )+import LPFPCore.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 LPFPCore.Current+ ( Current+ , CurrentDistribution(..)+ , crossedLineIntegral+ , totalCurrent+ , magneticDipoleMoment+ )+import LPFPCore.MagneticField+ ( bField+ )+import LPFPCore.Lorentz+ ( lorentzForce+ , newtonSecondPFS+ , defaultPFS+ )+import LPFPCore.Maxwell+ ( directionalDerivative+ , curl+ , FieldState+ )
+ src/LPFPCore/Charge.hs view
@@ -0,0 +1,123 @@+{-# OPTIONS -Wall #-}++{- | +Module : LPFPCore.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 LPFPCore.Charge where++import LPFPCore.SimpleVec ( R, Vec, vec, sumV, (*^), (^/), (<.>), magnitude, negateV )+import LPFPCore.Electricity ( elementaryCharge )+import LPFPCore.CoordinateSystems ( Position, ScalarField, origin, cart, sph+ , rVF, displacement, shiftPosition )+import LPFPCore.Geometry ( Curve(..), Surface(..), Volume(..)+ , straightLine, shiftSurface, disk )+import LPFPCore.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
+ src/LPFPCore/CoordinateSystems.hs view
@@ -0,0 +1,174 @@+{-# OPTIONS -Wall #-}++{- | +Module : LPFPCore.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 LPFPCore.CoordinateSystems where++import LPFPCore.SimpleVec+ ( R, Vec, (^/), vec, xComp, yComp, zComp, iHat, jHat, kHat+ , magnitude, sumV, zeroV )+import LPFPCore.MOExamples ( Table(..), Justification(..) )++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]++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]++magRad :: (R,R) -> (R,R)+magRad (x,y) = (sqrt (x*x + y*y), atan2 y x)++thetaSF :: ScalarField+thetaSF = undefined++thetaHat3D :: IO ()+thetaHat3D = undefined++phiHatGrad :: IO ()+phiHatGrad = undefined
+ src/LPFPCore/Current.hs view
@@ -0,0 +1,102 @@+{-# OPTIONS -Wall #-}++{- | +Module : LPFPCore.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 LPFPCore.Current where++import LPFPCore.SimpleVec+ ( R, Vec, sumV, (><), (*^) )+import LPFPCore.CoordinateSystems+ ( VectorField, rVF, cyl, phiHat )+import LPFPCore.Geometry+ ( Curve(..), Surface(..), Volume(..) )+import LPFPCore.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
+ src/LPFPCore/ElectricField.hs view
@@ -0,0 +1,281 @@+{-# OPTIONS -Wall #-}++{- | +Module : LPFPCore.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 LPFPCore.ElectricField where++import LPFPCore.SimpleVec+ ( R, Vec, (^+^), (^-^), (*^), (^*), (^/), (<.>), (><)+ , sumV, magnitude, vec )+import LPFPCore.CoordinateSystems+ ( Position, ScalarField, VectorField+ , displacement, shiftPosition, addVectorFields+ , origin, rVF )+import LPFPCore.Geometry ( Curve(..), Surface(..), Volume(..) )+import LPFPCore.Charge+ ( Charge, ChargeDistribution(..)+ , diskCap, 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++simpleDipoleSodiumChloride :: ChargeDistribution+simpleDipoleSodiumChloride = simpleDipole (vec 0 0 2.99e-29) 2.36e-10++eFieldSodiumChloride :: VectorField+eFieldSodiumChloride = eField simpleDipoleSodiumChloride++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)++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++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
+ src/LPFPCore/Electricity.hs view
@@ -0,0 +1,77 @@+{-# OPTIONS -Wall #-}++{- | +Module : LPFPCore.Electricity+Copyright : (c) Scott N. Walck 2023+License : BSD3 (see LICENSE)+Maintainer : Scott N. Walck <walck@lvc.edu>+Stability : stable++Code from chapter 21 of the book Learn Physics with Functional Programming+-}++module LPFPCore.Electricity where++import LPFPCore.SimpleVec+ ( Vec(..), R, (*^), iHat )+import LPFPCore.Mechanics3D+ ( ParticleState(..), defaultParticleState )+import LPFPCore.MultipleObjects+ ( TwoBodyForce, MultiParticleState(..), Force(..), statesMPS+ , eulerCromerMPS, centralForce )++type Charge = R++elementaryCharge :: Charge+elementaryCharge = 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++oneProtonPosition :: R -- dt+ -> R -- starting separation+ -> [(R,R)] -- (time,position) pairs+oneProtonPosition dt d+ = undefined dt d
+ src/LPFPCore/Geometry.hs view
@@ -0,0 +1,134 @@+{-# OPTIONS -Wall #-}++{- | +Module : LPFPCore.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 LPFPCore.Geometry where++import LPFPCore.SimpleVec ( R, Vec, (*^) )+import LPFPCore.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
+ src/LPFPCore/Integrals.hs view
@@ -0,0 +1,143 @@+{-# OPTIONS -Wall #-}++{- | +Module : LPFPCore.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 LPFPCore.Integrals where++import LPFPCore.SimpleVec+ ( R, Vec, (^+^), (*^), (^*), (^/), (<.>), (><), sumV, magnitude )+import LPFPCore.CoordinateSystems ( Position, ScalarField, VectorField+ , displacement, shiftPosition, origin, rVF )+import LPFPCore.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
+ src/LPFPCore/Lorentz.hs view
@@ -0,0 +1,109 @@+{-# OPTIONS -Wall #-}+{-# LANGUAGE MultiParamTypeClasses #-}++{- | +Module : LPFPCore.Lorentz+Copyright : (c) Scott N. Walck 2023+License : BSD3 (see LICENSE)+Maintainer : Scott N. Walck <walck@lvc.edu>+Stability : stable++Code from chapter 28 of the book Learn Physics with Functional Programming+-}++module LPFPCore.Lorentz where++import LPFPCore.SimpleVec ( R, Vec, (^+^), (*^), (^*), (^/), (><), zeroV )+import LPFPCore.Mechanics1D ( RealVectorSpace(..), Diff(..), rungeKutta4 )+import LPFPCore.Mechanics3D ( HasTime(..) )+import LPFPCore.CoordinateSystems ( Position(..), VectorField, origin+ , shiftPosition, addVectorFields )++data ParticleFieldState = ParticleFieldState { 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 }++scalePos :: R -> Position -> Position+scalePos metersPerVis (Cart x y z)+ = Cart (x/metersPerVis) (y/metersPerVis) (z/metersPerVis)++newtonSecondPFS' :: [ParticleFieldState -> Vec]+ -> ParticleFieldState -> DParticleFieldState+newtonSecondPFS' fs st = undefined fs st
+ src/LPFPCore/MOExamples.hs view
@@ -0,0 +1,297 @@+{-# OPTIONS -Wall #-}++{- | +Module : LPFPCore.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 LPFPCore.MOExamples where++import LPFPCore.SimpleVec+ ( R, Vec, (^+^), (^-^), (*^), vec, zeroV, magnitude+ , sumV, iHat, jHat, kHat, zComp )+import LPFPCore.Mechanics1D ( TimeStep, NumericalMethod, euler, rungeKutta4 )+import LPFPCore.Mechanics3D+ ( ParticleState(..), HasTime(..), defaultParticleState+ , earthSurfaceGravity )+import LPFPCore.MultipleObjects+ ( MultiParticleState(..), DMultiParticleState, Force(..), TwoBodyForce+ , newtonSecondMPS, updateMPS, statesMPS, eulerCromerMPS+ , linearSpring, fixedLinearSpring, billiardForce )++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++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]++-- 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
+ src/LPFPCore/MagneticField.hs view
@@ -0,0 +1,86 @@+{-# OPTIONS -Wall #-}++{- | +Module : LPFPCore.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 LPFPCore.MagneticField where++import LPFPCore.SimpleVec ( Vec(..), R+ , (^-^), (*^), (^/), (<.>), (><)+ , magnitude )+import LPFPCore.CoordinateSystems+ ( VectorField+ , rVF, displacement, addVectorFields )+import LPFPCore.Geometry ( Curve(..), Surface(..), Volume(..) )+import LPFPCore.ElectricField+ ( curveSample, surfaceSample, volumeSample+ , vectorSurfaceIntegral, vectorVolumeIntegral, mu0 )+import LPFPCore.Current+ ( Current, CurrentDistribution(..)+ , 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++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)++bFieldWireToroid :: VectorField+bFieldWireToroid = bField (wireToroid 0.3 1 50 10)++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
+ src/LPFPCore/Maxwell.hs view
@@ -0,0 +1,180 @@+{-# OPTIONS -Wall #-}++{- | +Module : LPFPCore.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 LPFPCore.Maxwell where++import LPFPCore.SimpleVec+ ( R, Vec(..), (^/), (^+^), (^-^), (*^)+ , vec, negateV, magnitude, xComp, yComp, zComp, iHat, jHat, kHat )+import LPFPCore.CoordinateSystems+ ( ScalarField, VectorField+ , cart, shiftPosition, rVF )+import LPFPCore.ElectricField ( cSI, mu0 )+import qualified Data.Map.Strict as M++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++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)
+ src/LPFPCore/Mechanics1D.hs view
@@ -0,0 +1,238 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}++{- | +Module : LPFPCore.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 LPFPCore.Mechanics1D where++import LPFPCore.Newton2 ( fAir )++import LPFPCore.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++pingpongPosition :: Time -> Velocity+pingpongPosition = positionFtxv 0.001 0.0027 (0,0.1,0) dampedHOForces++pingpongVelocity :: Time -> Velocity+pingpongVelocity = velocityFtxv 0.001 0.0027 (0,0.1,0) dampedHOForces++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)
+ src/LPFPCore/Mechanics3D.hs view
@@ -0,0 +1,397 @@+{-# OPTIONS -Wall #-}+{-# LANGUAGE MultiParamTypeClasses #-}++{- | +Module : LPFPCore.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 LPFPCore.Mechanics3D where++import LPFPCore.SimpleVec+ ( R, Vec, PosVec, (^+^), (^-^), (*^), (^*), (^/), (<.>), (><)+ , vec, sumV, magnitude, zeroV, xComp, yComp, zComp, iHat, jHat, kHat)+import LPFPCore.Mechanics1D+ ( RealVectorSpace(..), Diff(..), NumericalMethod+ , Time, TimeStep, rungeKutta4, solver )++-- | 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++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++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++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+ }++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++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+ }++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)++relativityPS' :: R -- c+ -> [OneBodyForce]+ -> ParticleState -> DParticleState+relativityPS' c fs st = undefined c fs st
+ src/LPFPCore/MultipleObjects.hs view
@@ -0,0 +1,180 @@+{-# OPTIONS -Wall #-}+{-# LANGUAGE MultiParamTypeClasses #-}++{- | +Module : LPFPCore.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 LPFPCore.MultipleObjects where++import LPFPCore.SimpleVec+ ( Vec, R, (^+^), (^-^), (*^), (^*), (^/), zeroV, magnitude )+import LPFPCore.Mechanics1D+ ( RealVectorSpace(..), Diff(..), NumericalMethod, Mass, TimeStep, euler )+import LPFPCore.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]
+ src/LPFPCore/Newton2.hs view
@@ -0,0 +1,192 @@+{-# OPTIONS -Wall #-}++{- | +Module : LPFPCore.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 LPFPCore.Newton2 where++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++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++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]++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]++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
+ src/LPFPCore/SimpleVec.hs view
@@ -0,0 +1,256 @@+{-# OPTIONS -Wall #-}++{- | +Module : LPFPCore.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 LPFPCore.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)