learn-physics 0.5 → 0.5.2
raw patch · 13 files changed
+453/−359 lines, 13 filesdep +learn-physicsdep ~basedep ~glossdep ~not-glosssetup-changednew-component:exe:learn-physics-BCircularLoopnew-component:exe:learn-physics-LorentzForceSimulationnew-component:exe:learn-physics-PlaneWavenew-component:exe:learn-physics-Projectilenew-component:exe:learn-physics-eFieldLine2Dnew-component:exe:learn-physics-eFieldLine3Dnew-component:exe:learn-physics-sunEarth
Dependencies added: learn-physics
Dependency ranges changed: base, gloss, not-gloss, vector-space
Files
- LICENSE +1/−1
- Setup.hs +2/−0
- examples/src/BCircularLoop.hs +30/−0
- examples/src/LorentzForceSimulation.hs +64/−0
- examples/src/PlaneWave.hs +48/−0
- examples/src/Projectile.hs +77/−0
- examples/src/eFieldLine2D.hs +45/−0
- examples/src/eFieldLine3D.hs +48/−0
- examples/src/sunEarthRK4.hs +86/−0
- learn-physics.cabal +51/−8
- src/Physics/Learn/AdaptiveQuadrature.hs +0/−294
- src/Physics/Learn/Visual/VisTools.hs +1/−1
- src/Tests.hs +0/−55
LICENSE view
@@ -1,4 +1,4 @@-Copyright (c) 2011-2014 Scott N. Walck <walck@lvc.edu>.+Copyright (c) 2011-2015 Scott N. Walck <walck@lvc.edu>. All rights reserved. Redistribution and use in source and binary forms, with or without
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ examples/src/BCircularLoop.hs view
@@ -0,0 +1,30 @@+{-# OPTIONS_GHC -Wall #-}++module Main where++import Physics.Learn+import Vis++loopCurve :: Curve+loopCurve = Curve (\phi -> cyl 1 phi 0) 0 (2*pi)++loop :: CurrentDistribution+loop = LineCurrent 20 loopCurve++samplePoints :: [Position]+samplePoints = [cyl s phi z |+ s <- [0.25,0.75..1.75]+ , phi <- [pi/6,pi/2..2*pi]+ , z <- [-1.5,-1..1.5]]++arrows :: VisObject Double+arrows = displayVectorField blue 5e-5 samplePoints (bField loop)++drawFun :: VisObject Double+drawFun = VisObjects [curveObject red loopCurve, arrows]++myOptions :: Options+myOptions = defaultOpts {optWindowName = "Magnetic Field from a Current Loop"}++main :: IO ()+main = display myOptions drawFun
+ examples/src/LorentzForceSimulation.hs view
@@ -0,0 +1,64 @@+module Main where++import Physics.Learn+import Vis+import SpatialMath+ ( Euler(..)+ )++drawFunction :: SimpleState -> VisObject Double+drawFunction (_t,r,_v)+ = RotEulerDeg (Euler 270 0 0) $ RotEulerDeg (Euler 0 180 0) $+ VisObjects [ Axes (0.5, 15)+ , Trans (v3FromPos r) (Sphere 0.1 Solid red)+ ]++statePropagationFunction :: Float -> SimpleState -> SimpleState+statePropagationFunction t' (t,r,v) = rungeKutta4 newton2 (realToFrac t' - t) (t,r,v)++-- Newton's Second Law+newton2 :: SimpleState -> Diff SimpleState+newton2 (t,r,v) = (1,v,force (t,r,v) ^/ m)++-- Lorentz Force Law+force :: SimpleState -> Vec+force (_t,r,v) = q *^ (electricField r ^+^ v >< magneticField r)++myOptions :: Options+myOptions = defaultOpts {optWindowName = "Particle Experiencing Electromagnetic Force"}++main :: IO ()+main = simulate+ myOptions+ 0.01+ (0,initialPosition,initialVelocity)+ drawFunction+ statePropagationFunction++-- particle mass+m :: Double+m = 1++-- particle charge+q :: Double+q = 1++-- Electric Field+electricField :: VectorField+electricField r = vec 0 2 0+ where+ (x,y,z) = cartesianCoordinates r++-- Magnetic Field+magneticField :: VectorField+magneticField r = vec 0 0 4+ where+ (x,y,z) = cartesianCoordinates r++-- Initial displacement+initialPosition :: Position+initialPosition = cart 0 0 0++-- Initial velocity+initialVelocity :: Vec+initialVelocity = vec 0 0 0
+ examples/src/PlaneWave.hs view
@@ -0,0 +1,48 @@+{-# OPTIONS_GHC -Wall #-}++module Main where++import Vis+ ( animate+ , VisObject(..)+ , red+ , blue+ , Options(..)+ , defaultOpts+ )+import Physics.Learn.CarrotVec+ ( vec+ )+import Physics.Learn+ ( Position+ , VectorField+ , displayVectorField+ , cart+ , cartesianCoordinates+ )++samplePoints :: [Position]+samplePoints = [cart x y z | x <- [-2,0,2], y <- [-2,0,2], z <- [-4,-3.6..4]]++drawFun :: Float -> VisObject Double+drawFun time = VisObjects [displayVectorField blue 1 samplePoints (eField t)+ ,displayVectorField red 1 samplePoints (bField t)+ ]+ where+ t = realToFrac time++eField :: Double -> VectorField+eField t r = vec (cos (z - t)) 0 0+ where+ (_,_,z) = cartesianCoordinates r++bField :: Double -> VectorField+bField t r = vec 0 (cos (z - t)) 0+ where+ (_,_,z) = cartesianCoordinates r++myOptions :: Options+myOptions = defaultOpts {optWindowName = "Plane Wave"}++main :: IO ()+main = animate myOptions drawFun
+ examples/src/Projectile.hs view
@@ -0,0 +1,77 @@+{-# OPTIONS_GHC -Wall #-}++module Main where++import Graphics.Gnuplot.Simple+import Physics.Learn++--type StateTuple = (Double,Vec,Vec)+--type AccelerationFunction = StateTuple -> Vec++--eulerCromerSolution :: Double -> AccelerationFunction+-- -> StateTuple -> StateTuple+--eulerCromerSolution++-- vertical direction is y direction+earthSurfaceGravity :: OneParticleAccelerationFunction+earthSurfaceGravity _state = vec 0 (-g) 0++g :: Double+g = 9.81++projectileTuples :: Double -> Double+ -> OneParticleAccelerationFunction+ -> [OneParticleSystemState]+projectileTuples v0 theta af+ = oneParticleRungeKuttaSolution af 0.01+ (0, St (cart 0 0 0) (vec vx0 vy0 0))+ where+ vx0 = v0 * cos theta+ vy0 = v0 * sin theta++yCoord :: Position -> Double+yCoord r = y+ where+ (_,y,_) = cartesianCoordinates r++inAir :: [OneParticleSystemState] -> [OneParticleSystemState]+inAir = takeWhile (\(_,St r _) -> yCoord r >= 0)++initialProjState :: Double -> Double -> OneParticleSystemState+initialProjState v0 theta+ = (0, St (cart 0 0 0) (vec vx0 vy0 0))+ where+ vx0 = v0 * cos theta+ vy0 = v0 * sin theta++-- air resistance quadratic in the speed+surfaceGravityAirResistance :: Double -> Double+ -> OneParticleAccelerationFunction+surfaceGravityAirResistance m b (_t,St _r v)+ = netForce ^/ m+ where+ netForce = gravity ^+^ airResistance+ gravity = vec 0 (-m * g) 0+ airResistance = ((-b) * magnitude v) *^ v++trajectory :: [OneParticleSystemState] -> [(Double,Double)]+trajectory sts = [(x,y) | (_,St r _) <- sts, let (x,y,_) = cartesianCoordinates r]++traj :: Double -> [(Double,Double)]+traj b = trajectory $ inAir+ $ oneParticleRungeKuttaSolution+ (surfaceGravityAirResistance 2 b)+ 0.01+ (initialProjState 30 (pi/6))++main :: IO ()+main = plotPathsStyle+ [Title "Trajectories of 2-kg object, initial speed 30 m/s, angle 30 degrees"+ ,XLabel "Range (m)"+ ,YLabel "Height (m)"+ ,PNG "learn-physics-Projectile.png"+ ] [(defaultStyle {lineSpec = CustomStyle [LineTitle "No air resistance"]}, traj 0)+ ,(defaultStyle {lineSpec = CustomStyle [LineTitle "Drag 0.01 kg/m"]}, traj 0.01)+ ,(defaultStyle {lineSpec = CustomStyle [LineTitle "Drag 0.02 kg/m"]}, traj 0.02)+ ] >> putStrLn "output sent to file learn-physics-Projectile.png"+
+ examples/src/eFieldLine2D.hs view
@@ -0,0 +1,45 @@+{-# OPTIONS_GHC -Wall #-}++module Main where++import Graphics.Gloss+import Physics.Learn.CarrotVec+import Physics.Learn.Position+import Physics.Learn.Curve+import Physics.Learn.Charge+import Physics.Learn.Visual.GlossTools++pixelsPerMeter :: Float+pixelsPerMeter = 40++pixelsPerVPM :: Float+pixelsPerVPM = 5.6++scalePoint :: Float -> Point -> Point+scalePoint m (x,y) = (m*x,m*y)++twoD :: Vec -> Point+twoD r = (realToFrac $ xComp r,realToFrac $ yComp r)++twoDp :: Position -> Point+twoDp r = (realToFrac x, realToFrac y)+ where+ (x,y,_) = cartesianCoordinates r++samplePoints :: [Position]+samplePoints = [cart x y 0 | x <- [-8,-6..8], y <- [-6,-4..6], abs y > 0.5 || abs x > 4.5]++curve1 :: Curve+curve1 = Curve (\t -> cart t 0 0) (-4) 4++eFields :: [(Position,Vec)]+eFields = [(r,eFieldFromLineCharge (const 1e-9) curve1 r) | r <- samplePoints]++arrows :: [Picture]+arrows = [thickArrow 5 (scalePoint pixelsPerMeter $ twoDp r)+ (scalePoint pixelsPerVPM $ twoD e) | (r,e) <- eFields]++main :: IO ()+main = display (InWindow "Electric Field from a Line Charge" (680,520) (10,10)) white $+ Pictures [(Color blue (Pictures arrows))+ ,Color orange $ Line [(-4*pixelsPerMeter,0),(4*pixelsPerMeter,0)]]
+ examples/src/eFieldLine3D.hs view
@@ -0,0 +1,48 @@+{-# OPTIONS_GHC -Wall #-}++module Main where++import Vis+ ( display+ , VisObject(..)+ , red+ , blue+ , Options(..)+ , defaultOpts+ )+import Physics.Learn.Visual.VisTools+ ( curveObject+ , displayVectorField+ )+import Physics.Learn.Position+ ( Position+ , cart+ )+import Physics.Learn.Curve+ ( Curve(..)+ )+import Physics.Learn.Charge+ ( ChargeDistribution(..)+ , eField+ )++curve1 :: Curve+curve1 = Curve (\t -> cart t 0 0) (-4) 4++lineCharge :: ChargeDistribution+lineCharge = LineCharge (const 1e-9) curve1++samplePoints :: [Position]+samplePoints = [cart x y z | x <- [-8,-6..8], y <- [-4,-2..4], z <- [-4,-2..4], abs y + abs z > 0.5 || abs x > 4.5]++arrows :: VisObject Double+arrows = displayVectorField blue 10 samplePoints (eField lineCharge)++drawFun :: VisObject Double+drawFun = VisObjects [curveObject red curve1, arrows]++myOptions :: Options+myOptions = defaultOpts {optWindowName = "Electric Field from a Line Charge"}++main :: IO ()+main = display myOptions drawFun
+ examples/src/sunEarthRK4.hs view
@@ -0,0 +1,86 @@+{-# OPTIONS_GHC -Wall #-}++-- Animation of Earth orbiting around a fixed Sun+-- Using SI units++module Main where++import Physics.Learn+import Graphics.Gloss+import Graphics.Gloss.Data.ViewPort++type Acceleration = Vec++gGrav :: Double+gGrav = 6.67e-11++massSun :: Double+massSun = 1.99e30++-- This is enlarged so we can see it.+radiusSun :: Double+radiusSun = 0.1 * earthSunDistance++-- This is enlarged so we can see it.+radiusEarth :: Double+radiusEarth = 0.05 * earthSunDistance++earthSunDistance :: Double+earthSunDistance = 1.496e11++year :: Double+year = 365.25*24*60*60++-- Derived constants++initialEarthSpeed :: Double+initialEarthSpeed = 2*pi*earthSunDistance/year++initialState :: SimpleState+initialState = (0+ ,cart earthSunDistance 0 0+ ,vec 0 initialEarthSpeed 0)++rS :: Position+rS = cart 0 0 0++earthGravity :: SimpleAccelerationFunction+earthGravity (_,rE,_)+ = ((-gGrav) * massSun) *^ disp ^/ magnitude disp ** 3+ where+ disp = displacement rS rE++diskPic :: Double -> Picture+diskPic r = ThickCircle (radius/2) radius+ where radius = realToFrac r++-- A yellow disk will represent the Sun+yellowDisk :: Picture+yellowDisk = Color yellow (diskPic radiusSun)++-- A blue disk will represent the Earth+blueDisk :: Picture+blueDisk = Color blue (diskPic radiusEarth)++worldToPicture :: SimpleState -> Picture+worldToPicture (_,rE,_)+ = scale scl scl $ pictures [yellowDisk+ ,translate xE yE blueDisk+ ]+ where+ xE = realToFrac x+ yE = realToFrac y+ scl = 200 / realToFrac (earthSunDistance)+ (x,y,_) = cartesianCoordinates rE++timeScale :: Double+timeScale = 0.25 * year++simStep :: ViewPort -> Float -> SimpleState -> SimpleState+simStep _ dt = simpleRungeKuttaStep earthGravity dtScaled+ where+ dtScaled = timeScale * realToFrac dt++main :: IO ()+main = simulate (InWindow "Sun-Earth Animation" (1024, 768) (0, 0))+ black 50 initialState worldToPicture simStep
learn-physics.cabal view
@@ -1,5 +1,5 @@ Name: learn-physics-Version: 0.5+Version: 0.5.2 Synopsis: Haskell code for learning physics Description: A library of functions for vector calculus, calculation of electric field, electric flux,@@ -12,7 +12,7 @@ Category: Physics Build-type: Simple Cabal-version: >=1.8-Tested-with: GHC == 7.8.2+Tested-with: GHC == 7.10.2 Library Exposed-modules: Physics.Learn.Charge Physics.Learn.Current@@ -34,14 +34,57 @@ Physics.Learn.Visual.PlotTools Physics.Learn.Visual.VisTools Physics.Learn.Visual.GlossTools- Build-depends: base >= 4.2 && < 4.8,- vector-space >= 0.8.4 && < 0.9,- not-gloss >= 0.6 && < 0.7,+ Build-depends: base >= 4.2 && < 4.9,+ vector-space >= 0.8.4 && < 0.11,+ not-gloss >= 0.7.4 && < 0.8, spatial-math >= 0.2 && < 0.3,- gloss >= 1.8 && < 1.9,+ gloss >= 1.8 && < 1.10, gnuplot >= 0.5 && < 0.6 Hs-source-dirs: src Source-repository head- type: darcs- location: http://hub.darcs.net/scottwalck/learn-physics+ type: git+ location: https://github.com/walck/learn-physics++Executable learn-physics-PlaneWave+ Main-is: examples/src/PlaneWave.hs+ Build-depends: not-gloss >= 0.7.4 && < 0.8,+ base >= 4.5 && < 4.9,+ learn-physics++Executable learn-physics-eFieldLine3D+ Main-is: examples/src/eFieldLine3D.hs+ Build-depends: not-gloss >= 0.7.4 && < 0.8,+ base >= 4.5 && < 4.9,+ learn-physics++Executable learn-physics-LorentzForceSimulation+ Main-is: examples/src/LorentzForceSimulation.hs+ Build-depends: not-gloss >= 0.7.4 && < 0.8,+ spatial-math >= 0.2 && < 0.3,+ base >= 4.5 && < 4.9,+ learn-physics++Executable learn-physics-BCircularLoop+ Main-is: examples/src/BCircularLoop.hs+ Build-depends: not-gloss >= 0.7.4 && < 0.8,+ base >= 4.5 && < 4.9,+ learn-physics++Executable learn-physics-sunEarth+ Main-is: examples/src/sunEarthRK4.hs+ Build-depends: gloss >= 1.8 && < 1.10,+ base >= 4.5 && < 4.9,+ learn-physics++Executable learn-physics-eFieldLine2D+ Main-is: examples/src/eFieldLine2D.hs+ Build-depends: gloss >= 1.8 && < 1.10,+ base >= 4.5 && < 4.9,+ learn-physics++Executable learn-physics-Projectile+ Main-is: examples/src/Projectile.hs+ Build-depends: gnuplot >= 0.5 && < 0.6,+ base >= 4.5 && < 4.9,+ learn-physics
− src/Physics/Learn/AdaptiveQuadrature.hs
@@ -1,294 +0,0 @@-{-# 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/Visual/VisTools.hs view
@@ -63,7 +63,7 @@ -- | A displayable VisObject for a curve. curveObject :: Color -> Curve -> VisObject Double curveObject color (Curve f a b)- = Line' [(v3FromPos (f t), color) | t <- [a,a+(b-a)/1000..b]]+ = Line' Nothing [(v3FromPos (f t), color) | t <- [a,a+(b-a)/1000..b]] -- | Place a vector at a particular position. oneVector :: Color -> Position -> Vec -> VisObject Double
− src/Tests.hs
@@ -1,55 +0,0 @@-{-# OPTIONS_GHC -Wall #-}--module Main where--import Physics.Learn-import Test.QuickCheck--propGaussLaw1 :: (Double,Double,Double) -> Bool-propGaussLaw1 (x,y,z) = abs (eFlux - q/epsilon0) < 0.01- where- eFlux = fluxThroughLargeCenteredSphere (x,y,z) q- epsilon0 = 1 / (4 * pi * 9e9)- q = epsilon0--fluxThroughLargeCenteredSphere :: (Double,Double,Double) -> Double -> Double-fluxThroughLargeCenteredSphere (x,y,z) q- = electricFlux (centeredSphere radius) (PointCharge q (cart x y z))- where- radius = 2 * sqrt(x*x + y*y + z*z) + 1--currentLoop :: Double -> Current -> CurrentDistribution-currentLoop radius i- = LineCurrent i (Curve (\phi -> cyl radius phi 0) 0 (2*pi))--amperianLoop :: Double -> Curve-amperianLoop radius- = Curve (\t -> cart (radius + radius * sin t) 0 (radius * cos t)) 0 (2*pi)--magCirculation :: Double -> Current -> Double-magCirculation radius i- = dottedLineIntegral 20- (bFieldFromCurrentLoop i (Curve (\phi -> cyl radius phi 0) 0 (2*pi)))- (amperianLoop radius)--bFieldFromCurrentLoop :: Current -> Curve -> VectorField-bFieldFromCurrentLoop i c r- = k *^ crossedLineIntegral 20 integrand c- where- k = 1e-7 -- mu0 / (4 * pi)- integrand r' = (-i) *^ d ^/ magnitude d ** 3- where- d = displacement r' r--propAmpere1 :: Double -> Property-propAmpere1 radius- = radius > 0 ==> abs (magCirculation radius i - 4*pi*1e-7 * i) < 0.01- where- i = 1 / (4*pi*1e-7)--main :: IO ()-main = putStrLn "Gauss's law test:" >>- quickCheck propGaussLaw1 >>- putStrLn "Ampere's law test:" >>- quickCheck propAmpere1-