{-# OPTIONS -Wall #-}

{- | 
Module      :  LPFPCore.Integrals
Copyright   :  (c) Scott N. Walck 2023
License     :  BSD3 (see LICENSE)
Maintainer  :  Scott N. Walck <walck@lvc.edu>
Stability   :  stable

Code needed in chapter 24 of the book Learn Physics with Functional Programming
-}

module LPFPCore.Integrals where

import LPFPCore.SimpleVec
    ( R, Vec, (^+^), (*^), (^*), (^/), (<.>), (><), sumV, magnitude )
import LPFPCore.CoordinateSystems ( Position, ScalarField, VectorField
                         , displacement, shiftPosition, origin, rVF )
import LPFPCore.Geometry ( Curve(..), Surface(..), Volume(..) )

type CurveApprox = Curve -> [(Position,Vec)]

type SurfaceApprox = Surface -> [(Position,Vec)]

type VolumeApprox = Volume -> [(Position,R)]

scalarLineIntegral :: CurveApprox -> ScalarField -> Curve -> R
scalarLineIntegral :: CurveApprox -> ScalarField -> Curve -> R
scalarLineIntegral CurveApprox
approx ScalarField
f Curve
c
    = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ScalarField
f Position
r' forall a. Num a => a -> a -> a
* Vec -> R
magnitude Vec
dl' | (Position
r',Vec
dl') <- CurveApprox
approx Curve
c]

scalarSurfaceIntegral :: SurfaceApprox -> ScalarField -> Surface -> R
scalarSurfaceIntegral :: SurfaceApprox -> ScalarField -> Surface -> R
scalarSurfaceIntegral SurfaceApprox
approx ScalarField
f Surface
s
    = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ScalarField
f Position
r' forall a. Num a => a -> a -> a
* Vec -> R
magnitude Vec
da' | (Position
r',Vec
da') <- SurfaceApprox
approx Surface
s]

scalarVolumeIntegral :: VolumeApprox -> ScalarField -> Volume -> R
scalarVolumeIntegral :: VolumeApprox -> ScalarField -> Volume -> R
scalarVolumeIntegral VolumeApprox
approx ScalarField
f Volume
vol
    = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ScalarField
f Position
r' forall a. Num a => a -> a -> a
* R
dv' | (Position
r',R
dv') <- VolumeApprox
approx Volume
vol]

vectorLineIntegral :: CurveApprox -> VectorField -> Curve -> Vec
vectorLineIntegral :: CurveApprox -> VectorField -> Curve -> Vec
vectorLineIntegral CurveApprox
approx VectorField
vF Curve
c
    = [Vec] -> Vec
sumV [VectorField
vF Position
r' Vec -> R -> Vec
^* Vec -> R
magnitude Vec
dl' | (Position
r',Vec
dl') <- CurveApprox
approx Curve
c]

vectorSurfaceIntegral :: SurfaceApprox -> VectorField -> Surface -> Vec
vectorSurfaceIntegral :: SurfaceApprox -> VectorField -> Surface -> Vec
vectorSurfaceIntegral SurfaceApprox
approx VectorField
vF Surface
s
    = [Vec] -> Vec
sumV [VectorField
vF Position
r' Vec -> R -> Vec
^* Vec -> R
magnitude Vec
da' | (Position
r',Vec
da') <- SurfaceApprox
approx Surface
s]

vectorVolumeIntegral :: VolumeApprox -> VectorField -> Volume -> Vec
vectorVolumeIntegral :: VolumeApprox -> VectorField -> Volume -> Vec
vectorVolumeIntegral VolumeApprox
approx VectorField
vF Volume
vol
    = [Vec] -> Vec
sumV [VectorField
vF Position
r' Vec -> R -> Vec
^* R
dv' | (Position
r',R
dv') <- VolumeApprox
approx Volume
vol]

curveSample :: Int -> Curve -> [(Position,Vec)]
curveSample :: Int -> CurveApprox
curveSample Int
n Curve
c
    = let segCent :: Segment -> Position
          segCent :: Segment -> Position
