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

{- | 
Module      :  Physics.Learn.Surface
Copyright   :  (c) Scott N. Walck 2012-2019
License     :  BSD3 (see LICENSE)
Maintainer  :  Scott N. Walck <walck@lvc.edu>
Stability   :  experimental

This module contains functions for working with 'Surface's
and surface integrals over 'Surface's.
-}

module Physics.Learn.Surface
    (
    -- * Surfaces
      Surface(..)
    , unitSphere
    , centeredSphere
    , sphere
    , northernHemisphere
    , disk
    , shiftSurface
    -- * Surface Integrals
    , surfaceIntegral
    , dottedSurfaceIntegral
    )
    where

import Data.VectorSpace
    ( VectorSpace
    , InnerSpace
    , Scalar
    )
import Physics.Learn.CarrotVec
    ( vec
    , (^+^)
    , (^-^)
    , (^*)
    , (^/)
    , (<.>)
    , (><)
    , magnitude
    , sumV
    )
import Physics.Learn.Position
    ( Position
    , Displacement
    , VectorField
    , Field
    , cart
    , cyl
    , shiftPosition
    , displacement
    )

-- | Surface is a parametrized function from two parameters to space,
--   lower and upper limits on the first parameter, and
--   lower and upper limits for the second parameter
--   (expressed as functions of the first parameter).
data Surface = Surface { Surface -> (Double, Double) -> Position
surfaceFunc :: (Double,Double) -> Position  -- ^ function from two parameters (s,t) into space
                       , Surface -> Double
lowerLimit :: Double            -- ^ s_l
                       , Surface -> Double
upperLimit :: Double            -- ^ s_u
                       , Surface -> Double -> Double
lowerCurve :: Double -> Double  -- ^ t_l(s)
                       , Surface -> Double -> Double
upperCurve :: Double -> Double  -- ^ t_u(s)
                       }

-- | A unit sphere, centered at the origin.
unitSphere :: Surface
unitSphere :: Surface
unitSphere = ((Double, Double) -> Position)
-> Double
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> Surface
Surface (\(Double
th,Double
phi) -> Double -> Double -> Double -> Position
cart (Double -> Double
forall a. Floating a => a -> a
sin Double
th Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
phi) (Double -> Double
forall a. Floating a => a -> a
sin Double
th Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
phi) (Double -> Double
forall a. Floating a => a -> a
cos Double
th))
             Double
0 Double
forall a. Floating a => a
pi (Double -> Double -> Double
forall a b. a -> b -> a
const Double
0) (Double -> Double -> Double
forall a b. a -> b -> a
const (Double -> Double -> Double) -> Double -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi)

-- | A sphere with given radius centered at the origin.
centeredSphere :: Double -> Surface
centeredSphere :: Double -> Surface
centeredSphere Double
r = ((Double, Double) -> Position)
-> Double
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> Surface
Surface (\(Double
th,Double
phi) -> Double -> Double -> Double -> Position
cart (Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
th Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
phi) (Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
th Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
phi) (Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
th))
                   Double
0 Double
forall a. Floating a => a
pi (Double -> Double -> Double
forall a b. a -> b -> a
const Double
0) (Double -> Double -> Double
forall a b. a -> b -> a
const (Double -> Double -> Double) -> Double -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi)

-- | Sphere with given radius and center.
sphere :: Double -> Position -> Surface
sphere :: Double -> Position -> Surface
sphere Double
r Position
c = ((Double, Double) -> Position)
-> Double
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> Surface
Surface (\(Double
th,Double
phi) -> Vec -> Position -> Position
shiftPosition (Double -> Double -> Double -> Vec
vec (Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
th Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
phi) (Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
th Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
phi) (Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
th)) Position
c)
             Double
0 Double
forall a. Floating a => a
pi (Double -> Double -> Double
forall a b. a -> b -> a
const Double
0) (Double -> Double -> Double
forall a b. a -> b -> a
const (Double -> Double -> Double) -> Double -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi)

