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.2 +Version: 0.3 Synopsis: Haskell code for learning physics Description: A library of functions for vector calculus, calculation of electric field, electric flux, @@ -29,6 +29,11 @@ Physics.Learn.RungeKutta Physics.Learn.CompositeQuadrature Physics.Learn.RootFinding + Physics.Learn.Mechanics + Physics.Learn + Physics.Learn.Visual.VisTools Build-depends: base >= 4.2 && < 4.8, - vector-space >= 0.8.4 && < 0.9 + vector-space >= 0.8.4 && < 0.9, + not-gloss >= 0.5.0.4 && < 0.7, + spatial-math >= 0.1.7 && < 0.3 Hs-source-dirs: src diff --git a/src/Physics/Learn.hs b/src/Physics/Learn.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn.hs @@ -0,0 +1,307 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn +Copyright : (c) Scott N. Walck 2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +Functions for learning physics. +-} + +module Physics.Learn + ( + -- * Mechanics + TheTime + , TimeStep + , Velocity + -- ** Simple one-particle state + , SimpleState + , SimpleAccelerationFunction + , simpleStateDeriv + , simpleRungeKuttaStep + -- ** One-particle state + , St(..) + , DSt(..) + , OneParticleSystemState + , OneParticleAccelerationFunction + , oneParticleStateDeriv + , oneParticleRungeKuttaStep + , oneParticleRungeKuttaSolution + -- ** Two-particle state + , TwoParticleSystemState + , TwoParticleAccelerationFunction + , twoParticleStateDeriv + , twoParticleRungeKuttaStep + -- ** Many-particle state + , ManyParticleSystemState + , ManyParticleAccelerationFunction + , manyParticleStateDeriv + , manyParticleRungeKuttaStep + -- * E&M + -- ** Charge + , Charge + , ChargeDistribution(..) + , totalCharge + -- ** Current + , Current + , CurrentDistribution(..) + -- ** Electric Field + , eField + -- ** Electric Flux + , electricFlux + -- ** Electric Potential + , electricPotentialFromField + , electricPotentialFromCharge + -- ** Magnetic Field + , bField + -- ** Magnetic Flux + , magneticFlux + -- * Geometry + -- ** Vectors + , Vec + , xComp + , yComp + , zComp + , vec + , (^+^) + , (^-^) + , (*^) + , (^*) + , (^/) + , (<.>) + , (><) + , magnitude + , zeroV + , negateV + , sumV + , iHat + , jHat + , kHat + -- ** 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 + -- ** Curves + , Curve(..) + , normalizeCurve + , concatCurves + , concatenateCurves + , reverseCurve + , evalCurve + , shiftCurve + , straightLine + -- ** Line Integrals + , simpleLineIntegral + , dottedLineIntegral + , crossedLineIntegral + , compositeSimpsonDottedLineIntegral + , compositeSimpsonCrossedLineIntegral + -- ** Surfaces + , Surface(..) + , unitSphere + , centeredSphere + , sphere + , northernHemisphere + , disk + , shiftSurface + -- ** Surface Integrals + , surfaceIntegral + , dottedSurfaceIntegral + -- ** Volumes + , Volume(..) + , unitBall + , unitBallCartesian + , centeredBall + , ball + , northernHalfBall + , centeredCylinder + , shiftVolume + -- ** Volume Integral + , volumeIntegral + -- * Utilities + , StateSpace(..) +-- , (.-^) +-- , Time + , rungeKutta4 + , integrateSystem + -- * Visualization + -- ** Vis library + , xyzFromVec + , xyzFromPos + , visVec + , oneVector + , displayVectorField + , curveObject + ) + where + +import Physics.Learn.Charge + ( Charge + , ChargeDistribution(..) + , totalCharge + , eField + , electricFlux + , electricPotentialFromField + , electricPotentialFromCharge + ) +import Physics.Learn.Current + ( Current + , CurrentDistribution(..) + , bField + , magneticFlux + ) +import Physics.Learn.CarrotVec + ( Vec + , xComp + , yComp + , zComp + , vec + , (^+^) + , (^-^) + , (*^) + , (^*) + , (^/) + , (<.>) + , (><) + , magnitude + , zeroV + , negateV + , sumV + , iHat + , jHat + , kHat + ) +import 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 + ) +import Physics.Learn.Curve + ( Curve(..) + , normalizeCurve + , concatCurves + , concatenateCurves + , reverseCurve + , evalCurve + , shiftCurve + , straightLine + , simpleLineIntegral + , dottedLineIntegral + , crossedLineIntegral + , compositeSimpsonDottedLineIntegral + , compositeSimpsonCrossedLineIntegral + ) +import Physics.Learn.Surface + ( Surface(..) + , unitSphere + , centeredSphere + , sphere + , northernHemisphere + , disk + , shiftSurface + , surfaceIntegral + , dottedSurfaceIntegral + ) +import Physics.Learn.Volume + ( Volume(..) + , unitBall + , unitBallCartesian + , centeredBall + , ball + , northernHalfBall + , centeredCylinder + , shiftVolume + , volumeIntegral + ) +import Physics.Learn.Visual.VisTools + ( xyzFromVec + , xyzFromPos + , visVec + , oneVector + , displayVectorField + , curveObject + ) +import Physics.Learn.StateSpace + ( StateSpace(..) +-- , (.-^) +-- , Time + ) +import Physics.Learn.RungeKutta + ( rungeKutta4 + , integrateSystem + ) +import Physics.Learn.Mechanics + ( TheTime + , TimeStep + , Velocity + , SimpleState + , SimpleAccelerationFunction + , simpleStateDeriv + , simpleRungeKuttaStep + , St(..) + , DSt(..) + , OneParticleSystemState + , OneParticleAccelerationFunction + , oneParticleStateDeriv + , oneParticleRungeKuttaStep + , oneParticleRungeKuttaSolution + , TwoParticleSystemState + , TwoParticleAccelerationFunction + , twoParticleStateDeriv + , twoParticleRungeKuttaStep + , ManyParticleSystemState + , ManyParticleAccelerationFunction + , manyParticleStateDeriv + , manyParticleRungeKuttaStep + ) diff --git a/src/Physics/Learn/Charge.hs b/src/Physics/Learn/Charge.hs --- a/src/Physics/Learn/Charge.hs +++ b/src/Physics/Learn/Charge.hs @@ -17,6 +17,7 @@ -- * Charge Charge , ChargeDistribution(..) + , totalCharge -- * Electric Field , eField , eFieldFromPointCharge @@ -59,7 +60,7 @@ , volumeIntegral ) --- | 'Charge' is just a synonym for a double-precision floating point number. +-- | Electric charge, in units of Coulombs (C) type Charge = Double -- | A charge distribution is a point charge, a line charge, a surface charge, @@ -67,11 +68,19 @@ -- 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 + | LineCharge ScalarField Curve -- ^ 'ScalarField' is linear charge density (C/m) + | SurfaceCharge ScalarField Surface -- ^ 'ScalarField' is surface charge density (C/m^2) + | VolumeCharge ScalarField Volume -- ^ 'ScalarField' is volume charge density (C/m^3) + | MultipleCharges [ChargeDistribution] -- ^ combination of charge distributions +-- | Total charge (in C) of a charge distribution. +totalCharge :: ChargeDistribution -> Charge +totalCharge (PointCharge q _) = q +totalCharge (LineCharge lambda c) = simpleLineIntegral 1000 lambda c +totalCharge (SurfaceCharge sigma s) = surfaceIntegral 100 100 sigma s +totalCharge (VolumeCharge rho v) = volumeIntegral 50 50 50 rho v +totalCharge (MultipleCharges ds) = sum [totalCharge d | d <- ds] + {- shiftChargeDistribution :: Displacement -> ChargeDistribution -> ChargeDistribution shiftChargeDistribution d (Point @@ -83,7 +92,7 @@ eFieldFromPointCharge :: Charge -- ^ charge (in Coulombs) -> Position -- ^ of point charge - -> VectorField -- ^ electric field + -> VectorField -- ^ electric field (in V/m) eFieldFromPointCharge q r' r = (k * q) *^ d ^/ magnitude d ** 3 where @@ -96,7 +105,7 @@ eFieldFromLineCharge :: ScalarField -- ^ linear charge density lambda -> Curve -- ^ geometry of the line charge - -> VectorField -- ^ electric field + -> VectorField -- ^ electric field (in V/m) eFieldFromLineCharge lambda c r = k *^ simpleLineIntegral 1000 integrand c where @@ -111,7 +120,7 @@ eFieldFromSurfaceCharge :: ScalarField -- ^ surface charge density sigma -> Surface -- ^ geometry of the surface charge - -> VectorField -- ^ electric field + -> VectorField -- ^ electric field (in V/m) eFieldFromSurfaceCharge sigma s r = k *^ surfaceIntegral 100 100 integrand s where @@ -126,7 +135,7 @@ eFieldFromVolumeCharge :: ScalarField -- ^ volume charge density rho -> Volume -- ^ geometry of the volume charge - -> VectorField -- ^ electric field + -> VectorField -- ^ electric field (in V/m) eFieldFromVolumeCharge rho v r = k *^ volumeIntegral 50 50 50 integrand v where @@ -143,7 +152,7 @@ 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 +eField (MultipleCharges cds) = addFields $ map eField cds ------------------- -- Electric Flux -- @@ -171,7 +180,7 @@ 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 +electricPotentialFromCharge (MultipleCharges cds) = addFields $ map electricPotentialFromCharge cds ePotFromPointCharge :: Charge -- ^ charge (in Coulombs) diff --git a/src/Physics/Learn/Current.hs b/src/Physics/Learn/Current.hs --- a/src/Physics/Learn/Current.hs +++ b/src/Physics/Learn/Current.hs @@ -52,7 +52,7 @@ , volumeIntegral ) --- | 'Current' is just a synonym for a double-precision floating point number. +-- | Electric current, in units of Amperes (A) type Current = Double -- | A current distribution is a line current (current through a wire), a surface current, @@ -60,8 +60,8 @@ -- 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 + | SurfaceCurrent VectorField Surface -- ^ 'VectorField' is surface current density (A/m) + | VolumeCurrent VectorField Volume -- ^ 'VectorField' is volume current density (A/m^2) | MultipleCurrents [CurrentDistribution] -- ^ combination of current distributions -- | Magnetic field produced by a line current (current through a wire). @@ -70,7 +70,7 @@ bFieldFromLineCurrent :: Current -- ^ current (in Amps) -> Curve -- ^ geometry of the line current - -> VectorField -- ^ magnetic field + -> VectorField -- ^ magnetic field (in Tesla) bFieldFromLineCurrent i c r = k *^ crossedLineIntegral 1000 integrand c where @@ -88,7 +88,7 @@ bFieldFromSurfaceCurrent :: VectorField -- ^ surface current density -> Surface -- ^ geometry of the surface current - -> VectorField -- ^ magnetic field + -> VectorField -- ^ magnetic field (in T) bFieldFromSurfaceCurrent kCurrent c r = k *^ surfaceIntegral 100 100 integrand c where @@ -103,7 +103,7 @@ bFieldFromVolumeCurrent :: VectorField -- ^ volume current density -> Volume -- ^ geometry of the volume current - -> VectorField -- ^ magnetic field + -> VectorField -- ^ magnetic field (in T) bFieldFromVolumeCurrent j c r = k *^ volumeIntegral 50 50 50 integrand c where diff --git a/src/Physics/Learn/Mechanics.hs b/src/Physics/Learn/Mechanics.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/Mechanics.hs @@ -0,0 +1,225 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.Mechanics +Copyright : (c) Scott N. Walck 2014 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +Newton's second law and all that +-} + +module Physics.Learn.Mechanics + ( TheTime + , TimeStep + , Velocity + -- * Simple one-particle state + , SimpleState + , SimpleAccelerationFunction + , simpleStateDeriv + , simpleRungeKuttaStep + -- * One-particle state + , St(..) + , DSt(..) + , OneParticleSystemState + , OneParticleAccelerationFunction + , oneParticleStateDeriv + , oneParticleRungeKuttaStep + , oneParticleRungeKuttaSolution + -- * Two-particle state + , TwoParticleSystemState + , TwoParticleAccelerationFunction + , twoParticleStateDeriv + , twoParticleRungeKuttaStep + -- * Many-particle state + , ManyParticleSystemState + , ManyParticleAccelerationFunction + , manyParticleStateDeriv + , manyParticleRungeKuttaStep + ) + where + +import Data.VectorSpace + ( AdditiveGroup(..) + , VectorSpace(..) + ) +import Physics.Learn.StateSpace + ( StateSpace(..) + , Diff + , TimeDerivative + ) +import Physics.Learn.RungeKutta + ( rungeKutta4 + , integrateSystem + ) +import Physics.Learn.Position + ( Position + ) +import Physics.Learn.CarrotVec + ( Vec + ) + +-- | Time (in s). +type TheTime = Double + +-- | A time step (in s). +type TimeStep = Double + +-- | Velocity of a particle (in m/s). +type Velocity = Vec + +------------------------------- +-- Simple one-particle state -- +------------------------------- + +-- | A simple one-particle state, +-- to get started quickly with mechanics of one particle. +type SimpleState = (TheTime,Position,Velocity) + +-- | An acceleration function gives the particle's acceleration as +-- a function of the particle's state. +-- The specification of this function is what makes one single-particle +-- mechanics problem different from another. +-- In order to write this function, add all of the forces +-- that act on the particle, and divide this net force by the particle's mass. +-- (Newton's second law). +type SimpleAccelerationFunction = SimpleState -> Vec + +-- | Time derivative of state for a single particle +-- with a constant mass. +simpleStateDeriv :: SimpleAccelerationFunction -- ^ acceleration function for the particle + -> TimeDerivative SimpleState -- ^ derivatives as a function of state +simpleStateDeriv a (t, r, v) = (1, v, a(t, r, v)) + +-- | Single Runge-Kutta step +simpleRungeKuttaStep :: SimpleAccelerationFunction -- ^ acceleration function for the particle + -> TimeStep -- ^ time step + -> SimpleState -- ^ initial state + -> SimpleState -- ^ state after one time step +simpleRungeKuttaStep = rungeKutta4 . simpleStateDeriv + +------------------------ +-- One-particle state -- +------------------------ + +-- | The state of a single particle is given by +-- the position of the particle and the velocity of the particle. +data St = St { position :: Position + , velocity :: Velocity + } + deriving (Show) + +-- | The associated vector space for the +-- state of a single particle. +data DSt = DSt Vec Vec + deriving (Show) + +instance AdditiveGroup DSt where + zeroV = DSt zeroV zeroV + negateV (DSt dr dv) = DSt (negateV dr) (negateV dv) + DSt dr1 dv1 ^+^ DSt dr2 dv2 = DSt (dr1 ^+^ dr2) (dv1 ^+^ dv2) + +instance VectorSpace DSt where + type Scalar DSt = Double + c *^ DSt dr dv = DSt (c*^dr) (c*^dv) + +instance StateSpace St where + type Diff St = DSt + St r1 v1 .-. St r2 v2 = DSt (r1 .-. r2) (v1 .-. v2) + St r1 v1 .+^ DSt dr dv = St (r1 .+^ dr) (v1 .+^ dv) + +-- | The state of a system of one particle is given by the current time, +-- the position of the particle, and the velocity of the particle. +-- Including time in the state like this allows us to +-- have time-dependent forces. +type OneParticleSystemState = (TheTime,St) + +-- | An acceleration function gives the particle's acceleration as +-- a function of the particle's state. +type OneParticleAccelerationFunction = OneParticleSystemState -> Vec + +-- | Time derivative of state for a single particle +-- with a constant mass. +oneParticleStateDeriv :: OneParticleAccelerationFunction -- ^ acceleration function for the particle + -> TimeDerivative OneParticleSystemState -- ^ derivatives as a function of state +oneParticleStateDeriv a st@(_t, St _r v) = (1, DSt v (a st)) + +-- | Single Runge-Kutta step +oneParticleRungeKuttaStep :: OneParticleAccelerationFunction -- ^ acceleration function for the particle + -> TimeStep -- ^ time step + -> OneParticleSystemState -- ^ initial state + -> OneParticleSystemState -- ^ state after one time step +oneParticleRungeKuttaStep = rungeKutta4 . oneParticleStateDeriv + +-- | List of system states +oneParticleRungeKuttaSolution :: OneParticleAccelerationFunction -- ^ acceleration function for the particle + -> TimeStep -- ^ time step + -> OneParticleSystemState -- ^ initial state + -> [OneParticleSystemState] -- ^ state after one time step +oneParticleRungeKuttaSolution = integrateSystem . oneParticleStateDeriv + +------------------------ +-- Two-particle state -- +------------------------ + +-- | The state of a system of two particles is given by the current time, +-- the position and velocity of particle 1, +-- and the position and velocity of particle 2. +type TwoParticleSystemState = (TheTime,St,St) + +-- | An acceleration function gives a pair of accelerations +-- (one for particle 1, one for particle 2) as +-- a function of the system's state. +type TwoParticleAccelerationFunction = TwoParticleSystemState -> (Vec,Vec) + +-- | Time derivative of state for two particles +-- with constant mass. +twoParticleStateDeriv :: TwoParticleAccelerationFunction -- ^ acceleration function for two particles + -> TimeDerivative TwoParticleSystemState -- ^ derivatives as a function of state +twoParticleStateDeriv af2 st2@(_t, St _r1 v1, St _r2 v2) = (1, DSt v1 a1, DSt v2 a2) + where + (a1,a2) = af2 st2 + +-- | Single Runge-Kutta step for two-particle system +twoParticleRungeKuttaStep :: TwoParticleAccelerationFunction -- ^ acceleration function + -> TimeStep -- ^ time step + -> TwoParticleSystemState -- ^ initial state + -> TwoParticleSystemState -- ^ state after one time step +twoParticleRungeKuttaStep = rungeKutta4 . twoParticleStateDeriv + +------------------------- +-- Many-particle state -- +------------------------- + +-- | The state of a system of many particles is given by the current time +-- and a list of one-particle states. +type ManyParticleSystemState = (TheTime,[St]) + +-- | An acceleration function gives a list of accelerations +-- (one for each particle) as +-- a function of the system's state. +type ManyParticleAccelerationFunction = ManyParticleSystemState -> [Vec] + +-- | Time derivative of state for many particles +-- with constant mass. +manyParticleStateDeriv :: ManyParticleAccelerationFunction -- ^ acceleration function for many particles + -> TimeDerivative ManyParticleSystemState -- ^ derivatives as a function of state +manyParticleStateDeriv af st@(_t, sts) = (1, [DSt v a | (v,a) <- zip vs as]) + where + vs = map velocity sts + as = af st + +-- | Single Runge-Kutta step for many-particle system +manyParticleRungeKuttaStep :: ManyParticleAccelerationFunction -- ^ acceleration function + -> TimeStep -- ^ time step + -> ManyParticleSystemState -- ^ initial state + -> ManyParticleSystemState -- ^ state after one time step +manyParticleRungeKuttaStep = rungeKutta4 . manyParticleStateDeriv + + + +-- Can we automatically incorporate Newton's third law? + diff --git a/src/Physics/Learn/Position.hs b/src/Physics/Learn/Position.hs --- a/src/Physics/Learn/Position.hs +++ b/src/Physics/Learn/Position.hs @@ -63,6 +63,7 @@ -- | A type for position. -- Position is not a vector because it makes no sense to add positions. data Position = Cart Double Double Double + deriving (Show) -- | A displacement is a vector. type Displacement = Vec diff --git a/src/Physics/Learn/StateSpace.hs b/src/Physics/Learn/StateSpace.hs --- a/src/Physics/Learn/StateSpace.hs +++ b/src/Physics/Learn/StateSpace.hs @@ -1,4 +1,4 @@ -{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_GHC -Wall -fno-warn-orphans #-} {-# LANGUAGE FlexibleContexts, FlexibleInstances, TypeFamilies #-} {-# LANGUAGE Trustworthy #-} @@ -22,13 +22,16 @@ ( StateSpace(..) , (.-^) , Time + , TimeDerivative ) where +import Data.AdditiveGroup + ( AdditiveGroup(..) + ) import Data.VectorSpace - ( VectorSpace - , Scalar - , negateV + ( VectorSpace(..) +-- , Scalar ) import Physics.Learn.Position ( Position @@ -37,7 +40,7 @@ ) import Physics.Learn.CarrotVec ( Vec - , (^+^) +-- , (^+^) , (^-^) ) @@ -88,3 +91,28 @@ 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) + +inf :: a -> [a] +inf x = x : inf x + +instance AdditiveGroup v => AdditiveGroup [v] where + zeroV = inf zeroV + (^+^) = zipWith (^+^) + negateV = map negateV + +instance VectorSpace v => VectorSpace [v] where + type Scalar [v] = Scalar v + c *^ xs = [c *^ x | x <- xs] + +instance StateSpace p => StateSpace [p] where + type Diff [p] = [Diff p] + (.-.) = zipWith (.-.) + (.+^) = zipWith (.+^) + +-- | The time derivative of a state is an element of the associated vector space. +type TimeDerivative state = state -> Diff state + +{- +class HasTimeDerivative state where + timeDeriv :: state -> Diff state +-} diff --git a/src/Physics/Learn/Surface.hs b/src/Physics/Learn/Surface.hs --- a/src/Physics/Learn/Surface.hs +++ b/src/Physics/Learn/Surface.hs @@ -14,13 +14,16 @@ -} module Physics.Learn.Surface - ( Surface(..) + ( + -- * Surfaces + Surface(..) , unitSphere , centeredSphere , sphere , northernHemisphere , disk , shiftSurface + -- * Surface Integrals , surfaceIntegral , dottedSurfaceIntegral ) @@ -87,6 +90,8 @@ -- | 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)) + +-- To do : boundaryOfSurface :: Surface -> Curve -- | A plane surface integral, in which area element is a scalar. surfaceIntegral :: (VectorSpace v, Scalar v ~ Double) => diff --git a/src/Physics/Learn/Visual/VisTools.hs b/src/Physics/Learn/Visual/VisTools.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/Visual/VisTools.hs @@ -0,0 +1,126 @@ +{-# OPTIONS_GHC -Wall #-} + +-- | Some tools related to the not-gloss 3D graphics and animation library. + +module Physics.Learn.Visual.VisTools + ( xyzFromVec + , xyzFromPos + , visVec + , oneVector + , displayVectorField + , curveObject + ) + where + +import SpatialMath + ( Xyz(..) + , Euler(..) + ) +import Vis + ( VisObject(..) + , Color + ) +import Physics.Learn.CarrotVec + ( Vec + , xComp + , yComp + , zComp +-- , magnitude + , (^/) + ) +import Physics.Learn.Position + ( Position + , cartesianCoordinates + , VectorField + ) +import Physics.Learn.Curve + ( Curve(..) + ) + +-- | Make an 'Xyz' object from a 'Vec'. +xyzFromVec :: Vec -> Xyz Double +xyzFromVec v = Xyz x y z + where + x = xComp v + y = yComp v + z = zComp v + +-- | Make an 'Xyz' object from a 'Position'. +xyzFromPos :: Position -> Xyz Double +xyzFromPos r = Xyz x y z + where + (x,y,z) = cartesianCoordinates r + +-- | Display a vector field. +displayVectorField :: Color -- ^ color for the vector field + -> Double -- ^ scale factor + -> [Position] -- ^ list of positions to show the field + -> VectorField -- ^ vector field to display + -> VisObject Double -- ^ the displayable object +displayVectorField col unitsPerMeter samplePts field + = VisObjects [Trans (xyzFromPos r) $ visVec col (e ^/ unitsPerMeter) | r <- samplePts, let e = field r] + +-- | A displayable VisObject for a curve. +curveObject :: Color -> Curve -> VisObject Double +curveObject color (Curve f a b) + = Line' [(xyzFromPos (f t), color) | t <- [a,a+(b-a)/1000..b]] + +-- | Place a vector at a particular position. +oneVector :: Color -> Position -> Vec -> VisObject Double +oneVector c r v = Trans (xyzFromPos r) $ visVec c v + +data Cart = Cart Double Double Double + deriving (Show) + +data Sph = Sph Double Double Double + deriving (Show) + +sphericalCoords :: Cart -> Sph +sphericalCoords (Cart x y z) = Sph r theta phi + where + r = sqrt (x*x + y*y + z*z) + s = sqrt (x*x + y*y) + theta = atan2 s z + phi = atan2 y x + +-- | A VisObject arrow from a vector +visVec :: Color -> Vec -> VisObject Double +visVec c v = rotZ phi $ rotY theta $ Arrow (r,20*r) (Xyz 0 0 1) c + where + x = xComp v + y = yComp v + z = zComp v + Sph r theta phi = sphericalCoords (Cart x y z) + +{- +rotX :: Double -- ^ in radians + -> VisObject Double + -> VisObject Double +rotX alpha = RotEulerRad (Euler 0 0 alpha) +-} + +rotY :: Double -- ^ in radians + -> VisObject Double + -> VisObject Double +rotY alpha = RotEulerRad (Euler 0 alpha 0) + +rotZ :: Double -- ^ in radians + -> VisObject Double + -> VisObject Double +rotZ alpha = RotEulerRad (Euler alpha 0 0) + + +{- +adjacentDistance :: [Position] -> Double +adjacentDistance [] = 0 +adjacentDistance rs'@(_:rs) = minimum (map magnitude $ zipWith displacement rs' rs) + +visVectorField :: Color -> [Position] -> VectorField -> VisObject Double +visVectorField c rs vf = let prs = [(r,vf r) | r <- rs] + bigV = maximum [magnitude (snd pr) | pr <- prs] + disp = adjacentDistance rs + scaleFactor = disp / bigV + newPrs = [(r, scaleFactor *^ v) | (r,v) <- prs] + vecs = [oneVector c r v' | (r,v') <- newPrs] + in VisObjects vecs +-} diff --git a/src/Physics/Learn/Volume.hs b/src/Physics/Learn/Volume.hs --- a/src/Physics/Learn/Volume.hs +++ b/src/Physics/Learn/Volume.hs @@ -14,7 +14,9 @@ -} module Physics.Learn.Volume - ( Volume(..) + ( + -- * Volumes + Volume(..) , unitBall , unitBallCartesian , centeredBall @@ -22,6 +24,7 @@ , northernHalfBall , centeredCylinder , shiftVolume + -- * Volume Integral , volumeIntegral ) where diff --git a/src/examples/BCircularLoop.hs b/src/examples/BCircularLoop.hs new file mode 100644 --- /dev/null +++ b/src/examples/BCircularLoop.hs @@ -0,0 +1,27 @@ +{-# 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] + +main :: IO () +main = display Nothing "Magnetic Field from a Current Loop" drawFun diff --git a/src/examples/LorentzForceSimulation.hs b/src/examples/LorentzForceSimulation.hs new file mode 100644 --- /dev/null +++ b/src/examples/LorentzForceSimulation.hs @@ -0,0 +1,64 @@ +{-# OPTIONS_GHC -Wall #-} + +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 (xyzFromPos 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) + +main :: IO () +main = simulate + Nothing + "Particle Experiencing Electromagnetic Force" + 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/src/examples/sunEarthRK4.hs b/src/examples/sunEarthRK4.hs new file mode 100644 --- /dev/null +++ b/src/examples/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 (2 * earthSunDistance) 0 0 + ,vec 0 (initialEarthSpeed / 2) 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" (600, 600) (10, 10)) + black 50 initialState worldToPicture simStep