diff --git a/Attrac.cabal b/Attrac.cabal --- a/Attrac.cabal +++ b/Attrac.cabal @@ -1,5 +1,5 @@ name: Attrac -version: 0.1.0 +version: 0.1.1 author: Ruben Henner Zilibowitz maintainer: Ruben Henner Zilibowitz homepage: http://patch-tag.com/r/rhz/StrangeAttractors diff --git a/Attrac.hs b/Attrac.hs --- a/Attrac.hs +++ b/Attrac.hs @@ -34,7 +34,8 @@ b = recip (sqrt (a12*a12 + a22*a22 + a32*a32)) c = recip (sqrt (a13*a13 + a23*a23 + a33*a33)) -fm cs (p,m) = (f cs p,normCols ((fjac cs p) <*> m)) +--fm cs (p,m) = (f cs p,normCols ((fjac cs p) <*> m)) +fm cs (p,m) = (f cs p,normCols $ (fjac cs p) <*> m) (<*>) :: Mat -> Mat -> Mat (<*>) (Mat a11 a12 a13 a21 a22 a23 a31 a32 a33) (Mat b11 b12 b13 b21 b22 b23 b31 b32 b33) @@ -148,3 +149,6 @@ peturb :: Value -> Pt -> Pt peturb e (Pt x y z) = Pt (x+e) y z + +(<:*>) :: Pt -> Pt -> Pt +(<:*>) (Pt ax ay az) (Pt bx by bz) = Pt (ay*bz - az*by) (bx*az - ax*bz) (ax*by - ay*bx) diff --git a/AttracCloudView.hs b/AttracCloudView.hs --- a/AttracCloudView.hs +++ b/AttracCloudView.hs @@ -32,13 +32,25 @@ coeffs :: IORef Coeffs, screen :: IORef (Int,Int), sprottCode :: IORef Int, - curTransformMat :: IORef PMat, + modelViewMat :: IORef PMat, + projectionMat :: IORef PMat, curBuffer :: IORef (Array (Int,Int) (Value,Integer)), framesElapsed :: IORef Integer } +-- position of light source lightPt = Pt 10 10 (-10) -eyePt = Pt 0 0 (-10) +{- +-- (constant,linear,quadratic) factors +lightAttenuationConstant, lightAttenuationLinear, lightAttenuationQuadratic :: Value +lightAttenuationConstant = 1 +lightAttenuationLinear = 0 +lightAttenuationQuadratic = 1 +-} + +-- initial displacement of object +locPt = Pt 0 0 10 + findOrbit :: Coeffs -> (Pt,Pt) findOrbit cs = (iterate (ffm cs) (almostZeroPt,Pt 1 0 0)) !! 128 @@ -56,7 +68,8 @@ coeffs_ <- newIORef cs screen_ <- newIORef dims sprottCode_ <- newIORef 0 - curTransformMat_ <- newIORef idPMat + modelViewMat_ <- newIORef idPMat + projectionMat_ <- newIORef (frustumMat (-1) 1 (-1) 1 5 60) curBuffer_ <- newIORef (listArray ((0,0),(width-1,height-1)) (repeat (0,0))) framesElapsed_ <- newIORef 0 return $ State { leftMouseButton = leftMouseButton_, @@ -68,24 +81,39 @@ coeffs = coeffs_, screen = screen_, sprottCode = sprottCode_, - curTransformMat = curTransformMat_, + modelViewMat = modelViewMat_, + projectionMat = projectionMat_, curBuffer = curBuffer_, framesElapsed = framesElapsed_ } +translationMat :: Pt -> PMat +translationMat (Pt x y z) = PMat + 1 0 0 x + 0 1 0 y + 0 0 1 z + 0 0 0 1 + +frustumMat :: Value -> Value -> Value -> Value -> Value -> Value -> PMat +frustumMat left right bottom top zNear zFar = PMat + (2*zNear/(right-left)) 0 a 0 + 0 (2*zNear/(top-bottom)) b 0 + 0 0 c d + 0 0 (-1) 0 + where + a = (right + left) / (right - left) + b = (top + bottom) / (top - bottom) + c = negate (zFar + zNear) / (zFar - zNear) + d = negate (2 * zFar * zNear) / (zFar - zNear) + --- --- display callback --- iterations = 1000 -iterateFuncFor 0 _ z = ([],[],z) -iterateFuncFor n f z = (a:w,b:x,y) - where (a,b) = f z - (w,x,y) = iterateFuncFor (n-1) f (a,b) - shade :: Pt -> Pt -> Value shade p@(Pt x y z) t@(Pt tx ty tz) - | ((lightPt <-> p) <.> (eyePt <-> p) < 0) = ka + | ((lightPt <-> p) <.> ((Pt 0 0 0) <-> p) < 0) = ka | gamma < 0 = error "sqrt (-1)" | diffuse < 0 = error "neg diffuse" | specular < 0 = error "neg specular" @@ -102,34 +130,46 @@ n_ = (alpha <#> l) <+> (beta <#> t) n = if (l <.> n_ < 0) then (-1) <#> n_ else n_ r = ((2 * (l <.> n)) <#> n) <-> l - v = normalise (eyePt <-> p) + v = normalise ((Pt 0 0 0) <-> p) diffuse = l <.> n shininess = 8 specular = (r <.> v) ^ shininess ka = 0.09 ks = 0.4 - kd = 1 - ka - ks + kd = 1 - ka - ks{- + pp = p <.> p + root_pp = sqrt pp + attenuation = 40 / (lightAttenuationConstant + lightAttenuationLinear*root_pp + lightAttenuationQuadratic*pp)-} -transform :: (Int,Int) -> PMat -> Pt -> (Int,Int) -transform (width,height) m p = (round $ (fromIntegral width) * 0.5 + 100*x, round $ (fromIntegral height) * 0.5 + 100*y) - where (Pt x y _) = m <*#> p +transform :: (Int,Int) -> Pt -> (Int,Int) +transform (width,height) (Pt x y _) = (round $ (fromIntegral width) * 0.5 + 100*x, round $ (fromIntegral height) * 0.5 + 100*y) colorSpeed = 1e-2 +(<*#%>) :: PMat -> Pt -> Pt +(<*#%>) m t = normalise ((m <*#> t) <-> (m <*#> zeroPt)) + display :: State -> IO () display state = do p <- readIORef (curP state) cs <- readIORef (coeffs state) - let (ps,ts,q) = iterateFuncFor iterations (ffm cs) p - m <- readIORef (curTransformMat state) + let (ps,q:_) = splitAt iterations (iterate (ffm cs) p) + {-let q = p + g <- newStdGen + let rs = take iterations (randomRs (-2,2.0) g) + let ps = [let u = Pt x (y/(sqrt$y*y+z*z)) (z/(sqrt$y*y+z*z)) in (u,(Pt 0 (-z/(sqrt$y*y+z*z)) (y/(sqrt$y*y+z*z)))) | Pt x y z <- zipWith3 (\a b c -> (Pt a b c)) rs (tail rs) (tail$tail rs)]-} + vm1 <- readIORef (modelViewMat state) + pm <- readIORef (projectionMat state) + let vm = (translationMat locPt) <*%> vm1 + let pvm = pm <*%> vm dims@(width,height) <- readIORef (screen state) let xs = filter (\((px,py),_) -> 0<=px&&px p),shade (vm <*#> p) (vm <*#%> t)) | (p,t) <- ps] modifyIORef (curBuffer state) (flip (Array.accum (\(y,c) x -> (x+y,c+1))) xs) buf <- readIORef (curBuffer state) renderPrimitive Points$sequence_ [let (x,c) = buf!(px,py) in let SRGB.RGB r g b = hsl (colorSpeed * fromIntegral c) 1 (x / (fromIntegral c)) - in (color$Color3 r g b) >> + in (color$Color4 r g b 1) >> (vertex$Vertex2 (fromIntegral px :: GLint) (fromIntegral py :: GLint)) | ((px,py),_) <- xs] writeIORef (curP state) q @@ -164,11 +204,11 @@ keyboardMouse :: State -> KeyboardMouseCallback keyboardMouse state (Char c) Down _ _ = case c of - 'm' -> print =<< (readIORef $ curTransformMat state) + 'm' -> print =<< (readIORef $ modelViewMat state) 'p' -> print =<< (readIORef $ curP state) 'f' -> print =<< (readIORef $ framesElapsed state) 'r' -> postRedisplay Nothing - ' ' -> writeIORef (curTransformMat state) idPMat >> clearAll state + ' ' -> writeIORef (modelViewMat state) idPMat >> clearAll state 'c' -> do g <- newStdGen let (cs,ly) = randomiseCoeffs g @@ -180,11 +220,20 @@ n <- readIORef (sprottCode state) let n' = (n + 1) `mod` (length codes) let cs = listArray (0,numCoeffs-1) (readCode (codes !! n')) - putStrLn $ "\nSprott code: " ++ codes!!n' ++ "\nLyapunov Exponent: " ++ show (maxLyapunovExponent cs) + putStrLn $ "\nCode index: " ++ show n' ++ "\nSprott code: " ++ codes!!n' ++ "\nLyapunov Exponent: " ++ show (maxLyapunovExponent cs) writeIORef (coeffs state) cs writeIORef (curP state) (findOrbit cs) writeIORef (sprottCode state) n' clearAll state + 'e' -> do + putStr "\nEnter Sprott code: " + hFlush stdout + code <- getLine + let cs = listArray (0,numCoeffs-1) (readCode code) + putStrLn $ "\nLyapunov Exponent: " ++ show (maxLyapunovExponent cs) + writeIORef (coeffs state) cs + writeIORef (curP state) (findOrbit cs) + clearAll state '\27' -> putStrLn "" >> exitWith ExitSuccess _ -> return () keyboardMouse state (MouseButton LeftButton) buttonState _ pos = @@ -220,19 +269,19 @@ mmb <- readIORef (middleMouseButton state) rmb <- readIORef (rightMouseButton state) case lmb of - Down -> modifyIORef (curTransformMat state) ((translateXYMat dp) <*%>) + Down -> modifyIORef (modelViewMat state) ((translateXYMat dp) <*%>) Up -> return () case mmb of - Down -> modifyIORef (curTransformMat state) ((scaleMat dp) <*%>) + Down -> modifyIORef (modelViewMat state) ((scaleMat dp) <*%>) Up -> return () case rmb of - Down -> modifyIORef (curTransformMat state) ((rotateXYMat dp) <*%>) + Down -> modifyIORef (modelViewMat state) ((rotateXYMat dp) <*%>) Up -> return () translateXYMat :: Position -> PMat translateXYMat (Position x y) = PMat - 1 0 0 (fromIntegral x * 0.01) - 0 1 0 (fromIntegral y * (-0.01)) + 1 0 0 (fromIntegral x * (-0.01)) + 0 1 0 (fromIntegral y * (0.01)) 0 0 1 0 0 0 0 1 @@ -253,10 +302,10 @@ 0 1 0 0 (-sy) 0 cy 0 0 0 0 1) - where cx = cos $ 0.001 * fromIntegral (negate y) - sx = sin $ 0.001 * fromIntegral (negate y) - cy = cos $ 0.001 * fromIntegral (negate x) - sy = sin $ 0.001 * fromIntegral (negate x) + where cx = cos $ 0.001 * fromIntegral (y) + sx = sin $ 0.001 * fromIntegral (y) + cy = cos $ 0.001 * fromIntegral (x) + sy = sin $ 0.001 * fromIntegral (x) --- --- main