segCent (Position
p1,Position
p2) = Vec -> Position -> Position
shiftPosition ((VectorField
rVF Position
p1 Vec -> Vec -> Vec
^+^ VectorField
rVF Position
p2) Vec -> R -> Vec
^/ R
2) Position
origin
          segDisp :: Segment -> Vec
          segDisp :: Segment -> Vec
segDisp = forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Position -> VectorField
displacement
      in [(Segment -> Position
segCent Segment
seg, Segment -> Vec
segDisp Segment
seg) | Segment
seg <- Int -> Curve -> [Segment]
segments Int
n Curve
c]

type Segment = (Position,Position)
segments :: Int -> Curve -> [Segment]
segments :: Int -> Curve -> [Segment]
segments Int
n (Curve R -> Position
g R
a R
b)
    = let ps :: [Position]
ps = forall a b. (a -> b) -> [a] -> [b]
map R -> Position
g forall a b. (a -> b) -> a -> b
$ Int -> R -> R -> [R]
linSpaced Int
n R
a R
b
      in forall a b. [a] -> [b] -> [(a, b)]
zip [Position]
ps (forall a. [a] -> [a]
tail [Position]
ps)

linSpaced :: Int -> R -> R -> [R]
linSpaced :: Int -> R -> R -> [R]
linSpaced Int
n R
x0 R
x1 = forall a. Int -> [a] -> [a]
take (Int
nforall a. Num a => a -> a -> a
+Int
1) [R
x0, R
x0forall a. Num a => a -> a -> a
+R
dx .. R
x1]
    where dx :: R
dx = (R
x1 forall a. Num a => a -> a -> a
- R
x0) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n

surfaceSample :: Int -> Surface -> [(Position,Vec)]
surfaceSample :: Int -> SurfaceApprox
surfaceSample Int
n Surface
s = [(Triangle -> Position
triCenter Triangle
tri, Triangle -> Vec
triArea Triangle
tri) | Triangle
tri <- Int -> Surface -> [Triangle]
triangles Int
n Surface
s]

data Triangle = Tri Position Position Position

triCenter :: Triangle -> Position
triCenter :: Triangle -> Position
triCenter (Tri Position
p1 Position
p2 Position
p3)
    = Vec -> Position -> Position
shiftPosition ((VectorField
rVF Position
p1 Vec -> Vec -> Vec
^+^ VectorField
rVF Position
p2 Vec -> Vec -> Vec
^+^ VectorField
rVF Position
p3) Vec -> R -> Vec
^/ R
3) Position
origin

triArea :: Triangle -> Vec  -- vector area
triArea :: Triangle -> Vec
triArea (Tri Position
p1 Position
p2 Position
p3) = R
0.5 R -> Vec -> Vec
*^ (Position -> VectorField
displacement Position
p1 Position
p2 Vec -> Vec -> Vec
>< Position -> VectorField
displacement Position
p2 Position
p3)

triangles :: Int -> Surface -> [Triangle]
triangles :: Int -> Surface -> [Triangle]
triangles Int
n (Surface (R, R) -> Position
g R
sl R
su R -> R
tl R -> R
tu)
    = let sts :: [[(R, R)]]
sts = [[(R
s,R
t) | R
t <- Int -> R -> R -> [R]
linSpaced Int
n (R -> R
tl R
s) (R -> R
tu R
s)]
                     | R
s <- Int -> R -> R -> [R]
linSpaced Int
n R
sl R
su]
          stSquares :: [((R, R), (R, R), (R, R), (R, R))]
