{-# 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 (forall a. Floating a => a -> a
sin Double
th forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos Double
phi) (forall a. Floating a => a -> a
sin Double
th forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin Double
phi) (forall a. Floating a => a -> a
cos Double
th))
             Double
0 forall a. Floating a => a
pi (forall a b. a -> b -> a
const Double
0) (forall a b. a -> b -> a
const forall a b. (a -> b) -> a -> b
$ Double
2forall a. Num a => a -> a -> a
*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 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin Double
th forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos Double
phi) (Double
r forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin Double
th forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin Double
phi) (Double
r forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos Double
th))
                   Double
0 forall a. Floating a => a
pi (forall a b. a -> b -> a
const Double
0) (forall a b. a -> b -> a
const forall a b. (a -> b) -> a -> b
$ Double
2forall a. Num a => a -> a -> a
*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 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin Double
th forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos Double
phi) (Double
r forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin Double
th forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin Double
phi) (Double
r forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos Double
th)) Position
c)
             Double
0 forall a. Floating a => a
pi (forall a b. a -> b -> a
const Double
0) (forall a b. a -> b -> a
const forall a b. (a -> b) -> a -> b
$ Double
2forall a. Num a => a -> a -> a
*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 (forall a. Floating a => a -> a
sin Double
th forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos Double
phi) (forall a. Floating a => a -> a
sin Double
th forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin Double
phi) (forall a. Floating a => a -> a
cos Double
th))
                     Double
0 (forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Double
2) (forall a b. a -> b -> a
const Double
0) (forall a b. a -> b -> a
const forall a b. (a -> b) -> a -> b
$ Double
2forall a. Num a => a -> a -> a
*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 (forall a b. a -> b -> a
const Double
0) (forall a b. a -> b -> a
const (Double
2forall a. Num a => a -> a -> a
*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)
    = forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
(^*)) [[v]]
aveVals (forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map 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 = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Vec -> Vec -> Vec
(><)) [[Vec]]
dus [[Vec]]
dvs
        dus :: [[Vec]]
dus = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Position -> Position -> Vec
displacement) [[Position]]
pts (forall a. [a] -> [a]
tail [[Position]]
pts)
        dvs :: [[Vec]]
dvs = forall a b. (a -> b) -> [a] -> [b]
map (\[Position]
row -> forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Position -> Position -> Vec
displacement [Position]
row (forall a. [a] -> [a]
tail [Position]
row)) [[Position]]
pts
        vals :: [[v]]
vals = forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map Field v
field) [[Position]]
pts
        halfAveVals :: [[v]]
halfAveVals = forall a b. (a -> b) -> [a] -> [b]
map (\[v]
row -> forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall v. (VectorSpace v, Scalar v ~ Double) => v -> v -> v
ave (forall a. [a] -> [a]
tail [v]
row) [v]
row) [[v]]
vals
        aveVals :: [[v]]
aveVals = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall v. (VectorSpace v, Scalar v ~ Double) => v -> v -> v
ave) (forall a. [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)
    = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith 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 = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Vec -> Vec -> Vec
(><)) [[Vec]]
dus [[Vec]]
dvs
        dus :: [[Vec]]
dus = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Position -> Position -> Vec
displacement) [[Position]]
pts (forall a. [a] -> [a]
tail [[Position]]
pts)
        dvs :: [[Vec]]
dvs = forall a b. (a -> b) -> [a] -> [b]
map (\[Position]
row -> forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Position -> Position -> Vec
displacement [Position]
row (forall a. [a] -> [a]
tail [Position]
row)) [[Position]]
pts
        vals :: [[Vec]]
vals = forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map Position -> Vec
vf) [[Position]]
pts
        halfAveVals :: [[Vec]]
halfAveVals = forall a b. (a -> b) -> [a] -> [b]
map (\[Vec]
row -> forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall v. (VectorSpace v, Scalar v ~ Double) => v -> v -> v
ave (forall a. [a] -> [a]
tail [Vec]
row) [Vec]
row) [[Vec]]
vals
        aveVals :: [[Vec]]
aveVals = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall v. (VectorSpace v, Scalar v ~ Double) => v -> v -> v
ave) (forall a. [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 forall a. Ord a => a -> a -> Bool
< Double
b      = let dx :: Double
dx = (Double
b forall a. Num a => a -> a -> a
- Double
a) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
                   in [Double
a,Double
aforall a. Num a => a -> a -> a
+Double
dx..Double
b]
    | Double
a forall v. (InnerSpace v, Scalar v ~ Double) => v -> v -> Bool
~~ Double
b     = [forall v. (VectorSpace v, Scalar v ~ Double) => v -> v -> v
ave Double
a Double
b]
    | Bool
otherwise  = forall a. HasCallStack => [Char] -> a
error forall a b. (a -> b) -> a -> b
$ [Char]
"linSpaced:  lower limit greater than upper limit:  (n,a,b) = " forall a. [a] -> [a] -> [a]
++ 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 = forall v s. (InnerSpace v, s ~ Scalar v, Floating s) => v -> s
magnitude (v
a forall v. AdditiveGroup v => v -> v -> v
^-^ v
b) 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 forall v. AdditiveGroup v => v -> v -> v
^+^ v
v2) 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 forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double, Double) -> Position
f) Double
sl Double
su Double -> Double
tl Double -> Double
tu