OpenGL で RANDU

FortranでOpenGL

F03GL + FreeGlut

Fortran でグラフィックスをやる方法として OpenGL を利用するという手があります。Intel Visual Fortran のサンプルにも例題がありますが、Windows プログラミングも兼ねているので若干敷居が高い所があります。またプラットフォーム依存性が出てしまうという問題もあります。

より簡便な OpenGL 利用法として GLUT というのが昔からあります。Fortran2003ではC言語との連携が強化されたため、GLUT 利用も以前より楽になったようです。

ここでは、ここの記事からノウハウを拝借して、F03GL + FreeGlut の組み合わせで OpenGL を利用して三次元のグラフをプロットすることにします。

乱れてない乱数 RANDU

コンピュータでの乱数生成は、大型計算機では専用ハードウェアで熱雑音などで行われたりしますが、多くの場合は特殊な漸化式による疑似乱数列が用いられます。疑似乱数列は決定論的に定まるので本当は乱数では無いのですが、利用する状況に応じて十分乱数とみなせれば良いとされます。

しかしながら、一見乱数と思われたものに明らかな規則性が現れる場合もあります。その有名な例が、1960年代に IBM で利用されていた RANDU 疑似乱数列($r_0=1,r_{n+1}=65539 * r_n\, mod\,\, 2^{31}$)で、今に至るまで教科書などによく出てきます。

実行結果

RANDU により生成された数列を三次元の点座標（ $x_n=r_{3n+1}/2^{31},$ $y_n=r_{3n+2}/2^{31},$ $z_n=r_{3n+3}/2^{31}$）としてプロットしたものを以下に示します。各点は ｘ、ｙ、ｚ方向に [0,1] の範囲の立方体の中におさまっています。一見ランダムに点が打たれているように見えますが、特定角度から見ると幾枚かの平面上に点が分布しているのが見て取れます。

キャプチャ例

プログラム

以下のプログラムは F03GL のサンプル中の pointsz.f90 を土台に、いくらか書き直したものです。カーソルによって視点方向を動かすことができます。

module m_gl  
  use OpenGL_GL
  use OpenGL_GLUT
  implicit none  
  real, save :: xrot = 0.0, yrot = 0.0
contains
  subroutine SetupRC() bind(c)
    ! Black background
    call glClearColor(0.0, 0.0, 0.6, 1.0 )
    ! Set drawing color to green
    call glColor3f(1.0, 0.0, 0.0)
  end subroutine SetupRC

  subroutine KeyPressFunc(key, x, y) bind(C)
    integer(GLbyte), intent(in), value :: key
    integer(GLint ), intent(in), value  :: x, y
    if ( key == 27 ) stop ! ESC
  end subroutine KeyPressFunc

  subroutine KeySpecialFunc(key, x, y) bind(C)
    integer(GLint), intent(in), value  :: key, x, y
    if ( key == GLUT_KEY_UP    ) xRot = xrot - 5.0
    if ( key == GLUT_KEY_DOWN  ) xRot = xrot + 5.0
    if ( key == GLUT_KEY_LEFT  ) yRot = yrot - 5.0
    if ( key == GLUT_KEY_RIGHT ) yRot = yrot + 5.0 
    if ( xrot > 356.0 )then
      xRot = 0.0
    else if ( xrot < -1.0 )then
      xRot = 355.0
    endif
    if ( yrot > 356.0 )then
      yRot = 0.0
    else if ( yrot < -1.0 )then
      yRot = 355.0
    endif
    ! Refresh the Window
    call glutPostRedisplay
  end subroutine KeySpecialFunc

  subroutine ChangeSize(win, hin) bind(C)
    integer(GLcint), intent(IN), value :: win, hin
    integer(GLcint) :: w, h
    real(GLdouble)  :: Zero, One, Range = 100.0, Aspect
    w = win
    h = hin
    ! Prevent a divide by zero, when window is too short
    ! (you cant make a window of zero width).
    if( h == 0 ) h = 1
    ! Set the viewport to be the entire window
    call glViewport(0, 0, w, h)
    ! Reset coordinate system
    call glMatrixMode(GL_PROJECTION)
    call glLoadIdentity
    ! Establish clipping volume (left, right, bottom, top, near, far)
    aspect = float(w)/float(h)
    ! Keep the square square
    if( w <= h )then
      call glOrtho( -Range, Range, &
                    -Range / Aspect, Range / Aspect,&
                    -Range, Range  )
    else
      call glOrtho( -Range * Aspect, Range * Aspect, &
                    -Range, Range, &
                    -Range, Range  )
    endif
    call glMatrixMode(GL_MODELVIEW)
    call glLoadIdentity
  end subroutine ChangeSize 

  subroutine RenderScene() bind(C)
    integer, parameter :: n = 10**4
    real(glfloat), save :: dx(n - 1), dy(n - 1), dz(n - 1)
    logical, save :: first = .true.
    integer :: i
    ! make Lorenz attractor data        

    if (first) then
      block  
        integer(8) :: m = 2_8**31
        integer(8) :: ix 
        real(glfloat) :: x(n), y(n), z(n)
        first = .false.
        ix = 1
        do i = 1, n
          ix = mod(65539 * ix, m)  ! ;  print *, ix
          x(i) = real(ix) / m 
          ix = mod(65539 * ix, m)  ! ;  print *, ix
          y(i) = real(ix) / m
          ix = mod(65539 * ix, m)  ! ;  print *, ix
          z(i) = real(ix) / m
        end do
        dx = (x(2:) - x(:n - 1)) * 80.0_glfloat ! scale data 
        dy = (y(2:) - y(:n - 1)) * 80.0_glfloat
        dz = (z(2:) - z(:n - 1)) * 80.0_glfloat
      end block  
    end if  
    ! Clear the window with current clearing color
    call glClear(GL_COLOR_BUFFER_BIT)
    ! Save matrix state and do the rotation
    call glPushMatrix
    call glRotatef(xRot, 1.0, 0.0, 0.0)
    call glRotatef(yRot, 0.0, 1.0, 0.0)
    ! draw Lorenz attractor
    do i = 1, n - 1
      call glTranslatef( dx(i), dy(i), dz(i) )
      call glutSolidCube(0.5_gldouble )
    end do
    !
    call glPopMatrix
    ! Flush drawing commands
    call glutSwapBuffers
  end subroutine RenderScene
end module m_gl

program RANDU
   use m_gl
   implicit none
   integer :: iwin
   call glutInit()
   call glutInitDisplayMode(ior(GLUT_DOUBLE,ior(GLUT_RGB,GLUT_DEPTH)))
   call glutInitWindowSize(800, 600)
   iwin = glutCreateWindow('RANDU'//char(0))
   call glutReshapeFunc ( ChangeSize     )
   call glutKeyboardFunc( KeyPressFunc   )
   call glutSpecialFunc ( KeySpecialFunc )
   call glutDisplayFunc ( RenderScene    )
   call SetupRC()
   call glutMainLoop()
end program RANDU