stSquares = [( [[(R, R)]]
sts forall a. [a] -> Int -> a
!! Int
j     forall a. [a] -> Int -> a
!! Int
k
                       , [[(R, R)]]
sts forall a. [a] -> Int -> a
!! (Int
jforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!! Int
k
                       , [[(R, R)]]
sts forall a. [a] -> Int -> a
!! (Int
jforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!! (Int
kforall a. Num a => a -> a -> a
+Int
1)
                       , [[(R, R)]]
sts forall a. [a] -> Int -> a
!! Int
j     forall a. [a] -> Int -> a
!! (Int
kforall a. Num a => a -> a -> a
+Int
1))
                      | Int
j <- [Int
0..Int
nforall a. Num a => a -> a -> a
-Int
1], Int
k <- [Int
0..Int
nforall a. Num a => a -> a -> a
-Int
1]]
          twoTriangles :: ((R, R), (R, R), (R, R), (R, R)) -> [Triangle]
twoTriangles ((R, R)
pp1,(R, R)
pp2,(R, R)
pp3,(R, R)
pp4)
              = [Position -> Position -> Position -> Triangle
Tri ((R, R) -> Position
g (R, R)
pp1) ((R, R) -> Position
g (R, R)
pp2) ((R, R) -> Position
g (R, R)
pp3),Position -> Position -> Position -> Triangle
Tri ((R, R) -> Position
g (R, R)
pp1) ((R, R) -> Position
g (R, R)
pp3) ((R, R) -> Position
g (R, R)
pp4)]
      in forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap ((R, R), (R, R), (R, R), (R, R)) -> [Triangle]
twoTriangles [((R, R), (R, R), (R, R), (R, R))]
stSquares

volumeSample :: Int -> Volume -> [(Position,R)]
volumeSample :: Int -> VolumeApprox
volumeSample Int
n Volume
v = [(Tet -> Position
tetCenter Tet
tet, Tet -> R
tetVolume Tet
tet) | Tet
tet <- Int -> Volume -> [Tet]
tetrahedrons Int
n Volume
v]

data Tet = Tet Position Position Position Position

tetCenter :: Tet -> Position
tetCenter :: Tet -> Position
tetCenter (Tet Position
p1 Position
p2 Position
p3 Position
p4)
    = Vec -> Position -> Position
shiftPosition ((VectorField
rVF Position
p1 Vec -> Vec -> Vec
^+^ VectorField
rVF Position
p2 Vec -> Vec -> Vec
^+^ VectorField
rVF Position
p3 Vec -> Vec -> Vec
^+^ VectorField
rVF Position
p4) Vec -> R -> Vec
^/ R
4) Position
origin

tetVolume :: Tet -> R
tetVolume :: Tet -> R
tetVolume (Tet Position
p1 Position
p2 Position
p3 Position
p4)
    = forall a. Num a => a -> a
abs forall a b. (a -> b) -> a -> b
$ (Vec
d1 Vec -> Vec -> R
<.> (Vec
d2 Vec -> Vec -> Vec
>< Vec
d3)) forall a. Fractional a => a -> a -> a
/ R
6
      where
        d1 :: Vec
d1 = Position -> VectorField
displacement Position
p1 Position
p4
        d2 :: Vec
d2 = Position -> VectorField
displacement Position
p2 Position
p4
        d3 :: Vec
d3 = Position -> VectorField
displacement Position
p3 Position
p4

data ParamCube
    = PC { ParamCube -> (R, R, R)
v000 :: (R,R,R)
         , ParamCube -> (R, R, R)
v001 :: (R,R,R)
         , ParamCube -> (R, R, R)
v010 :: (R,R,R)
         , ParamCube -> (R, R, R)
v011 :: (R,R,R)
         , ParamCube -> (R, R, R)
v100 :: (R,R,R)
         , ParamCube -> (R, R, R)
v101 :: (R,R,R)
         , ParamCube -> (R, R, R)
v110 :: (R,R,R)
         , ParamCube -> (R, R, R)
v111 :: (R,R,R)
         }

tetrahedrons :: Int -> Volume -> [Tet]
tetrahedrons :: Int -> Volume -> [Tet]
tetrahedrons Int
n (Volume (R, R, R) -> Position
g R
sl R
su R -> R
tl R -> R
tu R -> R -> R
ul R -> R -> R
uu)
    = let stus :: [[[(R, R, R)]]]
stus = [[[(R
s,R
t,R
u) | R
u <- Int -> R -> R -> [R]
linSpaced Int
n (R -> R -> R
ul R
s R
t) (R -> R -> R
uu R
s R
t)]
                            | R
t <- Int -> R -> R -> [R]
linSpaced Int
n (R -> R
tl R
s) (R -> R
tu R
s)]
                            | R
s <- Int -> R -> R -> [R]
linSpaced Int
n R
sl R
su]
          stCubes :: [ParamCube]
stCubes = [(R, R, R)
-> (R, R, R)
-> (R, R, R)
-> (R, R, R)
-> (R, R, R)
-> (R, R, R)
-> (R, R, R)
-> (R, R, R)
-> ParamCube
PC ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!!  Int
j    forall a. [a] -> Int -> a
!!  Int
k    forall a. [a] -> Int -> a
!!  Int
l   )
                        ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!!  Int
j    forall a. [a] -> Int -> a
!!  Int
k    forall a. [a] -> Int -> a
!! (Int
lforall a. Num a => a -> a -> a
+Int
1))
                        ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!!  Int
j    forall a. [a] -> Int -> a
!! (Int
kforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!!  Int
l   )
                        ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!!  Int
j    forall a. [a] -> Int -> a
!! (Int
kforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!! (Int
lforall a. Num a => a -> a -> a
+Int
1))
                        ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!! (Int
jforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!!  Int
k    forall a. [a] -> Int -> a
!!  Int
l   )
                        ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!! (Int
jforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!!  Int
k    forall a. [a] -> Int -> a
!! (Int
lforall a. Num a => a -> a -> a
+Int
1))
                        ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!! (Int
jforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!! (Int
kforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!!  Int
l   )
                        ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!! (Int
jforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!! (Int
kforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!! (Int
lforall a. Num a => a -> a -> a
+Int
1))
                    | Int
j <- [Int
0..Int
nforall a. Num a => a -> a -> a
-Int
1], Int
k <- [Int
0..Int
nforall a. Num a => a -> a -> a
-Int
1], Int
l <- [Int
0..Int
nforall a. Num a => a -> a -> a
-Int
1]]
          tets :: ParamCube -> [Tet]
tets (PC (R, R, R)
c000 (R, R, R)
c001 (R, R, R)
c010 (R, R, R)
c011 (R, R, R)
c100 (R, R, R)
c101 (R, R, R)
c110 (R, R, R)
c111)
              = [Position -> Position -> Position -> Position -> Tet
Tet ((R, R, R) -> Position
g (R, R, R)
c000) ((R, R, R) -> Position
g (R, R, R)
c100) ((R, R, R) -> Position
g (R, R, R)
c010) ((R, R, R) -> Position
g (R, R, R)
c001)
                ,Position -> Position -> Position -> Position -> Tet
Tet ((R, R, R) -> Position
g (R, R, R)
c011) ((R, R, R) -> Position
g (R, R, R)
c111) ((R, R, R) -> Position
g (R, R, R)
c001) ((R, R, R) -> Position
g (R, R, R)
c010)
                ,Position -> Position -> Position -> Position -> Tet
Tet ((R, R, R) -> Position
g (R, R, R)
c110) ((R, R, R) -> Position
g (R, R, R)
c010) ((R, R, R) -> Position
g (R, R, R)
c100) ((R, R, R) -> Position
g (R, R, R)
c111)
                ,Position -> Position -> Position -> Position -> Tet
Tet ((R, R, R) -> Position
g (R, R, R)
c101) ((R, R, R) -> Position
g (R, R, R)
c001) ((R, R, R) -> Position
g (R, R, R)
c111) ((R, R, R) -> Position
g (R, R, R)
c100)
                ,Position -> Position -> Position -> Position -> Tet
Tet ((R, R, R) -> Position
g (R, R, R)
c111) ((R, R, R) -> Position
g (R, R, R)
c100) ((R, R, R) -> Position
g (R, R, R)
c010) ((R, R, R) -> Position
g (R, R, R)
c001)
                ]
      in forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap ParamCube -> [Tet]
tets [ParamCube]
stCubes