diff --git a/learn-physics.cabal b/learn-physics.cabal --- a/learn-physics.cabal +++ b/learn-physics.cabal @@ -1,5 +1,5 @@ Name: learn-physics -Version: 0.4.3 +Version: 0.5 Synopsis: Haskell code for learning physics Description: A library of functions for vector calculus, calculation of electric field, electric flux, @@ -11,8 +11,8 @@ Maintainer: Scott N. Walck Category: Physics Build-type: Simple -Cabal-version: >=1.6 -Tested-with: GHC == 7.6.3 +Cabal-version: >=1.8 +Tested-with: GHC == 7.8.2 Library Exposed-modules: Physics.Learn.Charge Physics.Learn.Current @@ -36,8 +36,8 @@ Physics.Learn.Visual.GlossTools Build-depends: base >= 4.2 && < 4.8, vector-space >= 0.8.4 && < 0.9, - not-gloss >= 0.5.0.4 && < 0.6, - spatial-math >= 0.1.7 && < 0.2, + not-gloss >= 0.6 && < 0.7, + spatial-math >= 0.2 && < 0.3, gloss >= 1.8 && < 1.9, gnuplot >= 0.5 && < 0.6 Hs-source-dirs: src diff --git a/src/Physics/Learn.hs b/src/Physics/Learn.hs --- a/src/Physics/Learn.hs +++ b/src/Physics/Learn.hs @@ -121,8 +121,6 @@ , simpleLineIntegral , dottedLineIntegral , crossedLineIntegral - , compositeSimpsonDottedLineIntegral - , compositeSimpsonCrossedLineIntegral -- ** Surfaces , Surface(..) , unitSphere @@ -162,9 +160,14 @@ , label , postscript , psFile + -- ** Gloss library + , polarToCart + , cartToPolar + , arrow + , thickArrow -- ** Vis library - , xyzFromVec - , xyzFromPos + , v3FromVec + , v3FromPos , visVec , oneVector , displayVectorField @@ -249,8 +252,6 @@ , simpleLineIntegral , dottedLineIntegral , crossedLineIntegral - , compositeSimpsonDottedLineIntegral - , compositeSimpsonCrossedLineIntegral ) import Physics.Learn.Surface ( Surface(..) @@ -275,8 +276,8 @@ , volumeIntegral ) import Physics.Learn.Visual.VisTools - ( xyzFromVec - , xyzFromPos + ( v3FromVec + , v3FromPos , visVec , oneVector , displayVectorField @@ -325,4 +326,10 @@ ( label , postscript , psFile + ) +import Physics.Learn.Visual.GlossTools + ( polarToCart + , cartToPolar + , arrow + , thickArrow ) diff --git a/src/Physics/Learn/Charge.hs b/src/Physics/Learn/Charge.hs --- a/src/Physics/Learn/Charge.hs +++ b/src/Physics/Learn/Charge.hs @@ -77,7 +77,7 @@ totalCharge :: ChargeDistribution -> Charge totalCharge (PointCharge q _) = q totalCharge (LineCharge lambda c) = simpleLineIntegral 1000 lambda c -totalCharge (SurfaceCharge sigma s) = surfaceIntegral 100 100 sigma s +totalCharge (SurfaceCharge sigma s) = surfaceIntegral 200 200 sigma s totalCharge (VolumeCharge rho v) = volumeIntegral 50 50 50 rho v totalCharge (MultipleCharges ds) = sum [totalCharge d | d <- ds] @@ -122,7 +122,7 @@ -> Surface -- ^ geometry of the surface charge -> VectorField -- ^ electric field (in V/m) eFieldFromSurfaceCharge sigma s r - = k *^ surfaceIntegral 100 100 integrand s + = k *^ surfaceIntegral 200 200 integrand s where k = 9e9 -- 1 / (4 * pi * epsilon0) integrand r' = sigma r' *^ d ^/ magnitude d ** 3 @@ -160,7 +160,7 @@ -- | The electric flux through a surface produced by a charge distribution. electricFlux :: Surface -> ChargeDistribution -> Double -electricFlux surf dist = dottedSurfaceIntegral 100 100 (eField dist) surf +electricFlux surf dist = dottedSurfaceIntegral 200 200 (eField dist) surf ------------------------ -- Electric Potential -- @@ -209,7 +209,7 @@ -> Surface -- ^ geometry of the surface charge -> ScalarField -- ^ electric potential ePotFromSurfaceCharge sigma s r - = k *^ surfaceIntegral 100 100 integrand s + = k *^ surfaceIntegral 200 200 integrand s where k = 9e9 -- 1 / (4 * pi * epsilon0) integrand r' = sigma r' / magnitude d diff --git a/src/Physics/Learn/Curve.hs b/src/Physics/Learn/Curve.hs --- a/src/Physics/Learn/Curve.hs +++ b/src/Physics/Learn/Curve.hs @@ -28,6 +28,8 @@ , simpleLineIntegral , dottedLineIntegral , crossedLineIntegral + , compositeTrapezoidDottedLineIntegral + , compositeTrapezoidCrossedLineIntegral , compositeSimpsonDottedLineIntegral , compositeSimpsonCrossedLineIntegral ) @@ -68,12 +70,32 @@ } -- | A dotted line integral. +-- Convenience function for 'compositeSimpsonDottedLineIntegral'. dottedLineIntegral + :: Int -- ^ number of half-intervals + -- (one less than the number of function evaluations) + -> VectorField -- ^ vector field + -> Curve -- ^ curve to integrate over + -> Double -- ^ scalar result +dottedLineIntegral = compositeSimpsonDottedLineIntegral + +-- | Calculates integral vf x dl over curve. +-- Convenience function for 'compositeSimpsonCrossedLineIntegral'. +crossedLineIntegral + :: Int -- ^ number of half-intervals + -- (one less than the number of function evaluations) + -> VectorField -- ^ vector field + -> Curve -- ^ curve to integrate over + -> Vec -- ^ vector result +crossedLineIntegral = compositeSimpsonCrossedLineIntegral + +-- | A dotted line integral, performed in an unsophisticated way. +compositeTrapezoidDottedLineIntegral :: Int -- ^ number of intervals -> VectorField -- ^ vector field -> Curve -- ^ curve to integrate over -> Double -- ^ scalar result -dottedLineIntegral n vf (Curve f a b) +compositeTrapezoidDottedLineIntegral n vf (Curve f a b) = sum $ zipWith (<.>) aveVecs dls where dt = (b - a) / fromIntegral n @@ -82,13 +104,13 @@ aveVecs = zipWith average vecs (tail vecs) dls = zipWith displacement pts (tail pts) --- | Calculates integral vf x dl over curve. -crossedLineIntegral +-- | Calculates integral vf x dl over curve in an unsophisticated way. +compositeTrapezoidCrossedLineIntegral :: Int -- ^ number of intervals -> VectorField -- ^ vector field -> Curve -- ^ curve to integrate over -> Vec -- ^ vector result -crossedLineIntegral n vf (Curve f a b) +compositeTrapezoidCrossedLineIntegral n vf (Curve f a b) = sumV $ zipWith (><) aveVecs dls where dt = (b - a) / fromIntegral n @@ -113,35 +135,6 @@ aveVecs = zipWith average vecs (tail vecs) dls = zipWith displacement pts (tail pts) -{- -lineIntegral :: (InnerSpace v, Scalar v ~ Double) => Double - -> (Vec -> v) - -> Curve - -> v -lineIntegral tol field (Curve f a b) - = let ca = f a - cb = f b - fielda = field ca - fieldb = field cb - val = average fielda fieldb ^* magnitude (cb ^-^ ca) - in evalInterval tol 1 20 field (Curve f a b) ca cb fielda fieldb val - -evalInterval :: (InnerSpace v, Scalar v ~ Double) => Double -> Int -> Int - -> (Vec -> v) -> Curve -> Vec -> Vec -> v -> v -> v -> v -evalInterval tol level maxlevel field (Curve f a b) ca cb fielda fieldb val - = let t = (a + b) / 2 - ct = f t - fieldt = field ct - vall = average fielda fieldt ^* magnitude (ct ^-^ ca) - valr = average fieldt fieldb ^* magnitude (cb ^-^ ct) - newval = vall ^+^ valr - in if magnitude (newval ^-^ val) < tol then - newval - else - evalInterval (tol/2) (level+1) maxlevel field (Curve f a t) ca ct fielda fieldt vall ^+^ - evalInterval (tol/2) (level+1) maxlevel field (Curve f t b) ct cb fieldt fieldb valr --} - -- | Reparametrize a curve from 0 to 1. normalizeCurve :: Curve -> Curve normalizeCurve (Curve f a b) @@ -227,10 +220,12 @@ -- Quadratic approximation to curve. -- Composite strategy. -- Dotted line integral. -compositeSimpsonDottedLineIntegral :: Int -- ^ number of half-intervals (one less than the number of function evaluations - -> VectorField -- ^ vector field - -> Curve -- ^ curve to integrate over - -> Double -- ^ scalar result +compositeSimpsonDottedLineIntegral + :: Int -- ^ number of half-intervals + -- (one less than the number of function evaluations) + -> VectorField -- ^ vector field + -> Curve -- ^ curve to integrate over + -> Double -- ^ scalar result compositeSimpsonDottedLineIntegral n vf (Curve c a b) = let nEven = 2 * div n 2 dt = (b - a) / fromIntegral nEven @@ -259,10 +254,12 @@ -- Quadratic approximation to curve. -- Composite strategy. -- Crossed line integral. -compositeSimpsonCrossedLineIntegral :: Int -- ^ number of half-intervals (one less than the number of function evaluations - -> VectorField -- ^ vector field - -> Curve -- ^ curve to integrate over - -> Vec -- ^ vector result +compositeSimpsonCrossedLineIntegral + :: Int -- ^ number of half-intervals + -- (one less than the number of function evaluations) + -> VectorField -- ^ vector field + -> Curve -- ^ curve to integrate over + -> Vec -- ^ vector result compositeSimpsonCrossedLineIntegral n vf (Curve c a b) = let nEven = 2 * div n 2 dt = (b - a) / fromIntegral nEven diff --git a/src/Physics/Learn/Visual/VisTools.hs b/src/Physics/Learn/Visual/VisTools.hs --- a/src/Physics/Learn/Visual/VisTools.hs +++ b/src/Physics/Learn/Visual/VisTools.hs @@ -3,8 +3,8 @@ -- | Some tools related to the not-gloss 3D graphics and animation library. module Physics.Learn.Visual.VisTools - ( xyzFromVec - , xyzFromPos + ( v3FromVec + , v3FromPos , visVec , oneVector , displayVectorField @@ -13,7 +13,7 @@ where import SpatialMath - ( Xyz(..) + ( V3(..) , Euler(..) ) import Vis @@ -37,17 +37,17 @@ ( Curve(..) ) --- | Make an 'Xyz' object from a 'Vec'. -xyzFromVec :: Vec -> Xyz Double -xyzFromVec v = Xyz x y z +-- | Make a 'V3' object from a 'Vec'. +v3FromVec :: Vec -> V3 Double +v3FromVec v = V3 x y z where x = xComp v y = yComp v z = zComp v --- | Make an 'Xyz' object from a 'Position'. -xyzFromPos :: Position -> Xyz Double -xyzFromPos r = Xyz x y z +-- | Make a 'V3' object from a 'Position'. +v3FromPos :: Position -> V3 Double +v3FromPos r = V3 x y z where (x,y,z) = cartesianCoordinates r @@ -58,16 +58,16 @@ -> VectorField -- ^ vector field to display -> VisObject Double -- ^ the displayable object displayVectorField col unitsPerMeter samplePts field - = VisObjects [Trans (xyzFromPos r) $ visVec col (e ^/ unitsPerMeter) | r <- samplePts, let e = field r] + = VisObjects [Trans (v3FromPos r) $ visVec col (e ^/ unitsPerMeter) | r <- samplePts, let e = field r] -- | A displayable VisObject for a curve. curveObject :: Color -> Curve -> VisObject Double curveObject color (Curve f a b) - = Line' [(xyzFromPos (f t), color) | t <- [a,a+(b-a)/1000..b]] + = Line' [(v3FromPos (f t), color) | t <- [a,a+(b-a)/1000..b]] -- | Place a vector at a particular position. oneVector :: Color -> Position -> Vec -> VisObject Double -oneVector c r v = Trans (xyzFromPos r) $ visVec c v +oneVector c r v = Trans (v3FromPos r) $ visVec c v data Cart = Cart Double Double Double deriving (Show) @@ -85,7 +85,7 @@ -- | A VisObject arrow from a vector visVec :: Color -> Vec -> VisObject Double -visVec c v = rotZ phi $ rotY theta $ Arrow (r,20*r) (Xyz 0 0 1) c +visVec c v = rotZ phi $ rotY theta $ Arrow (r,20*r) (V3 0 0 1) c where x = xComp v y = yComp v diff --git a/src/Tests.hs b/src/Tests.hs new file mode 100644 --- /dev/null +++ b/src/Tests.hs @@ -0,0 +1,55 @@ +{-# OPTIONS_GHC -Wall #-} + +module Main where + +import Physics.Learn +import Test.QuickCheck + +propGaussLaw1 :: (Double,Double,Double) -> Bool +propGaussLaw1 (x,y,z) = abs (eFlux - q/epsilon0) < 0.01 + where + eFlux = fluxThroughLargeCenteredSphere (x,y,z) q + epsilon0 = 1 / (4 * pi * 9e9) + q = epsilon0 + +fluxThroughLargeCenteredSphere :: (Double,Double,Double) -> Double -> Double +fluxThroughLargeCenteredSphere (x,y,z) q + = electricFlux (centeredSphere radius) (PointCharge q (cart x y z)) + where + radius = 2 * sqrt(x*x + y*y + z*z) + 1 + +currentLoop :: Double -> Current -> CurrentDistribution +currentLoop radius i + = LineCurrent i (Curve (\phi -> cyl radius phi 0) 0 (2*pi)) + +amperianLoop :: Double -> Curve +amperianLoop radius + = Curve (\t -> cart (radius + radius * sin t) 0 (radius * cos t)) 0 (2*pi) + +magCirculation :: Double -> Current -> Double +magCirculation radius i + = dottedLineIntegral 20 + (bFieldFromCurrentLoop i (Curve (\phi -> cyl radius phi 0) 0 (2*pi))) + (amperianLoop radius) + +bFieldFromCurrentLoop :: Current -> Curve -> VectorField +bFieldFromCurrentLoop i c r + = k *^ crossedLineIntegral 20 integrand c + where + k = 1e-7 -- mu0 / (4 * pi) + integrand r' = (-i) *^ d ^/ magnitude d ** 3 + where + d = displacement r' r + +propAmpere1 :: Double -> Property +propAmpere1 radius + = radius > 0 ==> abs (magCirculation radius i - 4*pi*1e-7 * i) < 0.01 + where + i = 1 / (4*pi*1e-7) + +main :: IO () +main = putStrLn "Gauss's law test:" >> + quickCheck propGaussLaw1 >> + putStrLn "Ampere's law test:" >> + quickCheck propAmpere1 +