{-# LANGUAGE CPP #-}
-- Quasicrystals demo.
--
-- Based on code from:
-- http://hackage.haskell.org/package/repa-examples
-- http://mainisusuallyafunction.blogspot.com/2011/10/quasicrystals-as-sums-of-waves-in-plane.html
--
module Main where
import Prelude as P
import Config
import Monitoring
import ParseArgs
import Data.Label ( get )
import Criterion.Main ( defaultMainWith, bench, whnf )
import System.Environment
import Data.Array.Accelerate as A hiding ( size )
import Graphics.Gloss.Accelerate.Raster.Field as G
-- Types ----------------------------------------------------------------------
-- | Angle in radians.
type Angle = Float
-- | Angle offset used for animation.
type Phi = Float
-- | Number of waves to sum for each pixel.
type Degree = Int
-- | Feature size of visualisation.
type Scale = Float
-- | Time in seconds since the program started.
type Time = Float
-- Point ----------------------------------------------------------------------
-- | Compute a single point of the visualisation.
quasicrystal :: Scale -> Degree -> Exp Time -> Exp Point -> Exp Color
quasicrystal scale degree time p
= let -- Scale the time to be the phi value of the animation.
-- The action seems to slow down at increasing phi values, so we
-- increase phi faster as time moves on.
phi = 1 + (time ** 1.5) * 0.005
in rampColour
$ waves degree phi
$ point scale p
-- | Sum up all the waves at a particular point.
waves :: Degree -> Exp Phi -> Exp Point -> Exp Float
waves degree phi x = wrap $ waver degree 0
where
th = pi / phi
waver :: Int -> Exp Float -> Exp Float
waver n acc
| n == 0 = acc
| otherwise = waver (n - 1) (acc + wave (A.constant (P.fromIntegral n) * th) x)
wrap n
= let n_ = A.truncate n :: Exp Int
n' = n - A.fromIntegral n_
in
A.odd n_ A.? (1-n', n')
-- | Generate the value for a single wave.
wave :: Exp Angle -> Exp Point -> Exp Float
wave th pt = (cos (cth*x + sth*y) + 1) / 2
where
(x, y) = xyOfPoint pt
cth = cos th
sth = sin th
-- | Convert an image point to a point on our wave plane.
point :: Scale -> Exp Point -> Exp Point
point scale pt = makePoint (x * scale') (y * scale')
where
(x, y) = xyOfPoint pt
scale' = A.constant scale
-- | Colour ramp from red to white, convert into RGBA
rampColour :: Exp Float -> Exp Color
rampColour v = rawColor 1 (0.4 + (v * 0.6)) v 1
-- Main -----------------------------------------------------------------------
main :: IO ()
main = do
beginMonitoring
argv <- getArgs
(config, crit, nops) <- parseArgs optHelp optBackend options defaults header footer argv
let size = get optSize config
zoom = get optZoom config
scale = get optScale config
degree = get optDegree config
backend = get optBackend config
-- for benchmarking
force arr = A.indexArray arr (Z:.0:.0) `seq` arr
frame = run1 backend
$ makeField size size (\time -> quasicrystal scale degree (the time))
if get optBench config
then withArgs nops $ defaultMainWith crit (return ())
[ bench "crystal" $ whnf (force . frame) (A.fromList Z [1.0]) ]
else G.animateFieldWith
(run1 backend)
(InWindow "Quasicrystals" (size * zoom, size * zoom) (10, 10))
(zoom, zoom)
(quasicrystal scale degree)