{-# LANGUAGE TypeFamilies, FlexibleContexts #-}
{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE Safe #-}

{- | 
Module      :  Physics.Learn.Curve
Copyright   :  (c) Scott N. Walck 2012-2018
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
    , compositeTrapezoidDottedLineIntegral
    , compositeTrapezoidCrossedLineIntegral
    , 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 { Curve -> R -> Position
curveFunc          :: Double -> Position  -- ^ function from one parameter into space
                   , Curve -> R
startingCurveParam :: Double              -- ^ starting value of the parameter
                   , Curve -> R
endingCurveParam   :: Double              -- ^ ending value of the parameter
                   }

-- | A dotted line integral.
--   Convenience function for 'compositeSimpsonDottedLineIntegral'.
dottedLineIntegral
    :: Int          -- ^ number of half-intervals
                    --   (one less than the number of function evaluations)
    -> VectorField  -- ^ vector field
    -> Curve        -- ^ curve to integrate over
    -> Double       -- ^ scalar result
dottedLineIntegral :: Int -> VectorField -> Curve -> R
dottedLineIntegral = Int -> VectorField -> Curve -> R
compositeSimpsonDottedLineIntegral

-- | Calculates integral vf x dl over curve.
--   Convenience function for 'compositeSimpsonCrossedLineIntegral'.
crossedLineIntegral
    :: Int          -- ^ number of half-intervals
                    --   (one less than the number of function evaluations)
    -> VectorField  -- ^ vector field
    -> Curve        -- ^ curve to integrate over
    -> Vec          -- ^ vector result
crossedLineIntegral :: Int -> VectorField -> Curve -> Vec
crossedLineIntegral = Int -> VectorField -> Curve -> Vec
compositeSimpsonCrossedLineIntegral

-- | A dotted line integral, performed in an unsophisticated way.
compositeTrapezoidDottedLineIntegral
    :: Int          -- ^ number of intervals
    -> VectorField  -- ^ vector field
    -> Curve        -- ^ curve to integrate over
    -> Double       -- ^ scalar result
compositeTrapezoidDottedLineIntegral :: Int -> VectorField -> Curve -> R
compositeTrapezoidDottedLineIntegral Int
n VectorField
vf (Curve R -> Position
f R
a R
b)
    = [R] -> R
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([R] -> R) -> [R] -> R
forall a b. (a -> b) -> a -> b
$ (Vec -> Vec -> R) -> [Vec] -> [Vec] -> [R]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Vec -> Vec -> R
Vec -> Vec -> Scalar Vec
forall v. InnerSpace v => v -> v -> Scalar v
(<.>) [Vec]
aveVecs [Vec]
dls
      where
        dt :: R
dt = (R
b R -> R -> R
forall a. Num a => a -> a -> a
- R
a) R -> R -> R
forall a. Fractional a => a -> a -> a
/ Int -> R
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
        pts :: [Position]
pts = [R -> Position
f R
t | R
t <- [R
a,R
aR -> R -> R
forall a. Num a => a -> a -> a
+R
dt..R
b]]
        vecs :: [Vec]
vecs = [VectorField
vf Position
pt | Position
pt <- [Position]
pts]
        aveVecs :: [Vec]
aveVecs = (Vec -> Vec -> Vec) -> [Vec] -> [Vec] -> [Vec]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Vec -> Vec -> Vec
forall v. (VectorSpace v, Scalar v ~ R) => v -> v -> v
average [Vec]
vecs ([Vec] -> [Vec]
forall a. HasCallStack => [a] -> [a]
tail [Vec]
vecs)
        dls :: [Vec]
dls = (Position -> VectorField) -> [Position] -> [Position] -> [Vec]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Position -> VectorField
displacement [Position]
pts ([Position] -> [Position]
forall a. HasCallStack => [a] -> [a]
tail [Position]
pts)

-- | Calculates integral vf x dl over curve in an unsophisticated way.
compositeTrapezoidCrossedLineIntegral
    :: Int          -- ^ number of intervals
    -> VectorField  -- ^ vector field
    -> Curve        -- ^ curve to integrate over
    -> Vec          -- ^ vector result
compositeTrapezoidCrossedLineIntegral :: Int -> VectorField -> Curve -> Vec
compositeTrapezoidCrossedLineIntegral Int
n VectorField
vf (Curve R -> Position
f R
a R
b)
    = [Vec] -> Vec
forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV ([Vec] -> Vec) -> [Vec] -> Vec
forall a b. (a -> b) -> a -> b
$ (Vec -> Vec -> Vec) -> [Vec] -> [Vec] -> [Vec]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Vec -> Vec -> Vec
(><) [Vec]
aveVecs [Vec]
dls
      where
        dt :: R
dt = (R
b R -> R -> R
forall a. Num a => a -> a -> a
- R
a) R -> R -> R
forall a. Fractional a => a -> a -> a
/ Int -> R
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
        pts :: [Position]
pts = [R -> Position
f R
t | R
t <- [R
a,R
aR -> R -> R
forall a. Num a => a -> a -> a
+R
dt..R
b]]
        vecs :: [Vec]
vecs = [VectorField
vf Position
pt | Position
pt <- [Position]
pts]
        aveVecs :: [Vec]
aveVecs = (Vec -> Vec -> Vec) -> [Vec] -> [Vec] -> [Vec]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Vec -> Vec -> Vec
forall v. (VectorSpace v, Scalar v ~ R) => v -> v -> v
average [Vec]
vecs ([Vec] -> [Vec]
forall a. HasCallStack => [a] -> [a]
tail [Vec]
vecs)
        dls :: [Vec]
dls = (Position -> VectorField) -> [Position] -> [Position] -> [Vec]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Position -> VectorField
displacement [Position]
pts ([Position] -> [Position]
forall a. HasCallStack => [a] -> [a]
tail [Position]
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 :: forall v.
(InnerSpace v, Scalar v ~ R) =>
Int -> Field v -> Curve -> v
simpleLineIntegral Int
n Field v
vf (Curve R -> Position
f R
a R
b)
    = [v] -> v
forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV ([v] -> v) -> [v] -> v
forall a b. (a -> b) -> a -> b
$ (v -> R -> v) -> [v] -> [R] -> [v]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith v -> R -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
(^*) [v]
aveVecs ((Vec -> R) -> [Vec] -> [R]
forall a b. (a -> b) -> [a] -> [b]
map Vec -> R
forall v s. (InnerSpace v, s ~ Scalar v, Floating s) => v -> s
magnitude [Vec]
dls)
      where
        dt :: R
dt = (R
b R -> R -> R
forall a. Num a => a -> a -> a
- R
a) R -> R -> R
forall a. Fractional a => a -> a -> a
/ Int -> R
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
        pts :: [Position]
pts = [R -> Position
f R
t | R
t <- [R
a,R
aR -> R -> R
forall a. Num a => a -> a -> a
+R
dt..R
b]]
        vecs :: [v]
vecs = [Field v
vf Position
pt | Position
pt <- [Position]
pts]
        aveVecs :: [v]
aveVecs = (v -> v -> v) -> [v] -> [v] -> [v]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith v -> v -> v
forall v. (VectorSpace v, Scalar v ~ R) => v -> v -> v
average [v]
vecs ([v] -> [v]
forall a. HasCallStack => [a] -> [a]
tail [v]
vecs)
        dls :: [Vec]
dls = (Position -> VectorField) -> [Position] -> [Position] -> [Vec]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Position -> VectorField
displacement [Position]
pts ([Position] -> [Position]
forall a. HasCallStack => [a] -> [a]
tail [Position]
pts)

-- | Reparametrize a curve from 0 to 1.
normalizeCurve :: Curve -> Curve
normalizeCurve :: Curve -> Curve
normalizeCurve (Curve R -> Position
f R
a R
b)
    = (R -> Position) -> R -> R -> Curve
Curve (R -> Position
f (R -> Position) -> (R -> R) -> R -> Position
forall b c a. (b -> c) -> (a -> b) -> a -> c
. R -> R
scl) R
0 R
1
      where
        scl :: R -> R
scl R
t = R
a R -> R -> R
forall a. Num a => a -> a -> a
+ (R
b R -> R -> R
forall a. Num a => a -> a -> a
- R
a) R -> R -> R
forall a. Num a => a -> a -> a
* R
t

-- | Concatenate two curves.
concatCurves :: Curve  -- ^ go first along this curve
             -> Curve  -- ^ then along this curve
             -> Curve  -- ^ to produce this new curve
concatCurves :: Curve -> Curve -> Curve
concatCurves Curve
c1 Curve
c2
    = Curve -> Curve
normalizeCurve (Curve -> Curve) -> Curve -> Curve
forall a b. (a -> b) -> a -> b
$ (R -> Position) -> R -> R -> Curve
Curve R -> Position
f R
0 R
2
      where
        (Curve R -> Position
f1 R
_ R
_) = Curve -> Curve
normalizeCurve Curve
c1
        (Curve R -> Position
f2 R
_ R
_) = Curve -> Curve
normalizeCurve Curve
c2
        f :: R -> Position
f R
t | R
t R -> R -> Bool
forall a. Ord a => a -> a -> Bool
<= R
1     = R -> Position
f1 R
t
            | Bool
otherwise  = R -> Position
f2 (R
tR -> R -> R
forall a. Num a => a -> a -> a
-R
1)

-- | Concatenate a list of curves.
--   Parametrizes curves equally.
concatenateCurves :: [Curve] -> Curve
concatenateCurves :: [Curve] -> Curve
concatenateCurves []     = [Char] -> Curve
forall a. HasCallStack => [Char] -> a
error [Char]
"concatenateCurves:  cannot concatenate empty list"
concatenateCurves [Curve]
cs = Curve -> Curve
normalizeCurve (Curve -> Curve) -> Curve -> Curve
forall a b. (a -> b) -> a -> b
$ (R -> Position) -> R -> R -> Curve
Curve R -> Position
f R
0 (Int -> R
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
    where
      n :: Int
n   = [Curve] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Curve]
cs
      ncs :: [Curve]
ncs = (Curve -> Curve) -> [Curve] -> [Curve]
forall a b. (a -> b) -> [a] -> [b]
map Curve -> Curve
normalizeCurve [Curve]
cs
      f :: R -> Position
f R
t = Curve -> R -> Position
evalCurve ([Curve]
ncs [Curve] -> Int -> Curve
forall a. HasCallStack => [a] -> Int -> a
!! Int
m) (R
t R -> R -> R
forall a. Num a => a -> a -> a
- Int -> R
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
m)
          where m :: Int
m = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (R -> Int
forall b. Integral b => R -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor R
t)

-- | Reverse a curve.
reverseCurve :: Curve -> Curve
reverseCurve :: Curve -> Curve
reverseCurve (Curve R -> Position
f R
a R
b)
    = (R -> Position) -> R -> R -> Curve
Curve (R -> Position
f (R -> Position) -> (R -> R) -> R -> Position
forall b c a. (b -> c) -> (a -> b) -> a -> c
. R -> R
rev) R
a R
b
      where
        rev :: R -> R
rev R
t = R
a R -> R -> R
forall a. Num a => a -> a -> a
+ R
b R -> R -> R
forall a. Num a => a -> a -> a
- R
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 -> R -> Position
evalCurve (Curve R -> Position
f R
_ R
_) R
t = R -> Position
f R
t

-- | Shift a curve by a displacement.
shiftCurve :: Displacement  -- ^ amount to shift
           -> Curve         -- ^ original curve
           -> Curve         -- ^ shifted curve
shiftCurve :: Vec -> Curve -> Curve
shiftCurve Vec
d (Curve R -> Position
f R
sl R
su)
    = (R -> Position) -> R -> R -> Curve
Curve (Vec -> Position -> Position
shiftPosition Vec
d (Position -> Position) -> (R -> Position) -> R -> Position
forall b c a. (b -> c) -> (a -> b) -> a -> c
. R -> Position
f) R
sl R
su

-- | The straight-line curve from one position to another.
straightLine :: Position  -- ^ starting position
             -> Position  -- ^ ending position
             -> Curve     -- ^ straight-line curve
straightLine :: Position -> Position -> Curve
straightLine Position
r1 Position
r2 = (R -> Position) -> R -> R -> Curve
Curve R -> Position
f R
0 R
1
    where
      f :: R -> Position
f R
t = Vec -> Position -> Position
shiftPosition (R
Scalar Vec
t Scalar Vec -> Vec -> Vec
forall v. VectorSpace v => Scalar v -> v -> v
*^ Vec
d) Position
r1
      d :: Vec
d = Position -> VectorField
displacement Position
r1 Position
r2

-------------
-- Helpers --
-------------

average :: (VectorSpace v, Scalar v ~ Double) => v -> v -> v
average :: forall v. (VectorSpace v, Scalar v ~ R) => v -> v -> v
average v
v1 v
v2 = (v
v1 v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^+^ v
v2) v -> R -> v
forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/ R
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 :: forall v.
(InnerSpace v, Fractional (Scalar v)) =>
v -> v -> v -> v -> v -> Scalar v
dottedSimp v
f0 v
f1 v
f2 v
g10 v
g21
    = ((v
g21 v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^+^ v
g10) v -> Scalar v -> v
forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/ Scalar v
6) v -> v -> Scalar v
forall v. InnerSpace v => v -> v -> Scalar v
<.> (v
f0 v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^+^ Scalar v
4 Scalar v -> v -> v
forall v. VectorSpace v => Scalar v -> v -> v
*^ v
f1 v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^+^ v
f2)
      Scalar v -> Scalar v -> Scalar v
forall a. Num a => a -> a -> a
+ ((v
g21 v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^-^ v
g10) v -> Scalar v -> v
forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/ Scalar v
3) v -> v -> Scalar v
forall v. InnerSpace v => v -> v -> Scalar v
<.> (v
f2 v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^-^ v
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 :: Int -> VectorField -> Curve -> R
compositeSimpsonDottedLineIntegral Int
n VectorField
vf (Curve R -> Position
c R
a R
b)
    = let nEven :: Int
nEven = Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int -> Int -> Int
forall a. Integral a => a -> a -> a
div Int
n Int
2
          dt :: R
dt = (R
b R -> R -> R
forall a. Num a => a -> a -> a
- R
a) R -> R -> R
forall a. Fractional a => a -> a -> a
/ Int -> R
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
nEven
          ts :: [R]
ts = [R
a R -> R -> R
forall a. Num a => a -> a -> a
+ Int -> R
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
m R -> R -> R
forall a. Num a => a -> a -> a
* R
dt | Int
m <- [Int
0..Int
nEven]]
          pairs :: [(Position, Vec)]
pairs = [(Position
ct,VectorField
vf Position
ct) | R
t <- [R]
ts, let ct :: Position
ct = R -> Position
c R
t]
          combine :: [(Position, Vec)] -> R
combine [] = [Char] -> R
forall a. HasCallStack => [Char] -> a
error [Char]
"compositeSimpson: odd number of half-intervals" -- this should never happen
          combine [(Position, Vec)
_] = R
forall v. AdditiveGroup v => v
zeroV
          combine ((Position, Vec)
_:(Position, Vec)
_:[]) = [Char] -> R
forall a. HasCallStack => [Char] -> a
error [Char]
"compositeSimpson: odd number of half-intervals" -- this should never happen
          combine ((Position
c0,Vec
f0):(Position
c1,Vec
f1):(Position
c2,Vec
f2):[(Position, Vec)]
ps)
              = Vec -> Vec -> Vec -> Vec -> Vec -> Scalar Vec
forall v.
(InnerSpace v, Fractional (Scalar v)) =>
v -> v -> v -> v -> v -> Scalar v
dottedSimp Vec
f0 Vec
f1 Vec
f2 (Position -> VectorField
displacement Position
c0 Position
c1) (Position -> VectorField
displacement Position
c1 Position
c2)
                R -> R -> R
forall v. AdditiveGroup v => v -> v -> v
^+^ [(Position, Vec)] -> R
combine ((Position
c2,Vec
f2)(Position, Vec) -> [(Position, Vec)] -> [(Position, Vec)]
forall a. a -> [a] -> [a]
:[(Position, Vec)]
ps)
      in [(Position, Vec)] -> R
combine [(Position, Vec)]
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 :: Vec -> Vec -> Vec -> Vec -> Vec -> Vec
crossedSimp Vec
f0 Vec
f1 Vec
f2 Vec
g10 Vec
g21
    = Vec -> Vec
forall v. AdditiveGroup v => v -> v
negateV (Vec -> Vec) -> Vec -> Vec
forall a b. (a -> b) -> a -> b
$
      ((Vec
g21 Vec -> Vec -> Vec
forall v. AdditiveGroup v => v -> v -> v
^+^ Vec
g10) Vec -> R -> Vec
forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/ R
6) Vec -> Vec -> Vec
>< (Vec
f0 Vec -> Vec -> Vec
forall v. AdditiveGroup v => v -> v -> v
^+^ R
Scalar Vec
4 Scalar Vec -> Vec -> Vec
forall v. VectorSpace v => Scalar v -> v -> v
*^ Vec
f1 Vec -> Vec -> Vec
forall v. AdditiveGroup v => v -> v -> v
^+^ Vec
f2)
      Vec -> Vec -> Vec
forall v. AdditiveGroup v => v -> v -> v
^+^ ((Vec
g21 Vec -> Vec -> Vec
forall v. AdditiveGroup v => v -> v -> v
^-^ Vec
g10) Vec -> R -> Vec
forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/ R
3) Vec -> Vec -> Vec
>< (Vec
f2 Vec -> Vec -> Vec
forall v. AdditiveGroup v => v -> v -> v
^-^ Vec
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 :: Int -> VectorField -> Curve -> Vec
compositeSimpsonCrossedLineIntegral Int
n VectorField
vf (Curve R -> Position
c R
a R
b)
    = let nEven :: Int
nEven = Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int -> Int -> Int
forall a. Integral a => a -> a -> a
div Int
n Int
2
          dt :: R
dt = (R
b R -> R -> R
forall a. Num a => a -> a -> a
- R
a) R -> R -> R
forall a. Fractional a => a -> a -> a
/ Int -> R
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
nEven
          ts :: [R]
ts = [R
a R -> R -> R
forall a. Num a => a -> a -> a
+ Int -> R
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
m R -> R -> R
forall a. Num a => a -> a -> a
* R
dt | Int
m <- [Int
0..Int
nEven]]
          pairs :: [(Position, Vec)]
pairs = [(Position
ct,VectorField
vf Position
ct) | R
t <- [R]
ts, let ct :: Position
ct = R -> Position
c R
t]
          combine :: [(Position, Vec)] -> Vec
combine [] = [Char] -> Vec
forall a. HasCallStack => [Char] -> a
error [Char]
"compositeSimpson: odd number of half-intervals" -- this should never happen
          combine [(Position, Vec)
_] = Vec
forall v. AdditiveGroup v => v
zeroV
          combine ((Position, Vec)
_:(Position, Vec)
_:[]) = [Char] -> Vec
forall a. HasCallStack => [Char] -> a
error [Char]
"compositeSimpson: odd number of half-intervals" -- this should never happen
          combine ((Position
c0,Vec
f0):(Position
c1,Vec
f1):(Position
c2,Vec
f2):[(Position, Vec)]
ps)
              = Vec -> Vec -> Vec -> Vec -> Vec -> Vec
crossedSimp Vec
f0 Vec
f1 Vec
f2 (Position -> VectorField
displacement Position
c0 Position
c1) (Position -> VectorField
displacement Position
c1 Position
c2)
                Vec -> Vec -> Vec
forall v. AdditiveGroup v => v -> v -> v
^+^ [(Position, Vec)] -> Vec
combine ((Position
c2,Vec
f2)(Position, Vec) -> [(Position, Vec)] -> [(Position, Vec)]
forall a. a -> [a] -> [a]
:[(Position, Vec)]
ps)
      in [(Position, Vec)] -> Vec
combine [(Position, Vec)]
pairs