diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,29 @@ +Copyright (c) 2011-2014 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/learn-physics.cabal b/learn-physics.cabal new file mode 100644 --- /dev/null +++ b/learn-physics.cabal @@ -0,0 +1,34 @@ +Name: learn-physics +Version: 0.2 +Synopsis: Haskell code for learning physics +Description: A library of functions for vector calculus, + calculation of electric field, electric flux, + magnetic field, and other quantities in mechanics + and electromagnetic theory. +License: BSD3 +License-file: LICENSE +Author: Scott N. Walck +Maintainer: Scott N. Walck +Category: Physics +Build-type: Simple +Cabal-version: >=1.6 +Tested-with: GHC == 7.6.3 +Library + Exposed-modules: Physics.Learn.Charge + Physics.Learn.Current + Physics.Learn.Position + Physics.Learn.Curve + Physics.Learn.Surface + Physics.Learn.Volume + Physics.Learn.CarrotVec + Physics.Learn.SimpleVec + Physics.Learn.CommonVec + Physics.Learn.CoordinateFields + Physics.Learn.CoordinateSystem + Physics.Learn.StateSpace + Physics.Learn.RungeKutta + Physics.Learn.CompositeQuadrature + Physics.Learn.RootFinding + Build-depends: base >= 4.2 && < 4.8, + vector-space >= 0.8.4 && < 0.9 + Hs-source-dirs: src diff --git a/src/Physics/Learn/AdaptiveQuadrature.hs b/src/Physics/Learn/AdaptiveQuadrature.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/AdaptiveQuadrature.hs @@ -0,0 +1,294 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE TypeFamilies, FlexibleContexts #-} + +-- | Algorithm 4.2 of Burden and Faires, 5th edition + +module Physics.Learn.AdaptiveQuadrature +-- ( adaptiveQuad +-- ) + where + +import Data.VectorSpace + ( VectorSpace + , InnerSpace + , Scalar + , (^+^) + , (^-^) + , (*^) + , magnitude + , sumV + ) + +-- | Simplest, most elegant implementation. +-- Evaluates function at same spot multiple times. +adaptiveQuad :: Double -- ^ tolerance + -> Double -- ^ lower limit a + -> Double -- ^ upper limit b + -> (Double -> Double) -- ^ function f + -> Double -- ^ definite integral +adaptiveQuad tol a b f + = let s0 = simpson a b f + m = (a + b) / 2 + s1a = simpson a m f + s1b = simpson m b f + in if abs (s1a + s1b - s0) < 10 * tol + then s1a + s1b + else adaptiveQuad (tol/2) a m f + adaptiveQuad (tol/2) m b f + +simpson :: Double -- ^ lower limit a + -> Double -- ^ upper limit b + -> (Double -> Double) -- ^ function f + -> Double -- ^ Simpson approximation +simpson a b f = (b - a) / 6 * (f a + 4 * f ((a + b) / 2) + f b) + +-- | Version of adaptiveQuad for vectors. +-- Evaluates function at same spot multiple times. +adaptiveQuadVec :: (InnerSpace v, Scalar v ~ Double) => + Double -- ^ tolerance + -> Double -- ^ lower limit a + -> Double -- ^ upper limit b + -> (Double -> v) -- ^ function f + -> v -- ^ definite integral +adaptiveQuadVec tol a b f + = let s0 = simpsonVec a b f + m = (a + b) / 2 + s1a = simpsonVec a m f + s1b = simpsonVec m b f + in if magnitude (s1a ^+^ s1b ^-^ s0) < 10 * tol + then s1a ^+^ s1b + else adaptiveQuadVec (tol/2) a m f ^+^ adaptiveQuadVec (tol/2) m b f + +-- | Version of simpson for vectors. +simpsonVec :: (VectorSpace v, Scalar v ~ Double) => + Double -- ^ lower limit a + -> Double -- ^ upper limit b + -> (Double -> v) -- ^ function f + -> v -- ^ Simpson approximation +simpsonVec a b f = ((b - a) / 6) *^ (f a ^+^ 4 *^ f ((a + b) / 2) ^+^ f b) + +-- | Burden and Faires, Example 2 on page 197 +example2f :: Double -> Double +example2f x = (100 / x**2) * sin (10 / x) + +example2integral :: Double +example2integral = adaptiveQuad 1e-4 1 3 example2f + +-- *AdaptiveQuadrature> example2integral +-- -1.426014810049443 + +-- | Does no function evaluations itself. +simpleSimpson :: Double -- ^ lower limit a + -> Double -- ^ upper limit b + -> Double -- ^ value f(a) + -> Double -- ^ value f((a+b)/2) + -> Double -- ^ value f(b) + -> Double -- ^ Simpson approximation +simpleSimpson a b fa fm fb = (b - a) / 6 * (fa + 4 * fm + fb) + +-- The workhorse of the adaptive Simpson method. +-- Called by adaptiveSimpson +adaptiveSimpsonStep :: Double -- ^ tolerance + -> Double -- ^ lower limit a + -> Double -- ^ upper limit b + -> (Double -> Double) -- ^ function f + -> Double -- ^ value f(a) + -> Double -- ^ value f((a+b)/2) + -> Double -- ^ value f(b) + -> Double -- ^ definite integral +adaptiveSimpsonStep tol a b f fa fm fb + = let s0 = simpleSimpson a b fa fm fb + m = (a + b) / 2 + am = (a + m) / 2 + mb = (m + b) / 2 + fam = f am + fmb = f mb + s1a = simpleSimpson a m fa fam fm + s1b = simpleSimpson m b fm fmb fb + in if abs (s1a + s1b - s0) < 10 * tol + then s1a + s1b + else adaptiveSimpsonStep (tol/2) a m f fa fam fm + adaptiveSimpsonStep (tol/2) m b f fm fmb fb + +-- | This version is more efficient in that it does not +-- repeat function evaluations. +adaptiveSimpson :: Double -- ^ tolerance + -> Double -- ^ lower limit a + -> Double -- ^ upper limit b + -> (Double -> Double) -- ^ function f + -> Double -- ^ definite integral +adaptiveSimpson tol a b f + = let fa = f a + m = (a + b) / 2 + fm = f m + fb = f b + in adaptiveSimpsonStep tol a b f fa fm fb + +-- | Does no function evaluations itself. +-- For vector functions. +simpleSimpsonVec :: (VectorSpace v, Fractional (Scalar v)) => + Scalar v -- ^ lower limit a + -> Scalar v -- ^ upper limit b + -> v -- ^ value f(a) + -> v -- ^ value f((a+b)/2) + -> v -- ^ value f(b) + -> v -- ^ Simpson approximation +simpleSimpsonVec a b fa fm fb = ((b - a) / 6) *^ (fa ^+^ 4 *^ fm ^+^ fb) + +------------------------------------------ +-- Resource-limited adaptive quadrature -- +------------------------------------------ + +{- +Want a version that gives an error estimate, and can be used by +a scheduler for a resource-limited adaptive algorithm. +We won't achieve a desired precision, but rather we'll use +a fixed amount of resources in the best way possible. + +I think we'll need to create a data structure to hold the results +of evaluations so far so that they can be fed to the next step +if necessary. + +-- | This version does not repeat function evaluations. +-- It provides an error estimate. + + +-} + +-- data EvPair v = EvPair (Scalar v) v + +data SimpInterval3 v = SI3 { prLo :: (Scalar v, v) + , prMi :: (Scalar v, v) + , prHi :: (Scalar v, v) + , intEst3 :: v + } + +data SimpInterval5 v = SI5 { pr0 :: (Scalar v, v) + , pr1 :: (Scalar v, v) + , pr2 :: (Scalar v, v) + , pr3 :: (Scalar v, v) + , pr4 :: (Scalar v, v) + , intEst012 :: v + , intEst234 :: v + , intEst024 :: v + , integralEst :: v -- sum of intEst012 and intEst234 + , errorEst :: Scalar v + } + +divideInterval :: SimpInterval5 v -> (SimpInterval3 v, SimpInterval3 v) +divideInterval (SI5 xy0 xy1 xy2 xy3 xy4 ie012 ie234 _ie024 _ _) + = (SI3 xy0 xy1 xy2 ie012, SI3 xy2 xy3 xy4 ie234) + +refineInterval :: (InnerSpace v , Floating (Scalar v)) => + (Scalar v -> v) + -> SimpInterval3 v + -> SimpInterval5 v +refineInterval f (SI3 (x0,y0) (x2,y2) (x4,y4) ie024) + = let x1 = (x0 + x2) / 2 + x3 = (x2 + x4) / 2 + y1 = f x1 + y3 = f x3 + ie012 = simpleSimpsonVec x0 x2 y0 y1 y2 + ie234 = simpleSimpsonVec x2 x4 y2 y3 y4 + ie = ie012 ^+^ ie234 + errEst = 1/10 * magnitude (ie ^-^ ie024) -- 1/10 instead of 1/15 + in SI5 (x0,y0) (x1,y1) (x2,y2) (x3,y3) (x4,y4) ie012 ie234 ie024 ie errEst + +divideWorstInterval :: (InnerSpace v, Ord (Scalar v), Floating (Scalar v)) => + (Scalar v -> v) + -> [SimpInterval5 v] + -> [SimpInterval5 v] +divideWorstInterval _ [] = error "divideWorstInterval should never have been called on an empty list" +divideWorstInterval f (si:sis) + = let (si3a,si3b) = divideInterval si + si5a = refineInterval f si3a + si5b = refineInterval f si3b + in insertSorted si5a $ insertSorted si5b sis + +insertSorted :: Ord (Scalar v) => + SimpInterval5 v + -> [SimpInterval5 v] + -> [SimpInterval5 v] +insertSorted si5 [] = [si5] +insertSorted si5 (si:sis) = if errorEst si5 > errorEst si + then si5:si:sis + else si:insertSorted si5 sis + +adaptiveSimpEvalLimit :: (InnerSpace v, Ord (Scalar v), Floating (Scalar v)) => + Int -- ^ approximate number of function evals + -> Scalar v -- ^ lower limit + -> Scalar v -- ^ upper limit + -> (Scalar v -> v) -- ^ scalar or vector function + -> v -- ^ approximate integral +adaptiveSimpEvalLimit n a b f + = let m = (a + b) / 2 + fa = f a + fm = f m + fb = f b + ie = simpleSimpsonVec a b fa fm fb + si3 = SI3 (a,fa) (m,fm) (b,fb) ie + si5 = refineInterval f si3 + in sumV $ map integralEst $ last $ take (div n 4) $ iterate (divideWorstInterval f) [si5] + +{- +data SimpsonInterval5 v = SI5 { pLo :: Scalar v + , pHi :: Scalar v + , fLo :: v + , fLM :: v + , fM :: v + , fMH :: v + , fHi :: v + , integralEst :: v + , errorEst :: Scalar v + } +-} + +------------------------------- +-- Two-Dimensional integrals -- +------------------------------- + +adaptiveQuad2D :: Double -- ^ tolerance + -> Double -- ^ lower limit x_0 + -> Double -- ^ upper limit x_1 + -> (Double -> Double) -- ^ lower limit y_0(x) + -> (Double -> Double) -- ^ upper limit y_1(x) + -> (Double -> Double -> Double) -- ^ function f + -> Double -- ^ definite integral +adaptiveQuad2D tol x0 x1 y0 y1 f + = let f1 x = adaptiveQuad tol' (y0 x) (y1 x) (f x) + tol' = tol / abs (x1 - x0) + in adaptiveQuad tol x0 x1 f1 + +aq2dTest :: Double -> Double +aq2dTest tol = adaptiveQuad2D tol (-1) 1 (\y -> -sqrt(1 - y**2)) (\y -> sqrt(1-y**2)) (\_ _ -> 1) + +adaptiveSimpson2D :: Double -- ^ tolerance + -> Double -- ^ lower limit x_0 + -> Double -- ^ upper limit x_1 + -> (Double -> Double) -- ^ lower limit y_0(x) + -> (Double -> Double) -- ^ upper limit y_1(x) + -> (Double -> Double -> Double) -- ^ function f + -> Double -- ^ definite integral +adaptiveSimpson2D tol x0 x1 y0 y1 f + = let f1 x = adaptiveSimpson tol' (y0 x) (y1 x) (f x) + tol' = tol / abs (x1 - x0) + in adaptiveSimpson tol x0 x1 f1 + +adaptiveSimpson3D :: Double -- ^ tolerance + -> Double -- ^ lower limit x_0 + -> Double -- ^ upper limit x_1 + -> (Double -> Double) -- ^ lower limit y_0(x) + -> (Double -> Double) -- ^ upper limit y_1(x) + -> (Double -> Double -> Double) -- ^ lower limit z_0(x,y) + -> (Double -> Double -> Double) -- ^ upper limit z_1(x,y) + -> (Double -> Double -> Double -> Double) -- ^ function f + -> Double -- ^ definite integral +adaptiveSimpson3D tol x0 x1 y0 y1 z0 z1 f + = let f1 x = adaptiveSimpson2D tol' (y0 x) (y1 x) (z0 x) (z1 x) (f x) + tol' = tol / abs (x1 - x0) + in adaptiveSimpson tol x0 x1 f1 + +as3dTest :: Double -> Double +as3dTest tol = adaptiveSimpson3D tol (-1) 1 + (\y -> -sqrt(1 - y**2)) (\y -> sqrt(1-y**2)) + (\x y -> -sqrt(1 - x**2 - y**2)) (\x y -> sqrt(1 - x**2 - y**2)) + (\_ _ _ -> 1) + diff --git a/src/Physics/Learn/CarrotVec.hs b/src/Physics/Learn/CarrotVec.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/CarrotVec.hs @@ -0,0 +1,89 @@ +{-# LANGUAGE TypeFamilies #-} +{-# OPTIONS_GHC -Wall -fno-warn-orphans #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.CarrotVec +Copyright : (c) Scott N. Walck 2011-2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +This module defines some basic vector functionality. +It uses the same internal data representation as 'SimpleVec', +but declares 'Vec' to be an instance of 'VectorSpace'. +We import 'zeroV', 'negateV', 'sumV', '^+^', '^-^' +from 'AdditiveGroup', and +'*^', '^*', '^/', '<.>', 'magnitude' +from 'VectorSpace'. + +'CarrotVec' exports exactly the same symbols as 'SimpleVec'; +they are just defined differently. +-} + +-- 2011 Apr 10 +-- Definitions common to SimpleVec and CarrotVec have been put in CommonVec. + +module Physics.Learn.CarrotVec + ( Vec + , xComp + , yComp + , zComp + , vec + , (^+^) + , (^-^) + , (*^) + , (^*) + , (^/) + , (<.>) + , (><) + , magnitude + , zeroV + , negateV + , sumV + , iHat + , jHat + , kHat + ) + where + +import Data.VectorSpace + ( VectorSpace(..) + , InnerSpace(..) + , AdditiveGroup(..) + , Scalar + , (^+^) + , (^-^) + , (*^) + , (^*) + , (^/) + , (<.>) + , magnitude + , zeroV + , negateV + , sumV + ) +import Physics.Learn.CommonVec + ( Vec(..) + , xComp + , yComp + , zComp + , vec + , (><) + , iHat + , jHat + , kHat + ) + +instance AdditiveGroup Vec where + zeroV = vec 0 0 0 + negateV (Vec ax ay az) = Vec (-ax) (-ay) (-az) + Vec ax ay az ^+^ Vec bx by bz = Vec (ax+bx) (ay+by) (az+bz) + +instance VectorSpace Vec where + type Scalar Vec = Double + c *^ Vec ax ay az = Vec (c*ax) (c*ay) (c*az) + +instance InnerSpace Vec where + Vec ax ay az <.> Vec bx by bz = ax*bx + ay*by + az*bz + diff --git a/src/Physics/Learn/Charge.hs b/src/Physics/Learn/Charge.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/Charge.hs @@ -0,0 +1,234 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.Charge +Copyright : (c) Scott N. Walck 2011-2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +This module contains functions for working with charge, electric field, +electric flux, and electric potential. +-} + +module Physics.Learn.Charge + ( + -- * Charge + Charge + , ChargeDistribution(..) + -- * Electric Field + , eField + , eFieldFromPointCharge + , eFieldFromLineCharge + , eFieldFromSurfaceCharge + , eFieldFromVolumeCharge + -- * Electric Flux + , electricFlux + -- * Electric Potential + , electricPotentialFromField + , electricPotentialFromCharge + ) + where + +import Physics.Learn.CarrotVec + ( magnitude + , (*^) + , (^/) + ) +import Physics.Learn.Position + ( Position + , ScalarField + , VectorField + , displacement + , addFields + ) +import Physics.Learn.Curve + ( Curve(..) + , straightLine + , simpleLineIntegral + , dottedLineIntegral + ) +import Physics.Learn.Surface + ( Surface(..) + , surfaceIntegral + , dottedSurfaceIntegral + ) +import Physics.Learn.Volume + ( Volume(..) + , volumeIntegral + ) + +-- | 'Charge' is just a synonym for a double-precision floating point number. +type Charge = Double + +-- | A charge distribution is a point charge, a line charge, a surface charge, +-- a volume charge, or a combination of these. +-- The 'ScalarField' describes a linear charge density, a surface charge density, +-- or a volume charge density. +data ChargeDistribution = PointCharge Charge Position -- ^ point charge + | LineCharge ScalarField Curve -- ^ 'ScalarField' is linear charge density + | SurfaceCharge ScalarField Surface -- ^ 'ScalarField' is surface charge density + | VolumeCharge ScalarField Volume -- ^ 'ScalarField' is volume charge density + | Multiple [ChargeDistribution] -- ^ combination of charge distributions + +{- +shiftChargeDistribution :: Displacement -> ChargeDistribution -> ChargeDistribution +shiftChargeDistribution d (Point +-} + +-- | Electric field produced by a point charge. +-- The function 'eField' calls this function +-- to evaluate the electric field produced by a point charge. +eFieldFromPointCharge + :: Charge -- ^ charge (in Coulombs) + -> Position -- ^ of point charge + -> VectorField -- ^ electric field +eFieldFromPointCharge q r' r + = (k * q) *^ d ^/ magnitude d ** 3 + where + k = 9e9 -- 1 / (4 * pi * epsilon0) + d = displacement r' r + +-- | Electric field produced by a line charge. +-- The function 'eField' calls this function +-- to evaluate the electric field produced by a line charge. +eFieldFromLineCharge + :: ScalarField -- ^ linear charge density lambda + -> Curve -- ^ geometry of the line charge + -> VectorField -- ^ electric field +eFieldFromLineCharge lambda c r + = k *^ simpleLineIntegral 1000 integrand c + where + k = 9e9 -- 1 / (4 * pi * epsilon0) + integrand r' = lambda r' *^ d ^/ magnitude d ** 3 + where + d = displacement r' r + +-- | Electric field produced by a surface charge. +-- The function 'eField' calls this function +-- to evaluate the electric field produced by a surface charge. +eFieldFromSurfaceCharge + :: ScalarField -- ^ surface charge density sigma + -> Surface -- ^ geometry of the surface charge + -> VectorField -- ^ electric field +eFieldFromSurfaceCharge sigma s r + = k *^ surfaceIntegral 100 100 integrand s + where + k = 9e9 -- 1 / (4 * pi * epsilon0) + integrand r' = sigma r' *^ d ^/ magnitude d ** 3 + where + d = displacement r' r + +-- | Electric field produced by a volume charge. +-- The function 'eField' calls this function +-- to evaluate the electric field produced by a volume charge. +eFieldFromVolumeCharge + :: ScalarField -- ^ volume charge density rho + -> Volume -- ^ geometry of the volume charge + -> VectorField -- ^ electric field +eFieldFromVolumeCharge rho v r + = k *^ volumeIntegral 50 50 50 integrand v + where + k = 9e9 -- 1 / (4 * pi * epsilon0) + integrand r' = rho r' *^ d ^/ magnitude d ** 3 + where + d = displacement r' r + +-- | The electric field produced by a charge distribution. +-- This is the simplest way to find the electric field, because it +-- works for any charge distribution (point, line, surface, volume, or combination). +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 (Multiple cds) = addFields $ map eField cds + +------------------- +-- Electric Flux -- +------------------- + +-- | The electric flux through a surface produced by a charge distribution. +electricFlux :: Surface -> ChargeDistribution -> Double +electricFlux surf dist = dottedSurfaceIntegral 100 100 (eField dist) surf + +------------------------ +-- Electric Potential -- +------------------------ + +-- | Electric potential from electric field, given a position to be the zero +-- of electric potential. +electricPotentialFromField :: Position -- ^ position where electric potential is zero + -> VectorField -- ^ electric field + -> ScalarField -- ^ electric potential +electricPotentialFromField base ef r = -dottedLineIntegral 1000 ef (straightLine base r) + +-- | Electric potential produced by a charge distribution. +-- The position where the electric potential is zero is taken to be infinity. +electricPotentialFromCharge :: ChargeDistribution -> ScalarField +electricPotentialFromCharge (PointCharge q r') = ePotFromPointCharge q r' +electricPotentialFromCharge (LineCharge lam c) = ePotFromLineCharge lam c +electricPotentialFromCharge (SurfaceCharge sig s) = ePotFromSurfaceCharge sig s +electricPotentialFromCharge (VolumeCharge rho v) = ePotFromVolumeCharge rho v +electricPotentialFromCharge (Multiple cds) = addFields $ map electricPotentialFromCharge cds + +ePotFromPointCharge + :: Charge -- ^ charge (in Coulombs) + -> Position -- ^ of point charge + -> ScalarField -- ^ electric potential +ePotFromPointCharge q r' r + = (k * q) / magnitude d + where + k = 9e9 -- 1 / (4 * pi * epsilon0) + d = displacement r' r + +ePotFromLineCharge + :: ScalarField -- ^ linear charge density lambda + -> Curve -- ^ geometry of the line charge + -> ScalarField -- ^ electric potential +ePotFromLineCharge lambda c r + = k *^ simpleLineIntegral 1000 integrand c + where + k = 9e9 -- 1 / (4 * pi * epsilon0) + integrand r' = lambda r' / magnitude d + where + d = displacement r' r + +ePotFromSurfaceCharge + :: ScalarField -- ^ surface charge density sigma + -> Surface -- ^ geometry of the surface charge + -> ScalarField -- ^ electric potential +ePotFromSurfaceCharge sigma s r + = k *^ surfaceIntegral 100 100 integrand s + where + k = 9e9 -- 1 / (4 * pi * epsilon0) + integrand r' = sigma r' / magnitude d + where + d = displacement r' r + +ePotFromVolumeCharge + :: ScalarField -- ^ volume charge density rho + -> Volume -- ^ geometry of the volume charge + -> ScalarField -- ^ electric potential +ePotFromVolumeCharge rho v r + = k *^ volumeIntegral 50 50 50 integrand v + where + k = 9e9 -- 1 / (4 * pi * epsilon0) + integrand r' = rho r' / magnitude d + where + d = displacement r' r + +{- +Student Exercise: Write a function for electric potential difference. + +-- | The electric potential difference V(end) - V(beginning) between the endpoints +-- of a curve. +electricPotentialDifference :: Curve -> ChargeDistribution -> Double +electricPotentialDifference c dist = -dottedLineIntegral 1000 (eField dist) c +-} + +--------------------------------- +-- Common Charge Distributions -- +--------------------------------- + diff --git a/src/Physics/Learn/CommonVec.hs b/src/Physics/Learn/CommonVec.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/CommonVec.hs @@ -0,0 +1,69 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Safe #-} + +{- | +Module : Physics.Learn.CommonVec +Copyright : (c) Scott N. Walck 2012-2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +This module defines some common vector operations. +It is intended that this module not be imported directly, but that its +functionality be gained by importing either 'SimpleVec' or 'CarrotVec', +but not both. Choose 'SimpleVec' for vector operations +(such as vector addition) with simple concrete types, +which work only with the type 'Vec' of three-dimensional vectors. +Choose 'CarrotVec' for vector operations that work with any type in the +appropriate type class. +-} + +-- The definitions that are common to SimpleVec and CarrotVec. +-- We need to export the data constructor Vec for both SimpleVec and CarrotVec. + +module Physics.Learn.CommonVec + ( Vec(..) + , vec + , (><) + , iHat + , jHat + , kHat + ) + where + +infixl 7 >< + +-- | A type for vectors. +data Vec = Vec { xComp :: Double -- ^ x component + , yComp :: Double -- ^ y component + , zComp :: Double -- ^ z component + } deriving (Eq) + +instance Show Vec where + show (Vec x y z) = "vec " ++ showDouble x ++ " " + ++ showDouble y ++ " " + ++ showDouble z + +showDouble :: Double -> String +showDouble x + | x < 0 = "(" ++ show x ++ ")" + | otherwise = show x + +-- | Form a vector by giving its x, y, and z components. +vec :: Double -- ^ x component + -> Double -- ^ y component + -> Double -- ^ z component + -> Vec +vec = Vec + +-- | Cross product. +(><) :: Vec -> Vec -> Vec +Vec ax ay az >< Vec bx by bz = Vec (ay*bz - az*by) (az*bx - ax*bz) (ax*by - ay*bx) + +iHat, jHat, kHat :: Vec +-- | Unit vector in the x direction. +iHat = vec 1 0 0 +-- | Unit vector in the y direction. +jHat = vec 0 1 0 +-- | Unit vector in the z direction. +kHat = vec 0 0 1 diff --git a/src/Physics/Learn/CompositeQuadrature.hs b/src/Physics/Learn/CompositeQuadrature.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/CompositeQuadrature.hs @@ -0,0 +1,61 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.CompositeQuadrature +Copyright : (c) Scott N. Walck 2012-2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +Composite Trapezoid Rule and Composite Simpson's Rule +-} + +module Physics.Learn.CompositeQuadrature + ( compositeTrapezoid + , compositeSimpson + ) + where + +import Data.VectorSpace + ( VectorSpace + , Scalar + , (^+^) + , (*^) + , zeroV + ) + +-- | Composite Trapezoid Rule +compositeTrapezoid :: (VectorSpace v, Fractional (Scalar v)) => + Int -- ^ number of intervals (one less than the number of function evaluations) + -> Scalar v -- ^ lower limit + -> Scalar v -- ^ upper limit + -> (Scalar v -> v) -- ^ function to be integrated + -> v -- ^ definite integral +compositeTrapezoid n a b f + = let dt = (b - a) / fromIntegral n + ts = [a + fromIntegral m * dt | m <- [0..n]] + pairs = [(t,f t) | t <- ts] + combine [] = error "compositeSimpson: odd number of half-intervals" -- this should never happen + combine [_] = zeroV + combine ((t0,f0):(t1,f1):ps) = ((t1 - t0) / 2) *^ (f0 ^+^ f1) ^+^ combine ((t1,f1):ps) + in combine pairs + +-- | Composite Simpson's Rule +compositeSimpson :: (VectorSpace v, Fractional (Scalar v)) => + Int -- ^ number of half-intervals (one less than the number of function evaluations) + -> Scalar v -- ^ lower limit + -> Scalar v -- ^ upper limit + -> (Scalar v -> v) -- ^ function to be integrated + -> v -- ^ definite integral +compositeSimpson n a b f + = let nEven = 2 * div n 2 + dt = (b - a) / fromIntegral nEven + ts = [a + fromIntegral m * dt | m <- [0..nEven]] + pairs = [(t,f t) | t <- ts] + combine [] = error "compositeSimpson: odd number of half-intervals" -- this should never happen + combine [_] = zeroV + combine (_:_:[]) = error "compositeSimpson: odd number of half-intervals" -- this should never happen + combine ((t0,f0):(_,f1):(t2,f2):ps) = ((t2 - t0) / 6) *^ (f0 ^+^ 4 *^ f1 ^+^ f2) ^+^ combine ((t2,f2):ps) + in combine pairs diff --git a/src/Physics/Learn/CoordinateFields.hs b/src/Physics/Learn/CoordinateFields.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/CoordinateFields.hs @@ -0,0 +1,73 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.CoordinateFields +Copyright : (c) Scott N. Walck 2012-2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +Coordinate fields for Cartesian, cylindrical, and spherical coordinates. +-} + +module Physics.Learn.CoordinateFields + ( x + , y + , z + , s + , phi + , r + , theta + ) + where + +import Physics.Learn.Position + ( ScalarField + , cartesianCoordinates + , cylindricalCoordinates + , sphericalCoordinates + ) + +fst3 :: (a,b,c) -> a +fst3 (v,_,_) = v + +snd3 :: (a,b,c) -> b +snd3 (_,v,_) = v + +thd3 :: (a,b,c) -> c +thd3 (_,_,v) = v + +-- | The x Cartesian coordinate of a position. +x :: ScalarField +x = fst3 . cartesianCoordinates + +-- | The y Cartesian coordinate of a position. +y :: ScalarField +y = snd3 . cartesianCoordinates + +-- | The z Cartesian (or cylindrical) coordinate of a position. +z :: ScalarField +z = thd3 . cartesianCoordinates + +-- | The s cylindrical coordinate of a position. +-- This is the distance of the position from the z axis. +s :: ScalarField +s = fst3 . cylindricalCoordinates + +-- | The phi cylindrical (or spherical) coordinate of a position. +-- This is the angle from the positive x axis +-- to the projection of the position onto the xy plane. +phi :: ScalarField +phi = snd3 . cylindricalCoordinates + +-- | The r spherical coordinate of a position. +-- This is the distance of the position from the origin. +r :: ScalarField +r = fst3 . sphericalCoordinates + +-- | The theta spherical coordinate of a position. +-- This is the angle from the positive z axis to the position. +theta :: ScalarField +theta = snd3 . sphericalCoordinates + diff --git a/src/Physics/Learn/CoordinateSystem.hs b/src/Physics/Learn/CoordinateSystem.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/CoordinateSystem.hs @@ -0,0 +1,61 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.CoordinateSystem +Copyright : (c) Scott N. Walck 2012-2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +A module for working with coordinate systems. +-} + +module Physics.Learn.CoordinateSystem + ( CoordinateSystem(..) + , standardCartesian + , standardCylindrical + , standardSpherical + , newCoordinateSystem + ) + where + +import Physics.Learn.Position + ( Position + , cartesian + , cartesianCoordinates + , cylindrical + , cylindricalCoordinates + , spherical + , sphericalCoordinates + ) + +-- | Specification of a coordinate system requires +-- a map from coordinates into space, and +-- a map from space into coordinates. +data CoordinateSystem + = CoordinateSystem { toPosition :: (Double,Double,Double) -> Position -- ^ a map from coordinates into space + , fromPosition :: Position -> (Double,Double,Double) -- ^ a map from space into coordinates + } + +-- | The standard Cartesian coordinate system +standardCartesian :: CoordinateSystem +standardCartesian = CoordinateSystem cartesian cartesianCoordinates + +-- | The standard cylindrical coordinate system +standardCylindrical :: CoordinateSystem +standardCylindrical = CoordinateSystem cylindrical cylindricalCoordinates + +-- | The standard spherical coordinate system +standardSpherical :: CoordinateSystem +standardSpherical = CoordinateSystem spherical sphericalCoordinates + +-- | Define a new coordinate system in terms of an existing one. +-- First parameter is a map from old coordinates to new coordinates. +-- Second parameter is the inverse map from new coordinates to old coordinates. +newCoordinateSystem :: ((Double,Double,Double) -> (Double,Double,Double)) -- ^ (x',y',z') = f(x,y,z) + -> ((Double,Double,Double) -> (Double,Double,Double)) -- ^ (x,y,z) = g(x',y',z') + -> CoordinateSystem -- ^ old coordinate system + -> CoordinateSystem +newCoordinateSystem f g (CoordinateSystem tp fp) + = CoordinateSystem (tp . g) (f . fp) diff --git a/src/Physics/Learn/Current.hs b/src/Physics/Learn/Current.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/Current.hs @@ -0,0 +1,131 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.Current +Copyright : (c) Scott N. Walck 2012-2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +This module contains functions for working with current, magnetic field, +and magnetic flux. +-} + +module Physics.Learn.Current + ( + -- * Current + Current + , CurrentDistribution(..) + -- * Magnetic Field + , bField + , bFieldFromLineCurrent + , bFieldFromSurfaceCurrent + , bFieldFromVolumeCurrent + -- * Magnetic Flux + , magneticFlux + ) + where + +import Physics.Learn.CarrotVec + ( magnitude + , (*^) + , (^/) + , (><) + ) +import Physics.Learn.Position + ( VectorField + , displacement + , addFields + ) +import Physics.Learn.Curve + ( Curve(..) + , crossedLineIntegral + ) +import Physics.Learn.Surface + ( Surface(..) + , surfaceIntegral + , dottedSurfaceIntegral + ) +import Physics.Learn.Volume + ( Volume(..) + , volumeIntegral + ) + +-- | 'Current' is just a synonym for a double-precision floating point number. +type Current = Double + +-- | A current distribution is a line current (current through a wire), a surface current, +-- a volume current, or a combination of these. +-- The 'VectorField' describes a surface current density +-- or a volume current density. +data CurrentDistribution = LineCurrent Current Curve -- ^ current through a wire + | SurfaceCurrent VectorField Surface -- ^ 'VectorField' is surface current density + | VolumeCurrent VectorField Volume -- ^ 'VectorField' is volume current density + | MultipleCurrents [CurrentDistribution] -- ^ combination of current distributions + +-- | Magnetic field produced by a line current (current through a wire). +-- The function 'bField' calls this function +-- to evaluate the magnetic field produced by a line current. +bFieldFromLineCurrent + :: Current -- ^ current (in Amps) + -> Curve -- ^ geometry of the line current + -> VectorField -- ^ magnetic field +bFieldFromLineCurrent i c r + = k *^ crossedLineIntegral 1000 integrand c + where + k = 1e-7 -- mu0 / (4 * pi) + integrand r' = (-i) *^ d ^/ magnitude d ** 3 + where + d = displacement r' r + +-- | Magnetic field produced by a surface current. +-- The function 'bField' calls this function +-- to evaluate the magnetic field produced by a surface current. +-- This function assumes that surface current density +-- will be specified parallel to the surface, and does +-- not check if that is true. +bFieldFromSurfaceCurrent + :: VectorField -- ^ surface current density + -> Surface -- ^ geometry of the surface current + -> VectorField -- ^ magnetic field +bFieldFromSurfaceCurrent kCurrent c r + = k *^ surfaceIntegral 100 100 integrand c + where + k = 1e-7 -- mu0 / (4 * pi) + integrand r' = (kCurrent r' >< d) ^/ magnitude d ** 3 + where + d = displacement r' r + +-- | Magnetic field produced by a volume current. +-- The function 'bField' calls this function +-- to evaluate the magnetic field produced by a volume current. +bFieldFromVolumeCurrent + :: VectorField -- ^ volume current density + -> Volume -- ^ geometry of the volume current + -> VectorField -- ^ magnetic field +bFieldFromVolumeCurrent j c r + = k *^ volumeIntegral 50 50 50 integrand c + where + k = 1e-7 -- mu0 / (4 * pi) + integrand r' = (j r' >< d) ^/ magnitude d ** 3 + where + d = displacement r' r + +-- | The magnetic field produced by a current distribution. +-- This is the simplest way to find the magnetic field, because it +-- works for any current distribution (line, surface, volume, or combination). +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) = addFields $ map bField cds + +------------------- +-- Magnetic Flux -- +------------------- + +-- | The magnetic flux through a surface produced by a current distribution. +magneticFlux :: Surface -> CurrentDistribution -> Double +magneticFlux surf dist = dottedSurfaceIntegral 100 100 (bField dist) surf + diff --git a/src/Physics/Learn/Curve.hs b/src/Physics/Learn/Curve.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/Curve.hs @@ -0,0 +1,278 @@ +{-# LANGUAGE TypeFamilies, FlexibleContexts #-} +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.Curve +Copyright : (c) Scott N. Walck 2012-2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +This module contains functions for working with 'Curve's +and line integrals along 'Curve's. +-} + +module Physics.Learn.Curve + ( + -- * Curves + Curve(..) + , normalizeCurve + , concatCurves + , concatenateCurves + , reverseCurve + , evalCurve + , shiftCurve + , straightLine + -- * Line Integrals + , simpleLineIntegral + , dottedLineIntegral + , crossedLineIntegral + , compositeSimpsonDottedLineIntegral + , compositeSimpsonCrossedLineIntegral + ) + where + +import Data.VectorSpace + ( VectorSpace + , InnerSpace + , Scalar + ) +import Physics.Learn.CarrotVec + ( Vec + , (><) + , (<.>) + , sumV + , (^*) + , (^/) + , (^+^) + , (^-^) + , (*^) + , magnitude + , zeroV + , negateV + ) +import Physics.Learn.Position + ( Position + , Displacement + , displacement + , Field + , VectorField + , shiftPosition + ) + +-- | 'Curve' is a parametrized function into three-space, an initial limit, and a final limit. +data Curve = Curve { curveFunc :: (Double -> Position) -- ^ function from one parameter into space + , startingCurveParam :: Double -- ^ starting value of the parameter + , endingCurveParam :: Double -- ^ ending value of the parameter + } + +-- | A dotted line integral. +dottedLineIntegral + :: Int -- ^ number of intervals + -> VectorField -- ^ vector field + -> Curve -- ^ curve to integrate over + -> Double -- ^ scalar result +dottedLineIntegral n vf (Curve f a b) + = sum $ zipWith (<.>) aveVecs dls + where + dt = (b - a) / fromIntegral n + pts = [f t | t <- [a,a+dt..b]] + vecs = [vf pt | pt <- pts] + aveVecs = zipWith average vecs (tail vecs) + dls = zipWith displacement pts (tail pts) + +-- | Calculates integral vf x dl over curve. +crossedLineIntegral + :: Int -- ^ number of intervals + -> VectorField -- ^ vector field + -> Curve -- ^ curve to integrate over + -> Vec -- ^ vector result +crossedLineIntegral n vf (Curve f a b) + = sumV $ zipWith (><) aveVecs dls + where + dt = (b - a) / fromIntegral n + pts = [f t | t <- [a,a+dt..b]] + vecs = [vf pt | pt <- pts] + aveVecs = zipWith average vecs (tail vecs) + dls = zipWith displacement pts (tail pts) + +-- | Calculates integral f dl over curve, where dl is a scalar line element. +simpleLineIntegral + :: (InnerSpace v, Scalar v ~ Double) + => Int -- ^ number of intervals + -> Field v -- ^ scalar or vector field + -> Curve -- ^ curve to integrate over + -> v -- ^ scalar or vector result +simpleLineIntegral n vf (Curve f a b) + = sumV $ zipWith (^*) aveVecs (map magnitude dls) + where + dt = (b - a) / fromIntegral n + pts = [f t | t <- [a,a+dt..b]] + vecs = [vf pt | pt <- pts] + 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) + = Curve (f . scl) 0 1 + where + scl t = a + (b - a) * t + +-- | Concatenate two curves. +concatCurves :: Curve -- ^ go first along this curve + -> Curve -- ^ then along this curve + -> Curve -- ^ to produce this new curve +concatCurves c1 c2 + = normalizeCurve $ Curve f 0 2 + where + (Curve f1 _ _) = normalizeCurve c1 + (Curve f2 _ _) = normalizeCurve c2 + f t | t <= 1 = f1 t + | otherwise = f2 (t-1) + +-- | Concatenate a list of curves. +-- Parametrizes curves equally. +concatenateCurves :: [Curve] -> Curve +concatenateCurves [] = error "concatenateCurves: cannot concatenate empty list" +concatenateCurves cs = normalizeCurve $ Curve f 0 (fromIntegral n) + where + n = length cs + ncs = map normalizeCurve cs + f t = evalCurve (ncs !! m) (t - fromIntegral m) + where m = min (n-1) (floor t) + +-- | Reverse a curve. +reverseCurve :: Curve -> Curve +reverseCurve (Curve f a b) + = Curve (f . rev) a b + where + rev t = a + b - t + +-- | Evaluate the position of a curve at a parameter. +evalCurve :: Curve -- ^ the curve + -> Double -- ^ the parameter + -> Position -- ^ position of the point on the curve at that parameter +evalCurve (Curve f _ _) t = f t + +-- | Shift a curve by a displacement. +shiftCurve :: Displacement -- ^ amount to shift + -> Curve -- ^ original curve + -> Curve -- ^ shifted curve +shiftCurve d (Curve f sl su) + = Curve (shiftPosition d . f) sl su + +-- | The straight-line curve from one position to another. +straightLine :: Position -- ^ starting position + -> Position -- ^ ending position + -> Curve -- ^ straight-line curve +straightLine r1 r2 = Curve f 0 1 + where + f t = shiftPosition (t *^ d) r1 + d = displacement r1 r2 + +------------- +-- Helpers -- +------------- + +average :: (VectorSpace v, Scalar v ~ Double) => v -> v -> v +average v1 v2 = (v1 ^+^ v2) ^/ 2 + +---------------------------------------- +-- Quadratic (Simpson) Approximations -- +---------------------------------------- + +dottedSimp :: (InnerSpace v, Fractional (Scalar v)) => + v -- ^ vector field low + -> v -- ^ vector field mid + -> v -- ^ vector field high + -> v -- ^ dl low to mid + -> v -- ^ dl mid to high + -> Scalar v -- ^ quadratic approximation +dottedSimp f0 f1 f2 g10 g21 + = ((g21 ^+^ g10) ^/ 6) <.> (f0 ^+^ 4 *^ f1 ^+^ f2) + + ((g21 ^-^ g10) ^/ 3) <.> (f2 ^-^ f0) + +-- | Quadratic approximation to vector field. +-- 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 n vf (Curve c a b) + = let nEven = 2 * div n 2 + dt = (b - a) / fromIntegral nEven + ts = [a + fromIntegral m * dt | m <- [0..nEven]] + pairs = [(ct,vf ct) | t <- ts, let ct = c t] + combine [] = error "compositeSimpson: odd number of half-intervals" -- this should never happen + combine [_] = zeroV + combine (_:_:[]) = error "compositeSimpson: odd number of half-intervals" -- this should never happen + combine ((c0,f0):(c1,f1):(c2,f2):ps) + = dottedSimp f0 f1 f2 (displacement c0 c1) (displacement c1 c2) + ^+^ combine ((c2,f2):ps) + in combine pairs + +crossedSimp :: Vec -- ^ vector field low + -> Vec -- ^ vector field mid + -> Vec -- ^ vector field high + -> Vec -- ^ dl low to mid + -> Vec -- ^ dl mid to high + -> Vec -- ^ quadratic approximation +crossedSimp f0 f1 f2 g10 g21 + = negateV $ + ((g21 ^+^ g10) ^/ 6) >< (f0 ^+^ 4 *^ f1 ^+^ f2) + ^+^ ((g21 ^-^ g10) ^/ 3) >< (f2 ^-^ f0) + +-- | Quadratic approximation to vector field. +-- 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 n vf (Curve c a b) + = let nEven = 2 * div n 2 + dt = (b - a) / fromIntegral nEven + ts = [a + fromIntegral m * dt | m <- [0..nEven]] + pairs = [(ct,vf ct) | t <- ts, let ct = c t] + combine [] = error "compositeSimpson: odd number of half-intervals" -- this should never happen + combine [_] = zeroV + combine (_:_:[]) = error "compositeSimpson: odd number of half-intervals" -- this should never happen + combine ((c0,f0):(c1,f1):(c2,f2):ps) + = crossedSimp f0 f1 f2 (displacement c0 c1) (displacement c1 c2) + ^+^ combine ((c2,f2):ps) + in combine pairs + diff --git a/src/Physics/Learn/Position.hs b/src/Physics/Learn/Position.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/Position.hs @@ -0,0 +1,270 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE TypeSynonymInstances, FlexibleInstances #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.Position +Copyright : (c) Scott N. Walck 2012-2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +A module for working with the idea of position and coordinate systems. +-} + +module Physics.Learn.Position + ( Position + , Displacement + , ScalarField + , VectorField + , Field + , CoordinateSystem + , cartesian + , cylindrical + , spherical + , cart + , cyl + , sph + , cartesianCoordinates + , cylindricalCoordinates + , sphericalCoordinates + , displacement + , shiftPosition + , shiftObject + , shiftField + , addFields + , rHat + , thetaHat + , phiHat + , sHat + , xHat + , yHat + , zHat + ) + where + +import Data.VectorSpace + ( AdditiveGroup + ) +import Physics.Learn.CarrotVec + ( Vec + , vec + , xComp + , yComp + , zComp + , iHat + , jHat + , kHat + , sumV + , magnitude + , (^/) + ) + +-- | A type for position. +-- Position is not a vector because it makes no sense to add positions. +data Position = Cart Double Double Double + +-- | A displacement is a vector. +type Displacement = Vec + +-- | A scalar field associates a number with each position in space. +type ScalarField = Position -> Double + +{- +-- | Scalar fields can be added, subtracted, multiplied, and negated, +-- just like scalars themselves. +instance Num ScalarField where + (f + g) x = f x + g x + (f * g) x = f x * g x + (f - g) x = f x - g x + negate f x = negate (f x) + abs f x = abs (f x) + signum f x = signum (f x) + fromInteger n = const (fromInteger n) + +-- | Scalar fields can be divided, just like scalars themselves. +instance Fractional ScalarField where + (f / g) x = f x / g x + recip f x = recip (f x) + fromRational rat = const (fromRational rat) + +-- | Cosine of a scalar field, etc. +instance Floating ScalarField where + pi = const pi + exp f x = exp (f x) + sqrt f x = sqrt (f x) + log f x = log (f x) + (f ** g) x = f x ** g x + logBase f g x = logBase (f x) (g x) + sin f x = sin (f x) + cos f x = cos (f x) + tan f x = tan (f x) + asin f x = asin (f x) + acos f x = acos (f x) + atan f x = atan (f x) + sinh f x = sinh (f x) + cosh f x = cosh (f x) + tanh f x = tanh (f x) + asinh f x = asinh (f x) + acosh f x = acosh (f x) + atanh f x = atanh (f x) +-} + +-- | A vector field associates a vector with each position in space. +type VectorField = Position -> Vec + +-- | Sometimes we want to be able to talk about a field without saying +-- whether it is a scalar field or a vector field. +type Field v = Position -> v + +-- | A coordinate system is a function from three parameters to space. +type CoordinateSystem = (Double,Double,Double) -> Position + +-- | Add two scalar fields or two vector fields. +addFields :: AdditiveGroup v => [Field v] -> Field v +addFields flds r = sumV [fld r | fld <- flds] + +-- | The Cartesian coordinate system. Coordinates are (x,y,z). +cartesian :: CoordinateSystem +cartesian (x,y,z) = Cart x y z + +-- | The cylindrical coordinate system. Coordinates are (s,phi,z), +-- where s is the distance from the z axis and phi is the angle +-- with the x axis. +cylindrical :: CoordinateSystem +cylindrical (s,phi,z) = Cart (s * cos phi) (s * sin phi) z + +-- | The spherical coordinate system. Coordinates are (r,theta,phi), +-- where r is the distance from the origin, theta is the angle with +-- the z axis, and phi is the azimuthal angle. +spherical :: CoordinateSystem +spherical (r,th,phi) = Cart (r * sin th * cos phi) (r * sin th * sin phi) (r * cos th) + +-- | A helping function to take three numbers x, y, and z and form the +-- appropriate position using Cartesian coordinates. +cart :: Double -- ^ x coordinate + -> Double -- ^ y coordinate + -> Double -- ^ z coordinate + -> Position +cart = Cart + +-- | A helping function to take three numbers s, phi, and z and form the +-- appropriate position using cylindrical coordinates. +cyl :: Double -- ^ s coordinate + -> Double -- ^ phi coordinate + -> Double -- ^ z coordinate + -> Position +cyl s phi z = Cart (s * cos phi) (s * sin phi) z + +-- | A helping function to take three numbers r, theta, and phi and form the +-- appropriate position using spherical coordinates. +sph :: Double -- ^ r coordinate + -> Double -- ^ theta coordinate + -> Double -- ^ phi coordinate + -> Position +sph r theta phi = Cart (r * sin theta * cos phi) (r * sin theta * sin phi) (r * cos theta) + +-- | Returns the three Cartesian coordinates as a triple from a position. +cartesianCoordinates :: Position -> (Double,Double,Double) +cartesianCoordinates (Cart x y z) = (x,y,z) + +-- | Returns the three cylindrical coordinates as a triple from a position. +cylindricalCoordinates :: Position -> (Double,Double,Double) +cylindricalCoordinates (Cart x y z) = (s,phi,z) + where + s = sqrt(x**2 + y**2) + phi = atan2 y x + +-- | Returns the three spherical coordinates as a triple from a position. +sphericalCoordinates :: Position -> (Double,Double,Double) +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 + +-- | Displacement from source position to target position. +displacement :: Position -- ^ source position + -> Position -- ^ target position + -> Displacement +displacement (Cart x' y' z') (Cart x y z) = vec (x-x') (y-y') (z-z') + +-- | Shift a position by a displacement. +shiftPosition :: Displacement -> Position -> Position +shiftPosition v (Cart x y z) = Cart (x + xComp v) (y + yComp v) (z + zComp v) + +-- | An object is a map into 'Position'. +shiftObject :: Displacement -> (a -> Position) -> (a -> Position) +shiftObject d f = shiftPosition d . f + +-- | A field is a map from 'Position'. +shiftField :: Displacement -> (Position -> v) -> (Position -> v) +shiftField d f = f . shiftPosition d + +-- | The vector field in which each point in space is associated +-- with a unit vector in the direction of increasing spherical coordinate +-- r, while spherical coordinates theta and phi +-- are held constant. +-- Defined everywhere except at the origin. +-- The unit vector 'rHat' points in different directions at different points +-- in space. It is therefore better interpreted as a vector field, rather +-- than a vector. +rHat :: VectorField +rHat rv = d ^/ magnitude d + where + d = displacement (cart 0 0 0) rv + +-- | The vector field in which each point in space is associated +-- with a unit vector in the direction of increasing spherical coordinate +-- theta, while spherical coordinates r and phi are held constant. +-- Defined everywhere except on the z axis. +thetaHat :: VectorField +thetaHat r = vec (cos theta * cos phi) (cos theta * sin phi) (-sin theta) + where + (_,theta,phi) = sphericalCoordinates r + +-- | The vector field in which each point in space is associated +-- with a unit vector in the direction of increasing (cylindrical or spherical) coordinate +-- phi, while cylindrical coordinates s and z +-- (or spherical coordinates r and theta) are held constant. +-- Defined everywhere except on the z axis. +phiHat :: VectorField +phiHat r = vec (-sin phi) (cos phi) 0 + where + (_,phi,_) = cylindricalCoordinates r + +-- | The vector field in which each point in space is associated +-- with a unit vector in the direction of increasing cylindrical coordinate +-- s, while cylindrical coordinates phi and z +-- are held constant. +-- Defined everywhere except on the z axis. +sHat :: VectorField +sHat r = vec (cos phi) (sin phi) 0 + where + (_,phi,_) = cylindricalCoordinates r + +-- | The vector field in which each point in space is associated +-- with a unit vector in the direction of increasing Cartesian coordinate +-- x, while Cartesian coordinates y and z +-- are held constant. +-- Defined everywhere. +xHat :: VectorField +xHat = const iHat + +-- | The vector field in which each point in space is associated +-- with a unit vector in the direction of increasing Cartesian coordinate +-- y, while Cartesian coordinates x and z +-- are held constant. +-- Defined everywhere. +yHat :: VectorField +yHat = const jHat + +-- | The vector field in which each point in space is associated +-- with a unit vector in the direction of increasing Cartesian coordinate +-- z, while Cartesian coordinates x and y +-- are held constant. +-- Defined everywhere. +zHat :: VectorField +zHat = const kHat + diff --git a/src/Physics/Learn/RootFinding.hs b/src/Physics/Learn/RootFinding.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/RootFinding.hs @@ -0,0 +1,111 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Safe #-} + +{- | +Module : Physics.Learn.RootFinding +Copyright : (c) Scott N. Walck 2012-2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +Functions for approximately solving equations like f(x) = 0. +These functions proceed by assuming that f is continuous, +and that a root is bracketed. A bracket around a root consists +of numbers a, b such that f(a) f(b) <= 0. Since the product +changes sign, there must be an x with a < x < b such that f(x) = 0. +-} + +module Physics.Learn.RootFinding + ( findRoots + , findRootsN + , findRoot + , bracketRoot + , bracketRootStep + ) + where + +-- | Given an initial bracketing of a root +-- (an interval (a,b) for which f(a) f(b) <= 0), +-- produce a bracket of arbitrary smallness. +bracketRoot :: (Ord a, Fractional a) => + a -- ^ desired accuracy + -> (a -> a) -- ^ function + -> (a,a) -- ^ initial bracket + -> (a,a) -- ^ final bracket +bracketRoot dx f (a,b) + = let fa = f a + fb = f b + bRoot ((c,fc),(d,fd)) = let m = (c + d) / 2 + fm = f m + in if abs (c - d) < dx + then (c,d) + else if fc * fm <= 0 + then bRoot ((c,fc),(m,fm)) + else bRoot ((m,fm),(d,fd)) + in if fa * fb > 0 + then error "bracketRoot: initial interval is not a bracket" + else bRoot ((a,fa),(b,fb)) + +-- | Given a bracketed root, return a half-width bracket. +bracketRootStep :: (Ord a, Fractional a) => + (a -> a) -- ^ function + -> ((a,a),(a,a)) -- ^ original bracket + -> ((a,a),(a,a)) -- ^ new bracket +bracketRootStep f ((a,fa),(b,fb)) + = let m = (a + b) / 2 + fm = f m + in if fa * fm <= 0 + then ((a,fa),(m,fm)) + else ((m,fm),(b,fb)) + +findRootMachinePrecision :: (Double -> Double) + -> ((Double,Double),(Double,Double)) + -> Double +findRootMachinePrecision f ((c,fc),(d,fd)) + = let m = (c + d) / 2 + fm = f m + in if fc == 0 + then c + else if fd == 0 + then d + else if c == m + then c + else if m == d + then d + else if fc * fm <= 0 + then findRootMachinePrecision f ((c,fc),(m,fm)) + else findRootMachinePrecision f ((m,fm),(d,fd)) + +-- | Find a single root in a bracketed region. +-- The algorithm continues until it exhausts the +-- precision of a 'Double'. This could cause the function to hang. +findRoot :: (Double -> Double) -- ^ function + -> (Double,Double) -- ^ initial bracket + -> Double -- ^ approximate root +findRoot f (a,b) + = let fa = f a + fb = f b + in if fa * fb > 0 + then error "bracketRoot: initial interval is not a bracket" + else findRootMachinePrecision f ((a,fa),(b,fb)) + +-- | Find a list of roots for a function over a given range. +-- First parameter is the initial number of intervals to +-- use to find the roots. If roots are closely spaced, +-- this number of intervals may need to be large. +findRootsN :: Int -- ^ initial number of intervals to use + -> (Double -> Double) -- ^ function + -> (Double,Double) -- ^ range over which to search + -> [Double] -- ^ list of roots +findRootsN n f (a,b) + = let dx = (b - a) / fromIntegral n + xs = [a,a+dx..b] + in map (findRootMachinePrecision f) [((x0,fx0),(x1,fx1)) | (x0,x1) <- zip xs (tail xs), let fx0 = f x0, let fx1 = f x1, fx0 * fx1 <= 0] + +-- | Find a list of roots for a function over a given range. +-- There are no guarantees that all roots will be found. +-- Uses 'findRootsN' with 1000 intervals. +findRoots :: (Double -> Double) -- ^ function + -> (Double,Double) -- ^ range over which to search + -> [Double] -- ^ list of roots +findRoots = findRootsN 1000 diff --git a/src/Physics/Learn/RungeKutta.hs b/src/Physics/Learn/RungeKutta.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/RungeKutta.hs @@ -0,0 +1,63 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.RungeKutta +Copyright : (c) Scott N. Walck 2012-2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +Differential equation solving using 4th-order Runge-Kutta +-} + +module Physics.Learn.RungeKutta + ( rungeKutta4 + , integrateSystem + ) + where + +import Physics.Learn.StateSpace + ( StateSpace(..) + , Diff + , Time + , (.+^) + ) +import Data.VectorSpace + ( (^+^) + , (*^) + , (^/) + ) + +-- | Take a single 4th-order Runge-Kutta step +rungeKutta4 :: StateSpace p => (p -> Diff p) -> Time p -> p -> p +rungeKutta4 f dt y + = let k0 = dt *^ f y + k1 = dt *^ f (y .+^ k0 ^/ 2) + k2 = dt *^ f (y .+^ k1 ^/ 2) + k3 = dt *^ f (y .+^ k2) + in y .+^ (k0 ^+^ 2 *^ k1 ^+^ 2 *^ k2 ^+^ k3) ^/ 6 + +-- | Solve a first-order system of differential equations with 4th-order Runge-Kutta +integrateSystem :: StateSpace p => (p -> Diff p) -> Time p -> p -> [p] +integrateSystem systemDerivative dt + = iterate (rungeKutta4 systemDerivative dt) + + + +{- +-- | Take a single 4th-order Runge-Kutta step +rungeKutta4 :: (VectorSpace v, Fractional (Scalar v)) => (v -> v) -> Scalar v -> v -> v +rungeKutta4 f h y + = let k0 = h *^ f y + k1 = h *^ f (y ^+^ k0 ^/ 2) + k2 = h *^ f (y ^+^ k1 ^/ 2) + k3 = h *^ f (y ^+^ k2) + in y ^+^ (k0 ^+^ 2 *^ k1 ^+^ 2 *^ k2 ^+^ k3) ^/ 6 + +-- | Solve a first-order system of differential equations with 4th-order Runge-Kutta +integrateSystem :: (VectorSpace v, Fractional (Scalar v)) => (v -> v) -> Scalar v -> v -> [v] +integrateSystem systemDerivative dt + = iterate (rungeKutta4 systemDerivative dt) +-} diff --git a/src/Physics/Learn/SimpleVec.hs b/src/Physics/Learn/SimpleVec.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/SimpleVec.hs @@ -0,0 +1,117 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Safe #-} + +{- | +Module : Physics.Learn.SimpleVec +Copyright : (c) Scott N. Walck 2012-2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +Basic operations on the vector type 'Vec', such as vector addition +and scalar multiplication. +This module is simple in the sense that the operations +on vectors all have simple, concrete types, +without the need for type classes. +This makes using and reasoning about vector operations +easier for a person just learning Haskell. +-} + +-- 2011 Apr 10 +-- Placed the code common to SimpleVec and CarrotVec in CommonVec + +-- 2011 Mar 19 +-- Add support for sumV, so that the interface matches CarrotVec.hs + +-- This uses the same internal data representation as SimpleVector, +-- but uses an interface to match Conal Elliott's operators for +-- vectors. (A similar interface to CarrotVector and SimpleCarrotVector.) +-- The notation +-- zeroV, negateV, (^+^), (^-^) +-- is borrowed from Data.AdditiveGroup, and +-- (*^), (^*), (^/), (<.>), magnitude +-- is borrowed from Data.VectorSpace. +-- Cross product operator is my own. + +module Physics.Learn.SimpleVec + ( Vec + , xComp + , yComp + , zComp + , vec + , (^+^) + , (^-^) + , (*^) + , (^*) + , (^/) + , (<.>) + , (><) + , magnitude + , zeroV + , negateV + , sumV + , iHat + , jHat + , kHat + ) + where + +import Physics.Learn.CommonVec + ( Vec(..) + , vec + , iHat + , jHat + , kHat + , (><) + ) + +infixl 6 ^+^ +infixl 6 ^-^ +infixl 7 *^ +infixl 7 ^* +infixl 7 ^/ +infixl 7 <.> + +-- | The zero vector. +zeroV :: Vec +zeroV = vec 0 0 0 + +-- | The additive inverse of a vector. +negateV :: Vec -> Vec +negateV (Vec ax ay az) = Vec (-ax) (-ay) (-az) + +-- | Sum of a list of vectors. +sumV :: [Vec] -> Vec +sumV = foldr (^+^) zeroV + +-- | 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) + +-- | Scalar multiplication, where the scalar is on the left +-- and the vector is on the right. +(*^) :: Double -> Vec -> Vec +c *^ Vec ax ay az = Vec (c*ax) (c*ay) (c*az) + +-- | Scalar multiplication, where the scalar is on the right +-- and the vector is on the left. +(^*) :: Vec -> Double -> Vec +Vec ax ay az ^* c = Vec (c*ax) (c*ay) (c*az) + +-- | Division of a vector by a scalar. +(^/) :: Vec -> Double -> Vec +Vec ax ay az ^/ c = Vec (ax/c) (ay/c) (az/c) + +-- | Dot product of two vectors. +(<.>) :: Vec -> Vec -> Double +Vec ax ay az <.> Vec bx by bz = ax*bx + ay*by + az*bz + +-- | Magnitude of a vector. +magnitude :: Vec -> Double +magnitude v = sqrt(v <.> v) + diff --git a/src/Physics/Learn/StateSpace.hs b/src/Physics/Learn/StateSpace.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/StateSpace.hs @@ -0,0 +1,90 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE FlexibleContexts, FlexibleInstances, TypeFamilies #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.StateSpace +Copyright : (c) Scott N. Walck 2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +A 'StateSpace' is an affine space where the associated vector space +has scalars that are instances of 'Fractional'. +If p is an instance of 'StateSpace', then the associated vectorspace +'Diff' p is intended to represent the space of time derivatives +of paths in p. + +'StateSpace' is very similar to Conal Elliott's 'AffineSpace'. +-} + +module Physics.Learn.StateSpace + ( StateSpace(..) + , (.-^) + , Time + ) + where + +import Data.VectorSpace + ( VectorSpace + , Scalar + , negateV + ) +import Physics.Learn.Position + ( Position + , shiftPosition + , displacement + ) +import Physics.Learn.CarrotVec + ( Vec + , (^+^) + , (^-^) + ) + +infixl 6 .+^, .-^ +infix 6 .-. + +-- | A 'StateSpace' has an associated vector space, the vectors of which +-- can be multiplied or divided by scalars. +-- An example would be the set of positions of a particle. +-- Position is not a vector, but displacement (difference in position) is a vector. +class (VectorSpace (Diff p), Fractional (Scalar (Diff p))) => StateSpace p where + -- | Associated vector space + type Diff p + -- | Subtract points + (.-.) :: p -> p -> Diff p + -- | Point plus vector + (.+^) :: p -> Diff p -> p + +-- | The scalars of the associated vector space can be thought of as time intervals. +type Time p = Scalar (Diff p) + +-- | Point minus vector +(.-^) :: StateSpace p => p -> Diff p -> p +p .-^ v = p .+^ negateV v + +instance StateSpace Double where + type Diff Double = Double + (.-.) = (-) + (.+^) = (+) + +instance StateSpace Vec where + type Diff Vec = Vec + (.-.) = (^-^) + (.+^) = (^+^) + +instance StateSpace Position where + type Diff Position = Vec + (.-.) = flip displacement + (.+^) = flip shiftPosition + +instance (StateSpace p, StateSpace q, Time p ~ Time q) => StateSpace (p,q) where + type Diff (p,q) = (Diff p, Diff q) + (p,q) .-. (p',q') = (p .-. p', q .-. q') + (p,q) .+^ (u,v) = (p .+^ u, q .+^ v) + +instance (StateSpace p, StateSpace q, StateSpace r, Time p ~ Time q + ,Time q ~ Time r) => StateSpace (p,q,r) where + type Diff (p,q,r) = (Diff p, Diff q, Diff r) + (p,q,r) .-. (p',q',r') = (p .-. p', q .-. q', r .-. r') + (p,q,r) .+^ (u,v,w) = (p .+^ u, q .+^ v, r .+^ w) diff --git a/src/Physics/Learn/Surface.hs b/src/Physics/Learn/Surface.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/Surface.hs @@ -0,0 +1,290 @@ +{-# LANGUAGE TypeFamilies #-} +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.Surface +Copyright : (c) Scott N. Walck 2012-2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +This module contains functions for working with 'Surface's +and surface integrals over 'Surface's. +-} + +module Physics.Learn.Surface + ( Surface(..) + , unitSphere + , centeredSphere + , sphere + , northernHemisphere + , disk + , shiftSurface + , surfaceIntegral + , dottedSurfaceIntegral + ) + where + +import Data.VectorSpace + ( VectorSpace + , InnerSpace + , Scalar + ) +import Physics.Learn.CarrotVec + ( vec + , (^+^) + , (^-^) + , (^*) + , (^/) + , (<.>) + , (><) + , magnitude + , sumV + ) +import Physics.Learn.Position + ( Position + , Displacement + , VectorField + , Field + , cart + , cyl + , shiftPosition + , displacement + ) + +-- | Surface is a parametrized function from two parameters to space, +-- lower and upper limits on the first parameter, and +-- lower and upper limits for the second parameter +-- (expressed as functions of the first parameter). +data Surface = Surface { surfaceFunc :: (Double,Double) -> Position -- ^ function from two parameters (s,t) into space + , lowerLimit :: Double -- ^ s_l + , upperLimit :: Double -- ^ s_u + , lowerCurve :: Double -> Double -- ^ t_l(s) + , upperCurve :: Double -> Double -- ^ t_u(s) + } + +-- | A unit sphere, centered at the origin. +unitSphere :: Surface +unitSphere = Surface (\(th,phi) -> cart (sin th * cos phi) (sin th * sin phi) (cos th)) + 0 pi (const 0) (const $ 2*pi) + +-- | A sphere with given radius centered at the origin. +centeredSphere :: Double -> Surface +centeredSphere r = Surface (\(th,phi) -> cart (r * sin th * cos phi) (r * sin th * sin phi) (r * cos th)) + 0 pi (const 0) (const $ 2*pi) + +-- | Sphere with given radius and center. +sphere :: Double -> Position -> Surface +sphere r c = Surface (\(th,phi) -> shiftPosition (vec (r * sin th * cos phi) (r * sin th * sin phi) (r * cos th)) c) + 0 pi (const 0) (const $ 2*pi) + +-- | The upper half of a unit sphere, centered at the origin. +northernHemisphere :: Surface +northernHemisphere = Surface (\(th,phi) -> cart (sin th * cos phi) (sin th * sin phi) (cos th)) + 0 (pi/2) (const 0) (const $ 2*pi) + +-- | A disk with given radius, centered at the origin. +disk :: Double -> Surface +disk radius = Surface (\(s,phi) -> cyl s phi 0) 0 radius (const 0) (const (2*pi)) + +-- | A plane surface integral, in which area element is a scalar. +surfaceIntegral :: (VectorSpace v, Scalar v ~ Double) => + Int -- ^ number of intervals for first parameter, s + -> Int -- ^ number of intervals for second parameter, t + -> Field v -- ^ the scalar or vector field to integrate + -> Surface -- ^ the surface over which to integrate + -> v -- ^ the resulting scalar or vector +surfaceIntegral n1 n2 field (Surface f s_l s_u t_l t_u) + = sumV $ map sumV $ zipWith (zipWith (^*)) aveVals (map (map magnitude) areas) + where + pts = [[f (s,t) | t <- linSpaced n2 (t_l s) (t_u s)] | s <- linSpaced n1 s_l s_u] + areas = zipWith (zipWith (><)) dus dvs + dus = zipWith (zipWith displacement) pts (tail pts) + dvs = map (\row -> zipWith displacement row (tail row)) pts + vals = map (map field) pts + halfAveVals = map (\row -> zipWith ave (tail row) row) vals + aveVals = zipWith (zipWith ave) (tail halfAveVals) halfAveVals + +-- | A dotted surface integral, in which area element is a vector. +dottedSurfaceIntegral :: Int -- ^ number of intervals for first parameter, s + -> Int -- ^ number of intervals for second parameter, t + -> VectorField -- ^ the vector field to integrate + -> Surface -- ^ the surface over which to integrate + -> Double -- ^ the resulting scalar +dottedSurfaceIntegral n1 n2 vf (Surface f s_l s_u t_l t_u) + = sum $ map sum $ zipWith (zipWith (<.>)) aveVals areas + where + pts = [[f (s,t) | t <- linSpaced n2 (t_l s) (t_u s)] | s <- linSpaced n1 s_l s_u] + areas = zipWith (zipWith (><)) dus dvs + dus = zipWith (zipWith displacement) pts (tail pts) + dvs = map (\row -> zipWith displacement row (tail row)) pts + vals = map (map vf) pts + halfAveVals = map (\row -> zipWith ave (tail row) row) vals + aveVals = zipWith (zipWith ave) (tail halfAveVals) halfAveVals + +{- +evalSquare :: (InnerSpace v, Scalar v ~ Double) => Double -> Int -> Int + -> (Vec -> v) -> Surface + -> Vec -> Vec -> Vec -> Vec + -> v -> v -> v -> v -> v +evalSquare tol level maxlevel field (Surface f s_l s_u t_l t_u) + surfll surflu surful surfuu fieldll fieldlu fieldul fielduu val + = let s_m = (s_l + s_u) / 2 + t_m s = (t_l s + t_u s) / 2 + surflm = f (s_l,t_m s_l) + surfum = f (s_u,t_m s_u) + surfml = f (s_m,t_l s_m) + surfmu = f (s_m,t_u s_m) + surfmm = f (s_m,t_m s_m) + fieldlm = field surflm + fieldum = field surfum + fieldml = field surfml + fieldmu = field surfmu + fieldmm = field surfmm + dull = surfml ^-^ surfll + dulu = surfmm ^-^ surflm + duul = surful ^-^ surfml + duuu = surfum ^-^ surfmm + dvll = surflm ^-^ surfll + dvlu = surflu ^-^ surflm + dvul = surfmm ^-^ surfml + dvuu = surfmu ^-^ surfmm + areall = dull >< dvll + arealu = dulu >< dvlu + areaul = duul >< dvul + areauu = duuu >< dvuu + valll = average [fieldll,fieldlm,fieldml,fieldmm] <.> areall + vallu = average [fieldlm,fieldlu,fieldmm,fieldmu] <.> arealu + valul = average [fieldml,fieldmm,fieldul,fieldum] <.> areaul + valuu = average [fieldmm,fieldmu,fieldum,fielduu] <.> areauu + newval = valll ^+^ vallu ^+^ valul ^+^ valuu + in if magnitude (newval ^-^ val) < tol then + newval + else + evalSquare (tol/2) (level+1) maxlevel field (Surface f s_l s_m t_l t_m) + surfll surflm surfml surfmm fieldll fieldlm fieldml fieldmm valll ^+^ + evalSquare (tol/2) (level+1) maxlevel field (Surface f s_l s_m t_m t_u) + surflm surflu surfmm surfmu fieldlm fieldlu fieldmm fieldmu vallu ^+^ + evalSquare (tol/2) (level+1) maxlevel field (Surface f s_m s_u t_l t_m) + surfml surfmm surful surfum fieldml fieldmm fieldul fieldum valul ^+^ + evalSquare (tol/2) (level+1) maxlevel field (Surface f s_m s_u t_m t_u) + surfmm surfmu surfum surfuu fieldmm fieldmu fieldum fielduu valuu +-} + +{- +dottedSurfIntegral :: Double + -> (Vec -> Vec) -> Surface + -> Double +dottedSurfIntegral tol vf (Surface f s_l s_u t_l t_u) + = let surfll = f (s_l,t_l s_l) + surflu = f (s_l,t_u s_l) + surful = f (s_u,t_l s_u) + surfuu = f (s_u,t_u s_u) + fieldll = vf surfll + fieldlu = vf surflu + fieldul = vf surful + fielduu = vf surfuu + du = surful ^-^ surfll + dv = surflu ^-^ surfll + area = du >< dv + val = average [fieldll,fieldlu,fieldul,fielduu] <.> area + in evalSquare tol 1 2 20 vf (Surface f s_l s_u t_l t_u) + surfll surflu surful surfuu fieldll fieldlu fieldul fielduu val + +fullDottedSurfIntegral :: Double -> Int -> Int + -> (Vec -> Vec) -> Surface + -> Double +fullDottedSurfIntegral tol minlevel maxlevel vf (Surface f s_l s_u t_l t_u) + = let surfll = f (s_l,t_l s_l) + surflu = f (s_l,t_u s_l) + surful = f (s_u,t_l s_u) + surfuu = f (s_u,t_u s_u) + fieldll = vf surfll + fieldlu = vf surflu + fieldul = vf surful + fielduu = vf surfuu + du = surful ^-^ surfll + dv = surflu ^-^ surfll + area = du >< dv + val = average [fieldll,fieldlu,fieldul,fielduu] <.> area + in evalSquare tol 1 minlevel maxlevel vf (Surface f s_l s_u t_l t_u) + surfll surflu surful surfuu fieldll fieldlu fieldul fielduu val + +evalSquare :: Double -> Int -> Int -> Int + -> (Vec -> Vec) -> Surface + -> Vec -> Vec -> Vec -> Vec + -> Vec -> Vec -> Vec -> Vec -> Double -> Double +evalSquare tol level minlevel maxlevel field (Surface f s_l s_u t_l t_u) + surfll surflu surful surfuu fieldll fieldlu fieldul fielduu val + = let s_m = (s_l + s_u) / 2 + t_m s = (t_l s + t_u s) / 2 + surflm = f (s_l,t_m s_l) + surfum = f (s_u,t_m s_u) + surfml = f (s_m,t_l s_m) + surfmu = f (s_m,t_u s_m) + surfmm = f (s_m,t_m s_m) + fieldlm = field surflm + fieldum = field surfum + fieldml = field surfml + fieldmu = field surfmu + fieldmm = field surfmm + dull = surfml ^-^ surfll + dulu = surfmm ^-^ surflm + duul = surful ^-^ surfml + duuu = surfum ^-^ surfmm + dvll = surflm ^-^ surfll + dvlu = surflu ^-^ surflm + dvul = surfmm ^-^ surfml + dvuu = surfmu ^-^ surfmm + areall = dull >< dvll + arealu = dulu >< dvlu + areaul = duul >< dvul + areauu = duuu >< dvuu + valll = average [fieldll,fieldlm,fieldml,fieldmm] <.> areall + vallu = average [fieldlm,fieldlu,fieldmm,fieldmu] <.> arealu + valul = average [fieldml,fieldmm,fieldul,fieldum] <.> areaul + valuu = average [fieldmm,fieldmu,fieldum,fielduu] <.> areauu + newval = valll + vallu + valul + valuu + in if level >= maxlevel || level >= minlevel && abs (newval - val) < tol then + newval + else + evalSquare (tol/4) (level+1) minlevel maxlevel field (Surface f s_l s_m t_l t_m) + surfll surflm surfml surfmm fieldll fieldlm fieldml fieldmm valll + + evalSquare (tol/4) (level+1) minlevel maxlevel field (Surface f s_l s_m t_m t_u) + surflm surflu surfmm surfmu fieldlm fieldlu fieldmm fieldmu vallu + + evalSquare (tol/4) (level+1) minlevel maxlevel field (Surface f s_m s_u t_l t_m) + surfml surfmm surful surfum fieldml fieldmm fieldul fieldum valul + + evalSquare (tol/4) (level+1) minlevel maxlevel field (Surface f s_m s_u t_m t_u) + surfmm surfmu surfum surfuu fieldmm fieldmu fieldum fielduu valuu +-} + +-- n+1 points +linSpaced :: Int -> Double -> Double -> [Double] +linSpaced n a b + | a < b = let dx = (b - a) / fromIntegral n + in [a,a+dx..b] + | a ~~ b = [ave a b] + | otherwise = error $ "linSpaced: lower limit greater than upper limit: (n,a,b) = " ++ show (n,a,b) + +(~~) :: (InnerSpace v, Scalar v ~ Double) => v -> v -> Bool +a ~~ b = magnitude (a ^-^ b) < tolerance + +tolerance :: Double +tolerance = 1e-10 + +ave :: (VectorSpace v, Scalar v ~ Double) => v -> v -> v +ave v1 v2 = (v1 ^+^ v2) ^/ 2 + +{- +average :: (VectorSpace v, Scalar v ~ Double) => [v] -> v +average vs = sumV vs ^/ fromIntegral (length vs) + +areaOfSurface :: Surface -> Double +areaOfSurface = surfaceIntegral 100 100 (const 1) +-} + +-- | Shift a surface by a displacement. +shiftSurface :: Displacement -> Surface -> Surface +shiftSurface d (Surface f sl su tl tu) + = Surface (shiftPosition d . f) sl su tl tu diff --git a/src/Physics/Learn/Volume.hs b/src/Physics/Learn/Volume.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/Volume.hs @@ -0,0 +1,175 @@ +{-# LANGUAGE TypeFamilies #-} +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.Volume +Copyright : (c) Scott N. Walck 2012-2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +This module contains functions for working with 'Volume's +and volume integrals over 'Volume's. +-} + +module Physics.Learn.Volume + ( Volume(..) + , unitBall + , unitBallCartesian + , centeredBall + , ball + , northernHalfBall + , centeredCylinder + , shiftVolume + , volumeIntegral + ) + where + +import Data.VectorSpace + ( VectorSpace + , InnerSpace + , Scalar + ) +import Physics.Learn.CarrotVec + ( Vec + , vec + , sumV + , (^+^) + , (^-^) + , (^*) + , (^/) + , (<.>) + , (><) + , magnitude + ) +import Physics.Learn.Position + ( Position + , Displacement + , Field + , cartesian + , cylindrical + , spherical + , shiftPosition + , displacement + ) + +-- | Volume is a parametrized function from three parameters to space, +-- lower and upper limits on the first parameter, +-- lower and upper limits for the second parameter +-- (expressed as functions of the first parameter), +-- and lower and upper limits for the third parameter +-- (expressed as functions of the first and second parameters). +data Volume = Volume { volumeFunc :: (Double,Double,Double) -> Position -- ^ function from 3 parameters to space + , loLimit :: Double -- ^ s_a + , upLimit :: Double -- ^ s_b + , loCurve :: Double -> Double -- ^ t_a(s) + , upCurve :: Double -> Double -- ^ t_b(s) + , loSurf :: Double -> Double -> Double -- ^ u_a(s,t) + , upSurf :: Double -> Double -> Double -- ^ u_b(s,t) + } + +-- | A unit ball, centered at the origin. +unitBall :: Volume +unitBall = Volume spherical 0 1 (const 0) (const pi) (\_ _ -> 0) (\_ _ -> 2*pi) + +-- | A unit ball, centered at the origin. Specified in Cartesian coordinates. +unitBallCartesian :: Volume +unitBallCartesian = Volume cartesian (-1) 1 (\x -> -sqrtTol (1 - x*x)) (\x -> sqrtTol (1 - x*x)) + (\x y -> -sqrtTol (1 - x*x - y*y)) (\x y -> sqrtTol (1 - x*x - y*y)) + +-- | A ball with given radius, centered at the origin. +centeredBall :: Double -> Volume +centeredBall a = Volume spherical 0 a (const 0) (const pi) (\_ _ -> 0) (\_ _ -> 2*pi) + +-- | Ball with given radius and center. +ball :: Double -- ^ radius + -> Position -- ^ center + -> Volume -- ^ ball with given radius and center +ball a c = Volume (\(r,th,phi) -> shiftPosition (vec (r * sin th * cos phi) (r * sin th * sin phi) (r * cos th)) c) + 0 a (const 0) (const pi) (\_ _ -> 0) (\_ _ -> 2*pi) + +-- | Upper half ball, unit radius, centered at origin. +northernHalfBall :: Volume +northernHalfBall = Volume spherical 0 1 (const 0) (const (pi/2)) (\_ _ -> 0) (\_ _ -> 2*pi) + +-- | Cylinder with given radius and height. Circular base of the cylinder +-- is centered at the origin. Circular top of the cylinder lies in plane z = h. +centeredCylinder :: Double -- radius + -> Double -- height + -> Volume -- cylinder +centeredCylinder r h = Volume cylindrical 0 r (const 0) (const (2*pi)) (\_ _ -> 0) (\_ _ -> h) + +-- | A volume integral +volumeIntegral :: (VectorSpace v, Scalar v ~ Double) => + Int -- ^ number of intervals for first parameter (s) + -> Int -- ^ number of intervals for second parameter (t) + -> Int -- ^ number of intervals for third parameter (u) + -> Field v -- ^ scalar or vector field + -> Volume -- ^ the volume + -> v -- ^ scalar or vector result +volumeIntegral n1 n2 n3 field (Volume f s_l s_u t_l t_u u_l u_u) + = sumV $ map sumV $ map (map sumV) (zipCubeWith (^*) aveVals volumes) + where + pts = [[[f (s,t,u) | u <- linSpaced n3 (u_l s t) (u_u s t) ] | t <- linSpaced n2 (t_l s) (t_u s)] | s <- linSpaced n1 s_l s_u] + volumes = zipWith3 (zipWith3 (zipWith3 (\du dv dw -> du <.> (dv >< dw)))) dus dvs dws + dus = uncurry zipSub3 (shift1 pts) + dvs = uncurry zipSub3 (shift2 pts) + dws = uncurry zipSub3 (shift3 pts) + vals = map (map (map field)) pts + aveVals = ((uncurry zipAve3 . shift1) . (uncurry zipAve3 . shift2) . (uncurry zipAve3 . shift3)) vals + +-- zipSquareWith :: (a -> b -> c) -> [[a]] -> [[b]] -> [[c]] +-- zipSquareWith = zipWith . zipWith + +zipCubeWith :: (a -> b -> c) -> [[[a]]] -> [[[b]]] -> [[[c]]] +zipCubeWith = zipWith . zipWith . zipWith + +-- zipSub :: [Vec] -> [Vec] -> [Vec] +-- zipSub = zipWith (^-^) + +-- zipSub2 :: [[Vec]] -> [[Vec]] -> [[Vec]] +-- zipSub2 = zipWith $ zipWith (^-^) + +zipSub3 :: [[[Position]]] -> [[[Position]]] -> [[[Vec]]] +zipSub3 = zipCubeWith displacement + +zipAve3 :: (VectorSpace v, Scalar v ~ Double) => [[[v]]] -> [[[v]]] -> [[[v]]] +zipAve3 = zipCubeWith ave + +shift1 :: [a] -> ([a],[a]) +shift1 pts = (pts, tail pts) + +shift2 :: [[a]] -> ([[a]],[[a]]) +shift2 pts2d = (pts2d, map tail pts2d) + +shift3 :: [[[a]]] -> ([[[a]]],[[[a]]]) +shift3 pts3d = (pts3d, map (map tail) pts3d) + +-- | n+1 points +linSpaced :: Int -> Double -> Double -> [Double] +linSpaced n a b + | a < b = let dx = (b - a) / fromIntegral n + in [a,a+dx..b] + | a ~~ b = [ave a b] + | otherwise = error $ "linSpaced: lower limit greater than upper limit: (n,a,b) = " ++ show (n,a,b) + +(~~) :: (InnerSpace v, Scalar v ~ Double) => v -> v -> Bool +a ~~ b = magnitude (a ^-^ b) < tolerance + +tolerance :: Double +tolerance = 1e-10 + +ave :: (VectorSpace v, Scalar v ~ Double) => v -> v -> v +ave v1 v2 = (v1 ^+^ v2) ^/ 2 + +sqrtTol :: Double -> Double +sqrtTol x + | x >= 0 = sqrt x + | abs x <= tolerance = 0 + | otherwise = error ("sqrtTol: I can't take the sqrt of " ++ show x) + +-- | Shift a volume by a displacement. +shiftVolume :: Displacement -> Volume -> Volume +shiftVolume d (Volume f sl su tl tu ul uu) + = Volume (shiftPosition d . f) sl su tl tu ul uu