-- | The upper half of a unit sphere, centered at the origin.
northernHemisphere :: Surface
northernHemisphere :: Surface
northernHemisphere = ((Double, Double) -> Position)
-> Double
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> Surface
Surface (\(Double
th,Double
phi) -> Double -> Double -> Double -> Position
cart (Double -> Double
forall a. Floating a => a -> a
sin Double
th Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
phi) (Double -> Double
forall a. Floating a => a -> a
sin Double
th Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
phi) (Double -> Double
forall a. Floating a => a -> a
cos Double
th))
                     Double
0 (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) (Double -> Double -> Double
forall a b. a -> b -> a
const Double
0) (Double -> Double -> Double
forall a b. a -> b -> a
const (Double -> Double -> Double) -> Double -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi)

-- | A disk with given radius, centered at the origin.
disk :: Double -> Surface
disk :: Double -> Surface
disk Double
radius = ((Double, Double) -> Position)
-> Double
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> Surface
Surface (\(Double
s,Double
phi) -> Double -> Double -> Double -> Position
cyl Double
s Double
phi Double
0) Double
0 Double
radius (Double -> Double -> Double
forall a b. a -> b -> a
const Double
0) (Double -> Double -> Double
forall a b. a -> b -> a
const (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi))

-- To do : boundaryOfSurface :: Surface -> Curve

-- | A plane surface integral, in which area element is a scalar.
surfaceIntegral :: (VectorSpace v, Scalar v ~ Double) =>
                   Int      -- ^ number of intervals for first parameter, s
                -> Int      -- ^ number of intervals for second parameter, t
                -> Field v  -- ^ the scalar or vector field to integrate
                -> Surface  -- ^ the surface over which to integrate
                -> v        -- ^ the resulting scalar or vector
surfaceIntegral :: forall v.
(VectorSpace v, Scalar v ~ Double) =>
Int -> Int -> Field v -> Surface -> v
surfaceIntegral Int
n1 Int
n2 Field v
field (Surface (Double, Double) -> Position
f Double
s_l Double
s_u Double -> Double
t_l Double -> Double
t_u)
    = [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] -> v) -> [[v]] -> [v]
