二重振り子のシミュレーション

Main.hs

{-# LANGUAGE DataKinds #-}
{-# LANGUAGE ViewPatterns #-}

module Main where

import Data.Maybe
import Graphics.Gloss
import Numeric.Hamilton
import Numeric.LinearAlgebra.Static hiding ((<>))
import qualified Data.Vector.Sized as V


pendulum :: System 4 2
pendulum = mkSystem' mass transform potential
  where
    mass = vec4 1 1 1 1
    transform (V.toList -> [t1, t2]) =
      let (x1, y1) = (sin t1, -cos t1)
          (x2, y2) = (sin t1 + sin t2, -(cos t1 + cos t2))
       in fromJust $ V.fromList [x1, y1, x2, y2]
    potential (V.toList -> [_, y1, _, y2]) = y1 + y2

type Model = Phase 2

draw :: Model -> Picture
draw (Phs qs _) =
  let radian   = 180 / pi
      t1       = realToFrac $ qs <.> vec2 1 0
      t2       = realToFrac $ qs <.> vec2 0 1
      l        = 100                                    -- 振り子の見た目の長さ
      stick  = translate 0 (-l/2)  $ rectangleSolid 1 l -- 棒
      weight = translate 0 (-l)    $ circleSolid 10     -- 重り
      (x1, y1) = (l *  sin t1, l * (-cos t1))           -- １つ目の振り子の先端位置
      pendulum1 =                   rotate (-t1 * radian) $ stick <> weight
      pendulum2 = translate x1 y1 . rotate (-t2 * radian) $ stick <> weight
   in pendulum1 <> pendulum2

main :: IO ()
main = simulate inWindow white 24 initModel draw (const step)
  where
    inWindow  = InWindow "Haskell Day 2018" (640, 480) (100, 100)
    initModel = Phs (vec2 2 2) (vec2 0 0)
    step dt = stepHam (realToFrac dt) pendulum