{-# LANGUAGE Rank2Types, TypeOperators, FlexibleContexts, TypeFamilies #-}
{-# OPTIONS_GHC -Wall #-}
----------------------------------------------------------------------
-- |
-- Module      :  Graphics.FieldTrip.ParamSurf
-- Copyright   :  (c) Conal Elliott 2008
-- License     :  BSD3
-- 
-- Maintainer  :  conal@conal.net
-- Stability   :  experimental
-- 
-- Parametric surfaces with automatic normals
----------------------------------------------------------------------

module Graphics.FieldTrip.ParamSurf where

import Control.Applicative

import Data.NumInstances ()
import Data.VectorSpace
import Data.Cross


type HeightField s = Two s -> One s
type Surf        s = Two s -> Three s

type USurf = forall s. Floating s => Surf s


type Curve2 s = One s -> Two s
type Curve3 s = One s -> Three s

type Warp1 s = One   s -> One   s
type Warp2 s = Two   s -> Two   s
type Warp3 s = Three s -> Three s

mul2pi :: Floating s => s -> s
mul2pi = (* (2*pi))

-- | Trig functions with unit period ([-1/2,1/2])
cosU, sinU :: Floating s => s -> s
cosU = cos . mul2pi
sinU = sin . mul2pi


-- | Turn a height field into a surface
hfSurf :: HeightField s -> Surf s
hfSurf field = \ (u,v) -> (u, v, field (u,v))

-- | Like 'hfSurf' but for curve construction
fcurve :: Warp1 s -> Curve2 s
fcurve f = \ u -> (u, f u)

-- | Unit circle.
circle :: Floating s => Curve2 s
circle = liftA2 (,) cosU sinU

-- | Half semi circle, with theta in [-pi/2,pi/2]
semiCircle :: Floating s => Curve2 s
semiCircle = circle . (/ 2)

-- | Torus, given radius of sweep circle and cross section
torus :: (Floating s, VectorSpace s, Scalar s ~ s) => s -> s -> Surf s
-- torus sr cr = revolve (\ s -> (sr,0) ^+^ cr *^ circle s)
torus sr cr = revolve (const (sr,0) ^+^ const cr *^ circle)

-- Surface of revolution, formed by rotation around Z axis.  The curve is
-- parameterized by u, and the rotation by v.  In this generalized
-- version, we have not a single curve, but a function from v to curves.
revolveG :: Floating s => (s -> Curve2 s) -> Surf s
revolveG curveF = \ (u,v) -> onXY (rotate (-2*pi*v)) (addY (curveF v) u)

revolve :: Floating s => Curve2 s -> Surf s
revolve curve = revolveG (const curve)

-- A sphere is a revolved semi-circle
sphere1 :: Floating s => Surf s
sphere1 = revolve semiCircle


-- | Profile product.
profile :: Num s => Curve2 s -> Curve2 s -> Surf s
profile curve prof (u,v) = (cx*px,cy*px,py)
 where
   (cx,cy) = curve u
   (px,py) = prof  v

-- More spheres
sphere2,sphere3 :: Floating s => Surf s
sphere2 = profile circle semiCircle
sphere3 = profile semiCircle circle

-- | Frustum, given base & cap radii and height.
frustum :: (Floating s, VectorSpace s, Scalar s ~ s) => s -> s -> s -> Surf s
frustum baseR topR h = profile circle rad
 where
   rad t = (lerp baseR topR (t + 1/2), h*t)

-- | Unit cylinder.  Unit height and radii
ucylinder :: (Floating s, VectorSpace s) => Surf s
ucylinder = profile circle (const 1)


-- | Given a combining op and two curves, make a surface.  A sort of
-- Cartesian product with combination.
cartF :: (a -> b -> c) -> (u -> a) -> (v -> b) -> ((u,v) -> c)
cartF op f g = \ (u,v) -> f u `op` g v

-- Sweep a basis curve by a sweep curve.  Warning: does not reorient the
-- basis curve as cross-section.  TODO: Frenet frame.
sweep :: VectorSpace s => Curve3 s -> Curve3 s -> Surf s
sweep = cartF (^+^)


-- | One period, unit height eggcrate
eggcrateH :: Floating s => HeightField s
eggcrateH = cartF (*) cosU sinU

revolveH :: (Floating s, InnerSpace s, Scalar s ~ s) => Warp1 s -> HeightField s
revolveH = (. magnitude)

rippleH :: (Floating s, InnerSpace s, Scalar s ~ s) => HeightField s
rippleH = revolveH sinU

-- | Simple ripply pond shape
ripple :: Floating s => Surf s
ripple = -- onXY' (2 *^) $
         revolve (const (0.5,0) - fcurve sinU)

-- | Apply a displacement map at a value

displaceV :: (InnerSpace v, s ~ Scalar v, Floating s, HasNormal v) =>
             v -> Scalar v -> v
displaceV v s = v ^+^ s *^ normal v

-- | Apply a displacement map to a function (e.g., 'Curve2' or 'Surf') or
-- other container.
displace :: (InnerSpace v, Scalar v ~ s, Floating s, HasNormal v, Applicative f) =>
            f v -> f (Scalar v) -> f v
displace = liftA2 displaceV



---- Misc

rotate :: Floating s => s -> Warp2 s
rotate theta = \ (x,y) -> (x * c - y * s, y * c + x * s)
 where c = cos theta
       s = sin theta

addX, addY, addZ :: Num s => (a -> Two s) -> (a -> Three s)
addX = fmap (\ (y,z) -> (0,y,z))
addY = fmap (\ (x,z) -> (x,0,z))
addZ = fmap (\ (x,y) -> (x,y,0))

addYZ,addXZ,addXY :: Num s => (a -> One s) -> (a -> Three s)
addYZ = fmap (\ x -> (x,0,0))
addXZ = fmap (\ y -> (0,y,0))
addXY = fmap (\ z -> (0,0,z))

onX,onY,onZ :: Warp1 s -> Warp3 s
onX f (x,y,z) = (f x, y, z)
onY f (x,y,z) = (x, f y, z)
onZ f (x,y,z) = (x, y, f z)

onXY,onYZ,onXZ :: Warp2 s -> Warp3 s
onXY f (x,y,z) = (x',y',z ) where (x',y') = f (x,y)
onXZ f (x,y,z) = (x',y ,z') where (x',z') = f (x,z)
onYZ f (x,y,z) = (x ,y',z') where (y',z') = f (y,z)


onX',onY',onZ' :: Warp1 s -> (a -> Three s) -> (a -> Three s)
onX' = fmap fmap onX
onY' = fmap fmap onY
onZ' = fmap fmap onZ

onXY',onXZ',onYZ' :: Warp2 s -> (a -> Three s) -> (a -> Three s)
onXY' = fmap fmap onXY
onXZ' = fmap fmap onXZ
onYZ' = fmap fmap onYZ