diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) 2011-2014 Scott N. Walck . +Copyright (c) 2011-2015 Scott N. Walck . All rights reserved. Redistribution and use in source and binary forms, with or without diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/examples/src/BCircularLoop.hs b/examples/src/BCircularLoop.hs new file mode 100644 --- /dev/null +++ b/examples/src/BCircularLoop.hs @@ -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 diff --git a/examples/src/LorentzForceSimulation.hs b/examples/src/LorentzForceSimulation.hs new file mode 100644 --- /dev/null +++ b/examples/src/LorentzForceSimulation.hs @@ -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 diff --git a/examples/src/PlaneWave.hs b/examples/src/PlaneWave.hs new file mode 100644 --- /dev/null +++ b/examples/src/PlaneWave.hs @@ -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 diff --git a/examples/src/Projectile.hs b/examples/src/Projectile.hs new file mode 100644 --- /dev/null +++ b/examples/src/Projectile.hs @@ -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" + diff --git a/examples/src/eFieldLine2D.hs b/examples/src/eFieldLine2D.hs new file mode 100644 --- /dev/null +++ b/examples/src/eFieldLine2D.hs @@ -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)]] diff --git a/examples/src/eFieldLine3D.hs b/examples/src/eFieldLine3D.hs new file mode 100644 --- /dev/null +++ b/examples/src/eFieldLine3D.hs @@ -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 diff --git a/examples/src/sunEarthRK4.hs b/examples/src/sunEarthRK4.hs new file mode 100644 --- /dev/null +++ b/examples/src/sunEarthRK4.hs @@ -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 diff --git a/learn-physics.cabal b/learn-physics.cabal --- a/learn-physics.cabal +++ b/learn-physics.cabal @@ -1,5 +1,5 @@ Name: learn-physics -Version: 0.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 diff --git a/src/Physics/Learn/AdaptiveQuadrature.hs b/src/Physics/Learn/AdaptiveQuadrature.hs deleted file mode 100644 --- a/src/Physics/Learn/AdaptiveQuadrature.hs +++ /dev/null @@ -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) - diff --git a/src/Physics/Learn/Visual/VisTools.hs b/src/Physics/Learn/Visual/VisTools.hs --- a/src/Physics/Learn/Visual/VisTools.hs +++ b/src/Physics/Learn/Visual/VisTools.hs @@ -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 diff --git a/src/Tests.hs b/src/Tests.hs deleted file mode 100644 --- a/src/Tests.hs +++ /dev/null @@ -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 -