Haskell+OpenGLでローレンツアトラクタを描画してみた

前回の記事を公開したところ、友達にどうせならアトラクタを持ってる方程式を書けと言われたので前回に引き続き今度はローレンツ方程式を描画してみることにしました。

今回はWASDでズームと回転が出来るようにしています。

import Data.IORef
import Graphics.UI.GLUT

GLUTが必要なので入っていない場合は

$ cabal install GLUT

で予めインストールしておいてください☕

type Point = (GLfloat, GLfloat, GLfloat)
type Line  = [Point]

型宣言がごちゃごちゃするのを防ぐために点と線の型を定義しておきます

-- 初期値
start = (1.0, 1.0, 1.0)

-- ローレンツ方程式
lorenz :: Point -> Point
lorenz = \(x, y, z) ->
    let dt = 0.01
        p  = 10.0
        r  = 28.0
        b  = 8.0/3.0
    in  ( (\x' -> x + x' * dt) $ p * (y - x)
        , (\y' -> y + y' * dt) $ x * (r - z) - y
        , (\z' -> z + z' * dt) $ x * y - b * z
        )

ローレンツ方程式本体と初期値です。lorenzは点の情報を受け取ると微小時間だけ動いた次の点の情報を返します。微小時間はdt = 0.01と決めで実装してますが引数で取るように修正してもいいでしょう。

main :: IO ()
main = do
    -- OpenGLの初期化
    (_progName, _args) <- getArgsAndInitialize
    initialDisplayMode $= [DoubleBuffered]
    _window <- createWindow "Lorenz Attractor"
    -- 共有パラメータの設定
    rot <- newIORef 0
    sc  <- newIORef 0.02
    orbit <- newIORef [start]
    -- コールバックの設定
    keyboardMouseCallback $= Just (keyboardMouse rot sc)
    idleCallback $= Just (idle orbit)
    displayCallback $= display orbit rot sc
    -- 実行
    mainLoop

共有パラメータはrotは回転角をscは拡大率をorbitは計算した軌道を保存しておくのに使います。

idle :: IORef Line -> IdleCallback
idle orbit = do
    ps <- get orbit
    orbit $= (lorenz (head ps) : take 10000 ps)
    postRedisplay Nothing

待機時に実行される関数でローレンツ方程式の計算をしています。計算した値はorbitに保存しています。

display :: IORef Line -> IORef GLfloat -> IORef GLfloat -> DisplayCallback
display orbit rot sc = do 
    -- バッファをクリア
    clear [ColorBuffer]
    -- 共有パラメータの取得
    r <- get rot
    s <- get sc
    ps <- get orbit
    -- 変換行列の初期化
    loadIdentity
    preservingMatrix $ do
        scale s s s
        rotate r $ Vector3 0 1 0
        -- 軌道の描画
        let mkVertex (x, y, z) = vertex $ Vertex3 x y z
        renderPrimitive LineStrip $ mapM_ mkVertex ps
    -- 表示
    swapBuffers

軌道を実際に描画する部分です。preservingMatrixはとても便利で行列のpushとpopと勝手にやってくれて現在の状態を切り出して中で自由に使うことができます。

keyboardMouse :: IORef GLfloat -> IORef GLfloat -> KeyboardMouseCallback
keyboardMouse rot sc _key _state _ _ = do
    case _key of
        Char 'w' -> sc  $~! (*1.1)
        Char 'a' -> rot $~! (+5)
        Char 's' -> sc  $~! (*0.9)
        Char 'd' -> rot $~! (subtract 5)
        _        -> return ()

キーボード（とマウス）のイベントを制御する関数です。WASDが押されると対応するパラメータの値を変更しています。

実行するときは必ずコンパイルしてから実行してください。

$ ghc Main.hs
$ ./Main

描画されたグラフを眺めていると手前と奥が表示されなくなっていることに気が付きました。調べるとカメラの設定でどこまで表示するかを調整できるようでfrustumとかperspevtiveとかをいじってたのですが結局直すことができませんでした… OpenGLの勉強を真面目にしたいと思います。

コード全体を以下に載せておきます。

---

import Data.IORef
import Graphics.UI.GLUT

type Point = (GLfloat, GLfloat, GLfloat)
type Line  = [Point]

-- 初期値
start = (1.0, 1.0, 1.0)

-- ローレンツ方程式
lorenz :: Point -> Point
lorenz = \(x, y, z) ->
    let dt = 0.01
        p  = 10.0
        r  = 28.0
        b  = 8.0/3.0
    in  ( (\x' -> x + x' * dt) $ p * (y - x)
        , (\y' -> y + y' * dt) $ x * (r - z) - y
        , (\z' -> z + z' * dt) $ x * y - b * z
        )

main :: IO ()
main = do
    -- OpenGLの初期化
    (_progName, _args) <- getArgsAndInitialize
    initialDisplayMode $= [DoubleBuffered]
    _window <- createWindow "Lorenz Attractor"
    -- 共有パラメータの設定
    rot <- newIORef 0
    sc  <- newIORef 0.02
    orbit <- newIORef [start]
    -- コールバックの設定
    keyboardMouseCallback $= Just (keyboardMouse rot sc)
    idleCallback $= Just (idle orbit)
    displayCallback $= display orbit rot sc
    -- 実行
    mainLoop

idle :: IORef Line -> IdleCallback
idle orbit = do
    ps <- get orbit
    orbit $= (lorenz (head ps) : take 10000 ps)
    postRedisplay Nothing

display :: IORef Line -> IORef GLfloat -> IORef GLfloat -> DisplayCallback
display orbit rot sc = do 
    -- バッファをクリア
    clear [ColorBuffer]
    -- 共有パラメータの取得
    r <- get rot
    s <- get sc
    ps <- get orbit
    -- 変換行列の初期化
    loadIdentity
    preservingMatrix $ do
        scale s s s
        rotate r $ Vector3 0 1 0
        -- 軌道の描画
        let mkVertex (x, y, z) = vertex $ Vertex3 x y z
        renderPrimitive LineStrip $ mapM_ mkVertex ps
    -- 表示
    swapBuffers

keyboardMouse :: IORef GLfloat -> IORef GLfloat -> KeyboardMouseCallback
keyboardMouse rot sc _key _state _ _ = do
    case _key of
        Char 'w' -> sc  $~! (*1.1)
        Char 'a' -> rot $~! (+5)
        Char 's' -> sc  $~! (*0.9)
        Char 'd' -> rot $~! (subtract 5)
        _        -> return ()