diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,29 @@ +Copyright (c) 2022 Scott N. Walck . +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Scott N. Walck nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/LPFP-core.cabal b/LPFP-core.cabal new file mode 100644 --- /dev/null +++ b/LPFP-core.cabal @@ -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 diff --git a/src/LPFPCore.hs b/src/LPFPCore.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore.hs @@ -0,0 +1,419 @@ +{-# OPTIONS -Wall #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : LPFPCore +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from the book Learn Physics with Functional Programming +-} + +module 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 + ) diff --git a/src/LPFPCore/Charge.hs b/src/LPFPCore/Charge.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/Charge.hs @@ -0,0 +1,123 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFPCore.Charge +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 24 of the book Learn Physics with Functional Programming +-} + +module 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 diff --git a/src/LPFPCore/CoordinateSystems.hs b/src/LPFPCore/CoordinateSystems.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/CoordinateSystems.hs @@ -0,0 +1,174 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFPCore.CoordinateSystems +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 22 of the book Learn Physics with Functional Programming +-} + +module 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 diff --git a/src/LPFPCore/Current.hs b/src/LPFPCore/Current.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/Current.hs @@ -0,0 +1,102 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFPCore.Current +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 26 of the book Learn Physics with Functional Programming +-} + +module 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 diff --git a/src/LPFPCore/ElectricField.hs b/src/LPFPCore/ElectricField.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/ElectricField.hs @@ -0,0 +1,281 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFPCore.ElectricField +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 25 of the book Learn Physics with Functional Programming +-} + +module 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 diff --git a/src/LPFPCore/Electricity.hs b/src/LPFPCore/Electricity.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/Electricity.hs @@ -0,0 +1,77 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFPCore.Electricity +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 21 of the book Learn Physics with Functional Programming +-} + +module 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 diff --git a/src/LPFPCore/Geometry.hs b/src/LPFPCore/Geometry.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/Geometry.hs @@ -0,0 +1,134 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFPCore.Geometry +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 23 of the book Learn Physics with Functional Programming +-} + +module 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 diff --git a/src/LPFPCore/Integrals.hs b/src/LPFPCore/Integrals.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/Integrals.hs @@ -0,0 +1,143 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFPCore.Integrals +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code needed in chapter 24 of the book Learn Physics with Functional Programming +-} + +module 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 diff --git a/src/LPFPCore/Lorentz.hs b/src/LPFPCore/Lorentz.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/Lorentz.hs @@ -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 +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 diff --git a/src/LPFPCore/MOExamples.hs b/src/LPFPCore/MOExamples.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/MOExamples.hs @@ -0,0 +1,297 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFPCore.MOExamples +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 20 of the book Learn Physics with Functional Programming +-} + +module 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 diff --git a/src/LPFPCore/MagneticField.hs b/src/LPFPCore/MagneticField.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/MagneticField.hs @@ -0,0 +1,86 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFPCore.MagneticField +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 27 of the book Learn Physics with Functional Programming +-} + +module 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 diff --git a/src/LPFPCore/Maxwell.hs b/src/LPFPCore/Maxwell.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/Maxwell.hs @@ -0,0 +1,180 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFPCore.Maxwell +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 29 of the book Learn Physics with Functional Programming +-} + +module 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) diff --git a/src/LPFPCore/Mechanics1D.hs b/src/LPFPCore/Mechanics1D.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/Mechanics1D.hs @@ -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 +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) diff --git a/src/LPFPCore/Mechanics3D.hs b/src/LPFPCore/Mechanics3D.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/Mechanics3D.hs @@ -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 +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 diff --git a/src/LPFPCore/MultipleObjects.hs b/src/LPFPCore/MultipleObjects.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/MultipleObjects.hs @@ -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 +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] diff --git a/src/LPFPCore/Newton2.hs b/src/LPFPCore/Newton2.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/Newton2.hs @@ -0,0 +1,192 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFPCore.Newton2 +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 14 of the book Learn Physics with Functional Programming +-} + +module 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 diff --git a/src/LPFPCore/SimpleVec.hs b/src/LPFPCore/SimpleVec.hs new file mode 100644 --- /dev/null +++ b/src/LPFPCore/SimpleVec.hs @@ -0,0 +1,256 @@ +{-# OPTIONS -Wall #-} + +{- | +Module : LPFPCore.SimpleVec +Copyright : (c) Scott N. Walck 2023 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : stable + +Code from chapter 10 of the book Learn Physics with Functional Programming +-} + +module 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)