packages feed

learn-physics (empty) → 0.2

raw patch · 18 files changed

+2469/−0 lines, 18 filesdep +basedep +vector-space

Dependencies added: base, vector-space

Files

+ LICENSE view
@@ -0,0 +1,29 @@+Copyright (c) 2011-2014 Scott N. Walck <walck@lvc.edu>.+All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Scott N. Walck nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ learn-physics.cabal view
@@ -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 <walck@lvc.edu>+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
+ src/Physics/Learn/AdaptiveQuadrature.hs view
@@ -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)+
+ src/Physics/Learn/CarrotVec.hs view
@@ -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 <walck@lvc.edu>+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+
+ src/Physics/Learn/Charge.hs view
@@ -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 <walck@lvc.edu>+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 --+---------------------------------+
+ src/Physics/Learn/CommonVec.hs view
@@ -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 <walck@lvc.edu>+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
+ src/Physics/Learn/CompositeQuadrature.hs view
@@ -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 <walck@lvc.edu>+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
+ src/Physics/Learn/CoordinateFields.hs view
@@ -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 <walck@lvc.edu>+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+
+ src/Physics/Learn/CoordinateSystem.hs view
@@ -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 <walck@lvc.edu>+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)
+ src/Physics/Learn/Current.hs view
@@ -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 <walck@lvc.edu>+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+
+ src/Physics/Learn/Curve.hs view
@@ -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 <walck@lvc.edu>+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+
+ src/Physics/Learn/Position.hs view
@@ -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 <walck@lvc.edu>+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+
+ src/Physics/Learn/RootFinding.hs view
@@ -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 <walck@lvc.edu>+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
+ src/Physics/Learn/RungeKutta.hs view
@@ -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 <walck@lvc.edu>+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)+-}
+ src/Physics/Learn/SimpleVec.hs view
@@ -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 <walck@lvc.edu>+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)+
+ src/Physics/Learn/StateSpace.hs view
@@ -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 <walck@lvc.edu>+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)
+ src/Physics/Learn/Surface.hs view
@@ -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 <walck@lvc.edu>+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
+ src/Physics/Learn/Volume.hs view
@@ -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 <walck@lvc.edu>+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