forall a b. (a -> b) -> [a] -> [b]
map [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] -> [Double] -> [v]) -> [[v]] -> [[Double]] -> [[v]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((v -> Double -> v) -> [v] -> [Double] -> [v]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith v -> Double -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
(^*)) [[v]]
aveVals (([Vec] -> [Double]) -> [[Vec]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ((Vec -> Double) -> [Vec] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Vec -> Double
forall v s. (InnerSpace v, s ~ Scalar v, Floating s) => v -> s
magnitude) [[Vec]]
areas)
      where
        pts :: [[Position]]
pts = [[(Double, Double) -> Position
f (Double
s,Double
t) | Double
t <- Int -> Double -> Double -> [Double]
linSpaced Int
n2 (Double -> Double
t_l Double
s) (Double -> Double
t_u Double
s)] | Double
s <- Int -> Double -> Double -> [Double]
linSpaced Int
n1 Double
s_l Double
s_u]
        areas :: [[Vec]]
areas = ([Vec] -> [Vec] -> [Vec]) -> [[Vec]] -> [[Vec]] -> [[Vec]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((Vec -> Vec -> Vec) -> [Vec] -> [Vec] -> [Vec]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Vec -> Vec -> Vec
(><)) [[Vec]]
dus [[Vec]]
dvs
        dus :: [[Vec]]
dus = ([Position] -> [Position] -> [Vec])
-> [[Position]] -> [[Position]] -> [[Vec]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((Position -> Position -> Vec) -> [Position] -> [Position] -> [Vec]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Position -> Position -> Vec
displacement) [[Position]]
pts ([[Position]] -> [[Position]]
forall a. HasCallStack => [a] -> [a]
tail [[Position]]
pts)
        dvs :: [[Vec]]
dvs = ([Position] -> [Vec]) -> [[Position]] -> [[Vec]]
forall a b. (a -> b) -> [a] -> [b]
map (\[Position]
row -> (Position -> Position -> Vec) -> [Position] -> [Position] -> [Vec]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Position -> Position -> Vec
displacement [Position]
row ([Position] -> [Position]
forall a. HasCallStack => [a] -> [a]
tail [Position]
row)) [[Position]]
pts
        vals :: [[v]]
vals = ([Position] -> [v]) -> [[Position]] -> [[v]]
forall a b. (a -> b) -> [a] -> [b]
map (Field v -> [Position] -> [v]
forall a b. (a -> b) -> [a] -> [b]
map Field v
field) [[Position]]
pts
        halfAveVals :: [[v]]
halfAveVals = ([v] -> [v]) -> [[v]] -> [[v]]
forall a b. (a -> b) -> [a] -> [b]
map (\[v]
row -> (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 ~ Double) => v -> v -> v
ave ([v] -> [v]
forall a. HasCallStack => [a] -> [a]
tail [v]
row) [v]
row) [[v]]
vals
        aveVals :: [[v]]
aveVals = ([v] -> [v] -> [v]) -> [[v]] -> [[v]] -> [[v]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((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 ~ Double) => v -> v -> v
ave) ([[v]] -> [[v]]
forall a. HasCallStack => [a] -> [a]
tail [[v]]
halfAveVals) [[v]]
halfAveVals

-- | A dotted surface integral, in which area element is a vector.
dottedSurfaceIntegral :: Int          -- ^ number of intervals for first parameter, s
                      -> Int          -- ^ number of intervals for second parameter, t
                      -> VectorField  -- ^ the vector field to integrate
                      -> Surface      -- ^ the surface over which to integrate
                      -> Double       -- ^ the resulting scalar
dottedSurfaceIntegral :: Int -> Int -> (Position -> Vec) -> Surface -> Double
dottedSurfaceIntegral Int
n1 Int
n2 Position -> Vec
vf (Surface (Double, Double) -> Position
f Double
s_l Double
s_u Double -> Double
t_l Double -> Double
t_u)
    = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Double] -> Double) -> [Double] -> Double
forall a b. (a -> b) -> a -> b
$ ([Double] -> Double) -> [[Double]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([[Double]] -> [Double]) -> [[Double]] -> [Double]
forall a b. (a -> b) -> a -> b
$ ([Vec] -> [Vec] -> [Double]) -> [[Vec]] -> [[Vec]] -> [[Double]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((Vec -> Vec -> Double) -> [Vec] -> [Vec] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Vec -> Vec -> Double
Vec -> Vec -> Scalar Vec
forall v. InnerSpace v => v -> v -> Scalar v
(<.>)) [[Vec]]
aveVals [[Vec]]
areas
      where
        pts :: [[Position]]
pts = [[(Double, Double) -> Position
f (Double
s,Double
t) | Double
t <- Int -> Double -> Double -> [Double]
linSpaced Int
n2 (Double -> Double
t_l Double
s) (Double -> Double
t_u Double
s)] | Double
s <- Int -> Double -> Double -> [Double]
linSpaced Int
n1 Double
s_l Double
s_u]
        areas :: [[Vec]]
areas = ([Vec] -> [Vec] -> [Vec]) -> [[Vec]] -> [[Vec]] -> [[Vec]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((Vec -> Vec -> Vec) -> [Vec] -> [Vec] -> [Vec]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Vec -> Vec -> Vec
(><)) [[Vec]]
dus [[Vec]]
dvs
        dus :: [[Vec]]
dus = ([Position] -> [Position] -> [Vec])
-> [[Position]] -> [[Position]] -> [[Vec]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((Position -> Position -> Vec) -> [Position] -> [Position] -> [Vec]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Position -> Position -> Vec
displacement) [[Position]]
pts ([[Position]] -> [[Position]]
forall a. HasCallStack => [a] -> [a]
tail [[Position]]
pts)
        dvs :: [[Vec]]
dvs = ([Position] -> [Vec]) -> [[Position]] -> [[Vec]]
forall a b. (a -> b) -> [a] -> [b]
map (\[Position]
row -> (Position -> Position -> Vec) -> [Position] -> [Position] -> [Vec]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Position -> Position -> Vec
displacement [Position]
row ([Position] -> [Position]
forall a. HasCallStack => [a] -> [a]
tail [Position]
row)) [[Position]]
pts
        vals :: [[Vec]]
vals = ([Position] -> [Vec]) -> [[Position]] -> [[Vec]]
forall a b. (a -> b) -> [a] -> [b]
map ((Position -> Vec) -> [Position] -> [Vec]
forall a b. (a -> b) -> [a] -> [b]
map Position -> Vec
vf) [[Position]]
pts
        halfAveVals :: [[Vec]]
halfAveVals = ([Vec] -> [Vec]) -> [[Vec]] -> [[Vec]]
forall a b. (a -> b) -> [a] -> [b]
map (\[Vec]
row -> (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 ~ Double) => v -> v -> v
ave ([Vec] -> [Vec]
forall a. HasCallStack => [a] -> [a]
tail [Vec]
row) [Vec]
row) [[Vec]]
vals
        aveVals :: [[Vec]]
aveVals = ([Vec] -> [Vec] -> [Vec]) -> [[Vec]] -> [[Vec]] -> [[Vec]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((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 ~ Double) => v -> v -> v
ave) ([[Vec]] -> [[Vec]]
forall a. HasCallStack => [a] -> [a]
tail [[Vec]]
halfAveVals) [[Vec]]
halfAveVals

{-
evalSquare :: (InnerSpace v, Scalar v ~ Double) => Double -> Int -> Int
             -> (Vec -> v) -> Surface
             -> Vec -> Vec -> Vec -> Vec
             -> v -> v -> v -> v -> v
evalSquare tol level maxlevel field (Surface f s_l s_u t_l t_u)
           surfll surflu surful surfuu fieldll fieldlu fieldul fielduu val
    = let s_m = (s_l + s_u) / 2
          t_m s = (t_l s + t_u s) / 2
          surflm = f (s_l,t_m s_l)
          surfum = f (s_u,t_m s_u)
          surfml = f (s_m,t_l s_m)
          surfmu = f (s_m,t_u s_m)
          surfmm = f (s_m,t_m s_m)
          fieldlm = field surflm
          fieldum = field surfum
          fieldml = field surfml
          fieldmu = field surfmu
          fieldmm = field surfmm
          dull = surfml ^-^ surfll
          dulu = surfmm ^-^ surflm
          duul = surful ^-^ surfml
          duuu = surfum ^-^ surfmm
          dvll = surflm ^-^ surfll
          dvlu = surflu ^-^ surflm
          dvul = surfmm ^-^ surfml
          dvuu = surfmu ^-^ surfmm
          areall = dull >< dvll
          arealu = dulu >< dvlu
          areaul = duul >< dvul
          areauu = duuu >< dvuu
          valll = average [fieldll,fieldlm,fieldml,fieldmm] <.> areall
          vallu = average [fieldlm,fieldlu,fieldmm,fieldmu] <.> arealu
          valul = average [fieldml,fieldmm,fieldul,fieldum] <.> areaul
          valuu = average [fieldmm,fieldmu,fieldum,fielduu] <.> areauu
          newval = valll ^+^ vallu ^+^ valul ^+^ valuu
      in if magnitude (newval ^-^ val) < tol then
             newval
         else
             evalSquare (tol/2) (level+1) maxlevel field (Surface f s_l s_m t_l t_m)
                        surfll surflm surfml surfmm fieldll fieldlm fieldml fieldmm valll ^+^
             evalSquare (tol/2) (level+1) maxlevel field (Surface f s_l s_m t_m t_u)
                        surflm surflu surfmm surfmu fieldlm fieldlu fieldmm fieldmu vallu ^+^
             evalSquare (tol/2) (level+1) maxlevel field (Surface f s_m s_u t_l t_m)
                        surfml surfmm surful surfum fieldml fieldmm fieldul fieldum valul ^+^
             evalSquare (tol/2) (level+1) maxlevel field (Surface f s_m s_u t_m t_u)
                        surfmm surfmu surfum surfuu fieldmm fieldmu fieldum fielduu valuu
-}

{-
dottedSurfIntegral :: Double
                   -> (Vec -> Vec) -> Surface
                   -> Double
dottedSurfIntegral tol vf (Surface f s_l s_u t_l t_u)
    = let surfll = f (s_l,t_l s_l)
          surflu = f (s_l,t_u s_l)
          surful = f (s_u,t_l s_u)
          surfuu = f (s_u,t_u s_u)
          fieldll = vf surfll
          fieldlu = vf surflu
          fieldul = vf surful
          fielduu = vf surfuu
          du = surful ^-^ surfll
          dv = surflu ^-^ surfll
          area = du >< dv
          val = average [fieldll,fieldlu,fieldul,fielduu] <.> area
      in evalSquare tol 1 2 20 vf (Surface f s_l s_u t_l t_u)
         surfll surflu surful surfuu fieldll fieldlu fieldul fielduu val

fullDottedSurfIntegral :: Double -> Int -> Int
                       -> (Vec -> Vec) -> Surface
                       -> Double
fullDottedSurfIntegral tol minlevel maxlevel vf (Surface f s_l s_u t_l t_u)
    = let surfll = f (s_l,t_l s_l)
          surflu = f (s_l,t_u s_l)
          surful = f (s_u,t_l s_u)
          surfuu = f (s_u,t_u s_u)
          fieldll = vf surfll
          fieldlu = vf surflu
          fieldul = vf surful
          fielduu = vf surfuu
          du = surful ^-^ surfll
          dv = surflu ^-^ surfll
          area = du >< dv
          val = average [fieldll,fieldlu,fieldul,fielduu] <.> area
      in evalSquare tol 1 minlevel maxlevel vf (Surface f s_l s_u t_l t_u)
         surfll surflu surful surfuu fieldll fieldlu fieldul fielduu val

evalSquare :: Double -> Int -> Int -> Int
           -> (Vec -> Vec) -> Surface
           -> Vec -> Vec -> Vec -> Vec
           -> Vec -> Vec -> Vec -> Vec -> Double -> Double
evalSquare tol level minlevel maxlevel field (Surface f s_l s_u t_l t_u)
           surfll surflu surful surfuu fieldll fieldlu fieldul fielduu val
    = let s_m = (s_l + s_u) / 2
          t_m s = (t_l s + t_u s) / 2
          surflm = f (s_l,t_m s_l)
          surfum = f (s_u,t_m s_u)
          surfml = f (s_m,t_l s_m)
          surfmu = f (s_m,t_u s_m)
          surfmm = f (s_m,t_m s_m)
          fieldlm = field surflm
          fieldum = field surfum
          fieldml = field surfml
          fieldmu = field surfmu
          fieldmm = field surfmm
          dull = surfml ^-^ surfll
          dulu = surfmm ^-^ surflm
          duul = surful ^-^ surfml
          duuu = surfum ^-^ surfmm
          dvll = surflm ^-^ surfll
          dvlu = surflu ^-^ surflm
          dvul = surfmm ^-^ surfml
          dvuu = surfmu ^-^ surfmm
          areall = dull >< dvll
          arealu = dulu >< dvlu
          areaul = duul >< dvul
          areauu = duuu >< dvuu
          valll = average [fieldll,fieldlm,fieldml,fieldmm] <.> areall
          vallu = average [fieldlm,fieldlu,fieldmm,fieldmu] <.> arealu
          valul = average [fieldml,fieldmm,fieldul,fieldum] <.> areaul
          valuu = average [fieldmm,fieldmu,fieldum,fielduu] <.> areauu
          newval = valll + vallu + valul + valuu
      in if level >= maxlevel || level >= minlevel && abs (newval - val) < tol then
             newval
         else
             evalSquare (tol/4) (level+1) minlevel maxlevel field (Surface f s_l s_m t_l t_m)
                        surfll surflm surfml surfmm fieldll fieldlm fieldml fieldmm valll +
             evalSquare (tol/4) (level+1) minlevel maxlevel field (Surface f s_l s_m t_m t_u)
                        surflm surflu surfmm surfmu fieldlm fieldlu fieldmm fieldmu vallu +
             evalSquare (tol/4) (level+1) minlevel maxlevel field (Surface f s_m s_u t_l t_m)
                        surfml surfmm surful surfum fieldml fieldmm fieldul fieldum valul +
             evalSquare (tol/4) (level+1) minlevel maxlevel field (Surface f s_m s_u t_m t_u)
                        surfmm surfmu surfum surfuu fieldmm fieldmu fieldum fielduu valuu
-}

-- n+1 points
linSpaced :: Int -> Double -> Double -> [Double]
linSpaced :: Int -> Double -> Double -> [Double]
linSpaced Int
n Double
a Double
b
    | Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
b      = let dx :: Double
dx = (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
                   in [Double
a,Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
dx..Double
b]
    | Double
a Double -> Double -> Bool
forall v. (InnerSpace v, Scalar v ~ Double) => v -> v -> Bool
~~ Double
b     = [Double -> Double -> Double
forall v. (VectorSpace v, Scalar v ~ Double) => v -> v -> v
ave Double
a Double
b]
    | Bool
otherwise  = [Char] -> [Double]
forall a. HasCallStack => [Char] -> a
error ([Char] -> [Double]) -> [Char] -> [Double]
forall a b. (a -> b) -> a -> b
$ [Char]
"linSpaced:  lower limit greater than upper limit:  (n,a,b) = " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ (Int, Double, Double) -> [Char]
forall a. Show a => a -> [Char]
show (Int
n,Double
a,Double
b)

(~~) :: (InnerSpace v, Scalar v ~ Double) => v -> v -> Bool
v
a ~~ :: forall v. (InnerSpace v, Scalar v ~ Double) => v -> v -> Bool
~~ v
b = v -> Double
forall v s. (InnerSpace v, s ~ Scalar v, Floating s) => v -> s
magnitude (v
a v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^-^ v
b) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
tolerance

tolerance :: Double
tolerance :: Double
tolerance = Double
1e-10

ave :: (VectorSpace v, Scalar v ~ Double) => v -> v -> v
ave :: forall v. (VectorSpace v, Scalar v ~ Double) => v -> v -> v
ave v
v1 v
v2 = (v
v1 v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^+^ v
v2) v -> Double -> v
forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/ Double
2

{-
average :: (VectorSpace v, Scalar v ~ Double) => [v] -> v
average vs = sumV vs ^/ fromIntegral (length vs)

areaOfSurface :: Surface -> Double
areaOfSurface = surfaceIntegral 100 100 (const 1)
-}

-- | Shift a surface by a displacement.
shiftSurface :: Displacement -> Surface -> Surface
shiftSurface :: Vec -> Surface -> Surface
shiftSurface Vec
d (Surface (Double, Double) -> Position
f Double
sl Double
su Double -> Double
tl Double -> Double
tu)
    = ((Double, Double) -> Position)
-> Double
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> Surface
Surface (Vec -> Position -> Position
shiftPosition Vec
d (Position -> Position)
-> ((Double, Double) -> Position) -> (Double, Double) -> Position
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double, Double) -> Position
f) Double
sl Double
su Double -> Double
tl Double -> Double
tu