fortran66のブログ

fortran について書きます。

Rutherford Scattering 2

散乱に関するラザフォードの公式では、クーロン場の符号(斥力・引力)に依らずに同じ角度依存性が得られるわけですが、それについて計算結果を見てみることにします。
対称性に鑑みて、散乱体を含む二次元平面で切った計算をすることにします。散乱ポテンシャルを与える散乱体は固定されているとします。

計算結果

赤は引力ポテンシャルによる散乱、青は斥力ポテンシャルによる散乱です。クーロン・ポテンシャルの大きさは同じです。青の斥力の場合、前回も見たように斥力を受けるので、散乱体を避ける形で散乱体の後ろ側に散乱された粒子がこない領域があります。一方、赤色の引力ポテンシャルの場合、ポテンシャル中心に引き寄せられる分、散乱体の後ろ側に多く散乱されているように見えます。

同じ結果を横から見たもの。平面内での散乱だと分かります。

しかし、散乱体の中心近くを通る散乱をスケールを変えて遠くから見ると、斥力ポテンシャルと引力ポテンシャルの散乱はとてもよく似た形になります。引力の場合、いわば散乱体の側を通ったものが引力をうけて散乱体の裏側を通り抜けて反対側に抜けてゆく結果、斥力の場合と似た形になっているように見受けられます。

裏から見た図。

なお、引力ポテンシャルの場合に粒子がいくつか激しく散乱されていますが、これはポテンシャル中心に近づきすぎてオイラー法の適用範囲を外れ、おかしな結果を与えているものではないかと思います。本来は、入射方向へぐるりと一回転するような軌跡を描くべきものではないかと想像します。

散乱体から離れたところを通る(散乱パラメータの大きな)入射粒子も考慮した結果を以下に示します。ほとんどの粒子は、ほぼ直進して通り抜けています。


軌跡の出発座標と終端座標から、散乱角に関する分布をまとめて、ラザフォードの公式と比較すると面白いのではないかと思います。

スケールが無いので、まともな図になってません。自明な剛体球からはじめるべきだった(><;

Fortran ソース

!=====================================================
MODULE m_scatter
IMPLICIT NONE
PRIVATE :: dt, alpha_prime
REAL :: dt = 0.1
REAL :: alpha_prime = 1.0 ! Coulomb Potential parameter
!
TYPE :: t_vec
 REAL :: x, y, z
END TYPE
!
TYPE :: t_particle
 TYPE (t_vec) :: pos, vel
END TYPE
!
INTERFACE OPERATOR (+)
 MODULE PROCEDURE vec_plus
END INTERFACE
!
INTERFACE OPERATOR (*)
 MODULE PROCEDURE vec_a_mult
END INTERFACE
!
CONTAINS
!---------------------------------
SUBROUTINE set_time_interval(t)
REAL, INTENT(IN) :: t
dt = t
RETURN
END SUBROUTINE set_time_interval
!---------------------------------
SUBROUTINE set_Coulomb_parameter(a)
REAL, INTENT(IN) :: a
alpha_prime = a
RETURN
END SUBROUTINE set_Coulomb_parameter
!---------------------------------
TYPE (t_vec) FUNCTION vec_plus(v1, v2)
TYPE (t_vec), INTENT(IN) :: v1, v2
vec_plus%x = v1%x + v2%x
vec_plus%y = v1%y + v2%y
vec_plus%z = v1%z + v2%z
RETURN
END FUNCTION vec_plus
!---------------------------------
TYPE (t_vec) FUNCTION vec_a_mult(a, v)
REAL, INTENT(IN) :: a
TYPE (t_vec), INTENT(IN) :: v
vec_a_mult%x = a * v%x
vec_a_mult%y = a * v%y
vec_a_mult%z = a * v%z
RETURN
END FUNCTION vec_a_mult
!---------------------------------
TYPE (t_particle) FUNCTION init_particle(position, velocity)
TYPE (t_vec), INTENT(IN) :: position, velocity
init_particle = t_particle( position, velocity )
RETURN
END FUNCTION init_particle
!---------------------------------
REAL FUNCTION radius(v)
TYPE (t_vec), INTENT(IN) :: v
radius = SQRT(v%x ** 2 + v%y **2 + v%z ** 2)
RETURN
END FUNCTION radius
!---------------------------------
TYPE (t_vec) FUNCTION force(position)
TYPE (t_vec), INTENT(IN) :: position
REAL :: r
r = radius(position)
force%x = alpha_prime * position%x / r ** 3 ! derivative of the Coulomb potential
force%y = alpha_prime * position%y / r ** 3
force%z = alpha_prime * position%z / r ** 3
RETURN
END FUNCTION force
!---------------------------------
SUBROUTINE step(p) ! Euler's method
TYPE (t_particle), INTENT(IN OUT) :: p
TYPE (t_vec) :: dv, dx
p%vel = p%vel + dt * force(p%pos)
p%pos = p%pos + dt * p%vel
RETURN
END SUBROUTINE step
!---------------------------------
END MODULE m_scatter
!=====================================================
MODULE m_ART
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
!---- TYPE definitions ---------------------------------
TYPE, BIND(C) :: t_ARParam
 INTEGER(C_INT) :: ixsize, iysize
 REAL(C_DOUBLE) :: dmat(4, 3)
 REAL(C_DOUBLE) :: dist_factor(4)
END TYPE
!
TYPE, BIND(C) :: t_ARMarkerInfo 
 INTEGER(C_INT) :: iarea
 INTEGER(C_INT) :: id
 INTEGER(C_INT) :: dir
 REAL(C_DOUBLE)  :: cf
 REAL(C_DOUBLE)  :: pos(2)
 REAL(C_DOUBLE)  :: dline(3, 4)
 REAL(C_DOUBLE)  :: vertex(2, 4)
END TYPE
!---- global variables ---------------------------------
INTEGER(C_INT) :: ithresh = 100
INTEGER(C_INT) :: icount = 0
INTEGER(C_INT) :: ipatt_id
REAL(C_DOUBLE) :: patt_width = 80.0_c_double
REAL(C_DOUBLE) :: patt_center(2) = (/0.0_c_double, 0.0_c_double/)
REAL(C_DOUBLE) :: patt_trans(3, 4)
INTEGER(C_INT) :: ixsize, iysize
TYPE(t_ARParam) :: cparam
!---- INTERFACE BLOCK -----------------------------------
INTERFACE
 SUBROUTINE arVideoCapStop() BIND(C, NAME = 'arVideoCapStop')
 END SUBROUTINE arVideoCapStop
! 
 SUBROUTINE arVideoClose() BIND(C, NAME = 'arVideoClose')
 END SUBROUTINE arVideoClose
!
 SUBROUTINE argCleanup() BIND(C, NAME = 'argCleanup')
 END SUBROUTINE argCleanup
! 
 INTEGER(C_INT) FUNCTION arVideoGetImage() BIND(C, NAME = 'arVideoGetImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE arUtilSleep(i) BIND(C, NAME = 'arUtilSleep')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: i 
 END
!
 SUBROUTINE arUtilTimerReset() BIND(C, NAME = 'arUtilTimerReset')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE argDrawMode2D() BIND(C, NAME = 'argDrawMode2D')
 END 
!
 SUBROUTINE argDispImage( iptr, i, j ) BIND(C, NAME = 'argDispImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER, VALUE :: i, j
 END 
!
 INTEGER(C_INT) FUNCTION arDetectMarker( iptr, i, mark, k ) BIND(C, NAME = 'arDetectMarker')
 USE, INTRINSIC :: ISO_C_BINDING 
 IMPORT :: t_ARMarkerInfo
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER(C_INT)     , VALUE :: i
 TYPE(t_ARMarkerInfo), POINTER :: mark(:)
 INTEGER(C_INT) :: k
 END 
!
 SUBROUTINE arVideoCapNext() BIND(C, NAME = 'arVideoCapNext')
 END 
!
 SUBROUTINE argSwapBuffers() BIND(C, NAME = 'argSwapBuffers')
 END 
!
 SUBROUTINE arGetTransMat(ptr, d1, d, d2) BIND(C, NAME = 'arGetTransMat')
 USE, INTRINSIC :: ISO_C_BINDING 
 TYPE(C_PTR), VALUE :: ptr
 REAL(C_DOUBLE), VALUE :: d
 REAL(C_DOUBLE):: d1(2), d2(3, 4)
 END 
!
 SUBROUTINE argDrawMode3D() BIND(C, NAME = 'argDrawMode3D')
 END SUBROUTINE argDrawMode3D
 !
 SUBROUTINE argDraw3dCamera(ix, iy) BIND(C, NAME = 'argDraw3dCamera')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: ix, iy 
 END
 !
 SUBROUTINE argConvGlpara(x, y) BIND(C, NAME = 'argConvGlpara')
 USE, INTRINSIC :: ISO_C_BINDING 
 REAL(C_DOUBLE) :: x(4, 3), y(4, 4) 
 END
! 
 SUBROUTINE arVideoCapStart() BIND(C, NAME = 'arVideoCapStart')
 END
!
 SUBROUTINE argMainLoop( i, j, k ) BIND(C, NAME = 'argMainLoop')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER, VALUE :: i
 TYPE(C_FUNPTR), VALUE :: j, k
 END
!
 INTEGER(C_INT) FUNCTION arVideoOpen(text) BIND(C, NAME = 'arVideoOpen')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(C_CHAR) :: text(*)
 END
!
 INTEGER(C_INT) FUNCTION arVideoInqSize(ix, iy) BIND(C, NAME = 'arVideoInqSize')
 USE, INTRINSIC :: ISO_C_BINDING
 INTEGER(C_INT) :: ix, iy
 END
! 
 INTEGER(C_INT) FUNCTION arParamLoad(text, i, xparam) BIND(C, NAME = 'arParamLoad')
 USE, INTRINSIC :: ISO_C_BINDING
 IMPORT :: t_ARParam
 CHARACTER(C_CHAR) :: text(*)
 TYPE(t_ARParam) :: xparam
 INTEGER(C_INT), VALUE :: i
 END
!
 SUBROUTINE arParamChangeSize(xparam, i, j, yparam) BIND(C, NAME = 'arParamChangeSize')
 USE, INTRINSIC :: ISO_C_BINDING
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam, yparam
 INTEGER(C_INT), VALUE :: i, j
 END
! 
 SUBROUTINE arInitCparam(xparam) BIND(C, NAME = 'arInitCparam' )
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
! 
 SUBROUTINE arParamDisp(xparam) BIND(C, NAME = 'arParamDisp')
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
!
 INTEGER(C_INT) FUNCTION arLoadPatt(text) BIND(C, NAME = 'arLoadPatt')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(C_CHAR) :: text(*)
 END
!
 SUBROUTINE argInit(iptr, d, i1, i2, i3, i4) BIND(C, NAME = 'argInit')
 USE, INTRINSIC :: ISO_C_BINDING
 TYPE(C_PTR), VALUE :: iptr
 REAL(C_DOUBLE)      , VALUE :: d
 INTEGER(C_INT)     , VALUE :: i1, i2, i3, i4
 END
!
 REAL(C_DOUBLE) FUNCTION arUtilTimer() BIND(C, NAME = 'arUtilTimer')
 USE, INTRINSIC :: ISO_C_BINDING
 END FUNCTION arUtilTimer
END INTERFACE
CONTAINS
!-------------------------------------------------
SUBROUTINE init() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
TYPE(t_ARParam) :: wparam
CHARACTER (LEN = 35) :: vconf = "Data/WDM_camera_flipV.xml" // c_null_char
CHARACTER (LEN = 35) :: cparam_name = "Data/camera_para.dat" // c_null_char
CHARACTER (LEN = 35) :: patt_name   = "Data/patt.hiro"       // c_null_char
!   /* open the video path */
IF( arVideoOpen( vconf ) < 0 ) STOP
!    /* find the size of the window */
IF( arVideoInqSize( ixsize, iysize ) < 0 ) STOP
PRINT *, "Image size (x,y) = ", ixsize, iysize
!    /* set the initial camera parameters */
IF ( arParamLoad( cparam_name, 1, wparam ) < 0 ) THEN
 PRINT *, "Camera parameter load error !!"
 STOP
END IF
CALL arParamChangeSize( wparam, ixsize, iysize, cparam )
CALL arInitCparam( cparam )
PRINT *, "*** Camera Parameter ***"
CALL arParamDisp( cparam )
ipatt_id = arLoadPatt( patt_name )
IF ( ipatt_id < 0 ) THEN
 PRINT *, "pattern load error !!"
 STOP
END IF
!    /* open the graphics window */
CALL argInit( C_LOC(cparam), 1.0_c_double, 0, 0, 0, 0 )
RETURN
END SUBROUTINE init
!-------------------------------------------------
SUBROUTINE mainLoop() BIND(C, NAME = 'mainLoop')
IMPLICIT NONE
INTEGER(C_INTPTR_T) :: dataPtr
TYPE(t_ARMarkerInfo), POINTER:: marker_info(:)
INTEGER(C_INT) :: marker_num
INTEGER(C_INT) :: j, k
LOGICAL, SAVE :: qfirst = .TRUE.
IF (qfirst) THEN 
 ALLOCATE(marker_info(50)) ! ?
 qfirst = .FALSE.
END IF
dataPtr = arVideoGetImage()
IF (dataPtr == 0) THEN
 CALL arUtilSleep(2)
 RETURN
END IF
IF (icount == 0) CALL arUtilTimerReset()
icount = icount + 1
CALL argDrawMode2D()
CALL argDispImage( dataPtr, 0, 0 )
IF ( arDetectMarker( dataPtr, ithresh, marker_info, marker_num ) < 0 ) THEN
 CALL cleanup()
 RETURN
END IF
CALL arVideoCapNext()
k = -1
DO j = 1, marker_num
 IF (ipatt_id == marker_info(j)%id) THEN
  IF (k == -1) THEN 
   k = j
  ELSE 
   IF (marker_info(k)%cf < marker_info(j)%cf ) k = j
  END IF
 END IF
END DO
IF ( k == -1 ) THEN 
 CALL argSwapBuffers()
 RETURN
END IF
! draw-main 
CALL arGetTransMat( C_LOC(marker_info(k)), patt_center, patt_width, patt_trans )
CALL draw()
CALL argSwapBuffers()
RETURN
END SUBROUTINE mainLoop
!-------------------------------------------------
SUBROUTINE cleanup() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE 
CALL arVideoCapStop()
CALL arVideoClose()
CALL argCleanup()
RETURN
END SUBROUTINE cleanup
!-------------------------------------------------
SUBROUTINE keyEvent( key, ix, iy ) BIND(C, NAME = 'keyEvent')
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
INTEGER(C_SIGNED_CHAR), VALUE :: key
INTEGER(C_INT)        , VALUE :: ix, iy
IF (key == 27) THEN
  PRINT *, "*** ", icount / arUtilTimer(), "(frame/sec)"
  CALL cleanup()
  STOP
END IF
RETURN
END SUBROUTINE keyEvent
!---------------------------------------------------------------
SUBROUTINE draw() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
USE m_scatter
IMPLICIT NONE
REAL(gldouble) :: gl_para(4, 4)
REAL(glfloat ) :: mat_ambient(4) = [ 0.0_glfloat, 1.0_glfloat, 0.0_glfloat, 1.0_glfloat ]
REAL(glfloat ) :: mat_flash(4)   = [ 0.0_glfloat, 0.0_glfloat, 0.0_glfloat, 1.0_glfloat ]
REAL(glfloat ) :: mat_flash_shiny(1) = [ 50.0 ]
REAL(glfloat ) :: ambi(4)            = [ 0.1_glfloat, 0.1_glfloat, 0.1_glfloat, 0.1_glfloat ]
REAL(glfloat ) :: lightZeroColor(4)  = [ 0.0_glfloat, 1.0_glfloat, 1.0_glfloat, 0.1_glfloat ]
! shadow  
REAL(gldouble) :: shadowMatrix(4, 4), plane(4)
REAL(glfloat ), SAVE :: xl = 100.0_glfloat, yl = 0.0_glfloat, zl = 100.0_glfloat, wl = 0.0_glfloat, theta = 0.0
!
INTEGER, SAVE :: icount
LOGICAL, SAVE :: qfirst = .TRUE.
REAL :: x0, y0, z0
!
INTEGER, PARAMETER :: max_pos = 1000, max_traject = 50000
REAL(glfloat), SAVE :: traject(max_traject, 3, max_pos, 2) 
TYPE (t_particle) :: p
INTEGER :: i, ipos
! 
CALL argDrawMode3D()
CALL argDraw3dCamera( 0, 0 )
CALL glClearDepth( 1.0_gldouble )
CALL glClear(GL_DEPTH_BUFFER_BIT)
CALL glEnable(GL_DEPTH_TEST)
CALL glDepthFunc(GL_LEQUAL)
CALL argConvGlpara(patt_trans, gl_para) ! /* load the camera transformation matrix */
CALL glMatrixMode(GL_MODELVIEW)
CALL glLoadMatrixd( gl_para )

lightZeroColor  = [ 0.0_glfloat, 0.0_glfloat, 1.0_glfloat, 0.1_glfloat ]
ambi            = [ 0.0_glfloat, 0.1_glfloat, 1.0_glfloat, 0.1_glfloat ]
CALL glLightfv(GL_LIGHT0, GL_POSITION, [xl, yl, zl, wl]) ! light_position + light kind
CALL glLightfv(GL_LIGHT0, GL_AMBIENT , ambi)
CALL glLightfv(GL_LIGHT0, GL_DIFFUSE , lightZeroColor)

lightZeroColor = [ 1.0_glfloat, 0.0_glfloat, 0.0_glfloat, 0.1_glfloat ]
ambi           = [ 1.0_glfloat, 0.0_glfloat, 0.0_glfloat, 0.1_glfloat ]
CALL glLightfv(GL_LIGHT2, GL_POSITION, [-xl, yl, zl, wl]) ! light_position + light kind
CALL glLightfv(GL_LIGHT2, GL_AMBIENT , ambi)
CALL glLightfv(GL_LIGHT2, GL_DIFFUSE , lightZeroColor)

CALL glMaterialfv(GL_FRONT, GL_SPECULAR , mat_flash)
CALL glMaterialfv(GL_FRONT, GL_SHININESS, mat_flash_shiny)	
CALL glMaterialfv(GL_FRONT, GL_AMBIENT  , mat_ambient)
CALL glMatrixMode(GL_MODELVIEW)
!
IF (qfirst) THEN
 qfirst = .FALSE.
 CALL RANDOM_SEED()
 icount = 0
END IF
icount = icount + 1 
IF (icount == max_pos) stop
!
z0 = -200.0_glfloat
CALL RANDOM_NUMBER(x0)
CALL RANDOM_NUMBER(y0)
x0 = ( x0 - 0.5_glfloat ) * 0.0_glfloat
y0 = ( y0 - 0.5_glfloat ) * 20.0_glfloat
CALL set_time_interval(0.025)
CALL set_Coulomb_parameter(1.0)
p = init_particle( t_vec(x0, y0, z0), t_vec(0.0, 0.0, 1.0) )
DO i = 1, max_traject
 CALL step(p)
 traject(i, 1, icount, 1) = p%pos%x * 0.3_glfloat
 traject(i, 2, icount, 1) = p%pos%y * 0.3_glfloat
 traject(i, 3, icount, 1) = p%pos%z * 0.3_glfloat
END DO
!
CALL RANDOM_NUMBER(x0)
CALL RANDOM_NUMBER(y0)
x0 = ( x0 - 0.5_glfloat ) * 0.0_glfloat
y0 = ( y0 - 0.5_glfloat ) * 20.0_glfloat
CALL set_time_interval(0.025)
CALL set_Coulomb_parameter(-1.0)
p = init_particle( t_vec(x0, y0, z0), t_vec(0.0, 0.0, 1.0) )
DO i = 1, max_traject
 CALL step(p)
 traject(i, 1, icount, 2) = p%pos%x * 0.3_glfloat
 traject(i, 2, icount, 2) = p%pos%y * 0.3_glfloat
 traject(i, 3, icount, 2) = p%pos%z * 0.3_glfloat
END DO

!
CALL glTranslatef( 0.0_glfloat, 0.0_glfloat, 100.0 -z0 * 0.2_glfloat )
CALL glScalef(0.5_glfloat, 0.5_glfloat, 0.5_glfloat)
!
! light   position = [xl, yl, zl, wl]
!
CALL glPushMatrix()
CALL glPushAttrib(GL_LIGHTING_BIT)
CALL glColor3f(0.0_glfloat, 1.0_glfloat, 0.0_glfloat)
CALL glTranslatef( 0.0, 0.0, 0.0 )
CALL glutSolidSphere(5.0_gldouble, 10, 10)
CALL glPopAttrib()
CALL glPopMatrix()
!
! objects
!
CALL glPushMatrix()
!CALL glEnable(GL_LIGHTING)
!CALL glEnable(GL_LIGHT0)
!CALL glEnable(GL_LIGHT2)
!
CALL glLineWidth(1.0_glfloat)
DO ipos = 1, icount 
 CALL glPushAttrib(GL_LIGHTING_BIT)
 CALL glColor3f(0.0_glfloat, 0.0_glfloat, 1.0_glfloat)
 CALL glBegin(GL_LINE_STRIP)
 DO i = 1, max_traject
  CALL glVertex3f(-0.5 + traject(i, 1, ipos, 1), traject(i, 2, ipos, 1), traject(i, 3, ipos, 1) )
 END DO
 CALL glEnd()
 CALL glFlush()
 CALL glPopAttrib()
 
 CALL glPushAttrib(GL_LIGHTING_BIT)
 CALL glColor3f(1.0_glfloat, 0.0_glfloat, 0.0_glfloat)
 CALL glBegin(GL_LINE_STRIP)
 DO i = 1, max_traject
   CALL glVertex3f(0.5 + traject(i, 1, ipos, 2), traject(i, 2, ipos, 2), traject(i, 3, ipos, 2) )
 END DO
 CALL glEnd()
 CALL glFlush()
 CALL glPopAttrib()
 
 END DO 
CALL glPopMatrix()
!
CALL glDisable( GL_LIGHTING )
CALL glDisable( GL_DEPTH_TEST )
RETURN
END SUBROUTINE draw
!=====================================================
PROGRAM main1
USE m_art
IMPLICIT NONE
CALL glutInit()
CALL init()
CALL arVideoCapStart()
CALL argMainLoop( 0, C_FUNLOC(keyEvent), C_FUNLOC(mainLoop) )
STOP
END PROGRAM main1

Rutherford scattering

原子核の存在を示唆することになった、ラザフォードによるα粒子散乱実験の、モデル計算を試みます。(ゆくゆくは散乱ポテンシャルを変えて散乱断面積の理論値を散乱体の影を見て実感したいという、おぼろげな希望が背景にあります。)

固定されたクーロン型のポテンシャルにマーカー裏側から平行に荷電粒子が入射してきて、斥力型の力を受けて散乱されます。ラザフォードの実験では金薄膜にα線を入射したところ、ほとんどはスカスカと通り抜けるにもかかわらず、まれに大きな角度に曲げられて散乱する場合があり、さらには入射側にまで跳ね返るように戻ってくる場合もあることを見出し、これから原子核電荷が空間的にごく小さい領域に密集していることを示唆しました。

実行結果

赤い球体の位置に散乱体を置いています。入射側にはじき返されている粒子の軌跡が見えます。

入射粒子の出発座標は乱数で決めています。矩形の領域から出発します。

入射側と逆側から見てみました。ポテンシャルの対称性から散乱体を中心に円形に散乱されてゆきます。

スケールを変えて、全体像を見ています。

散乱体を含むある平面を切り出したものです。対称性から本来はこれだけで必要な情報を得るのに十分です。散乱された粒子の軌跡の包絡線は二次曲線を描いているらしいです。

横から見ると散乱が平面内で起きていることが確かめられます。

ラザフォード散乱は散乱ポテンシャルの符号によらず(斥力か引力かによらないで)に同じ散乱の角度依存性を示すはずなので、それを確かめてみるのも一興かと思います。

ここでは計算を簡単にするため、座標系をスケールしているのですが、その単位がよくわからないという手抜き状態。本当は赤玉の大きさを散乱断面積の理論値にしたかったのですが・・・w

Fortran Source

ちょっとオブジェクト指向風に散乱ルーチンを書いてみました。でも運動方程式積分オイラー法と言う手抜き。暫定版ということで清書とかしてません。

直してませんが、配列の添え字の順番に混乱がありました。パラメータによっては、配列はみだしが起きます。

!=====================================================
MODULE m_scatter
IMPLICIT NONE
REAL :: dt = 0.1
!
TYPE :: t_vec
 REAL :: x, y, z
END TYPE
!
TYPE :: t_particle
 TYPE (t_vec) :: pos, vel
END TYPE
!
INTERFACE OPERATOR (+)
 MODULE PROCEDURE vec_plus
END INTERFACE
!
INTERFACE OPERATOR (*)
 MODULE PROCEDURE vec_a_mult
END INTERFACE
!
CONTAINS
!---------------------------------
SUBROUTINE set_time_interval(t)
REAL, INTENT(IN) :: t
dt = t
RETURN
END SUBROUTINE set_time_interval
!---------------------------------
TYPE (t_vec) FUNCTION vec_plus(v1, v2)
TYPE (t_vec), INTENT(IN) :: v1, v2
vec_plus%x = v1%x + v2%x
vec_plus%y = v1%y + v2%y
vec_plus%z = v1%z + v2%z
RETURN
END FUNCTION vec_plus
!---------------------------------
TYPE (t_vec) FUNCTION vec_a_mult(a, v)
REAL, INTENT(IN) :: a
TYPE (t_vec), INTENT(IN) :: v
vec_a_mult%x = a * v%x
vec_a_mult%y = a * v%y
vec_a_mult%z = a * v%z
RETURN
END FUNCTION vec_a_mult
!---------------------------------
TYPE (t_particle) FUNCTION init_particle(pos_, vel_)
TYPE (t_vec), INTENT(IN) :: pos_, vel_
init_particle = t_particle( pos_, vel_ )
RETURN
END FUNCTION init_particle
!---------------------------------
REAL FUNCTION radius(v)
TYPE (t_vec), INTENT(IN) :: v
radius = SQRT(v%x ** 2 + v%y **2 + v%z ** 2)
RETURN
END FUNCTION radius
!---------------------------------
TYPE (t_vec) FUNCTION force(pos)
TYPE (t_vec), INTENT(IN) :: pos
REAL, PARAMETER :: alpha_prime = 1.0
REAL :: r
r = radius(pos)
force%x = alpha_prime * pos%x / r ** 3 ! derivative of the Coulomb potential
force%y = alpha_prime * pos%y / r ** 3
force%z = alpha_prime * pos%z / r ** 3
RETURN
END FUNCTION force
!---------------------------------
SUBROUTINE step(p) ! Euler's method
TYPE (t_particle), INTENT(IN OUT) :: p
TYPE (t_vec) :: dv, dx
p%vel = p%vel + dt * force(p%pos)
p%pos = p%pos + dt * p%vel
RETURN
END SUBROUTINE step
!---------------------------------
END MODULE m_scatter
!=====================================================
MODULE m_ART
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
!---- TYPE definitions ---------------------------------
TYPE, BIND(C) :: t_ARParam
 INTEGER(C_INT) :: ixsize, iysize
 REAL(C_DOUBLE) :: dmat(4, 3)
 REAL(C_DOUBLE) :: dist_factor(4)
END TYPE
!
TYPE, BIND(C) :: t_ARMarkerInfo 
 INTEGER(C_INT) :: iarea
 INTEGER(C_INT) :: id
 INTEGER(C_INT) :: dir
 REAL(C_DOUBLE)  :: cf
 REAL(C_DOUBLE)  :: pos(2)
 REAL(C_DOUBLE)  :: dline(3, 4)
 REAL(C_DOUBLE)  :: vertex(2, 4)
END TYPE
!---- global variables ---------------------------------
INTEGER(C_INT) :: ithresh = 100
INTEGER(C_INT) :: icount = 0
INTEGER(C_INT) :: ipatt_id
REAL(C_DOUBLE) :: patt_width = 80.0_c_double
REAL(C_DOUBLE) :: patt_center(2) = (/0.0_c_double, 0.0_c_double/)
REAL(C_DOUBLE) :: patt_trans(3, 4)
INTEGER(C_INT) :: ixsize, iysize
TYPE(t_ARParam) :: cparam
!---- INTERFACE BLOCK -----------------------------------
INTERFACE
 SUBROUTINE arVideoCapStop() BIND(C, NAME = 'arVideoCapStop')
 END SUBROUTINE arVideoCapStop
! 
 SUBROUTINE arVideoClose() BIND(C, NAME = 'arVideoClose')
 END SUBROUTINE arVideoClose
!
 SUBROUTINE argCleanup() BIND(C, NAME = 'argCleanup')
 END SUBROUTINE argCleanup
! 
 INTEGER(C_INT) FUNCTION arVideoGetImage() BIND(C, NAME = 'arVideoGetImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE arUtilSleep(i) BIND(C, NAME = 'arUtilSleep')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: i 
 END
!
 SUBROUTINE arUtilTimerReset() BIND(C, NAME = 'arUtilTimerReset')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE argDrawMode2D() BIND(C, NAME = 'argDrawMode2D')
 END 
!
 SUBROUTINE argDispImage( iptr, i, j ) BIND(C, NAME = 'argDispImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER, VALUE :: i, j
 END 
!
 INTEGER(C_INT) FUNCTION arDetectMarker( iptr, i, mark, k ) BIND(C, NAME = 'arDetectMarker')
 USE, INTRINSIC :: ISO_C_BINDING 
 IMPORT :: t_ARMarkerInfo
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER(C_INT)     , VALUE :: i
 TYPE(t_ARMarkerInfo), POINTER :: mark(:)
 INTEGER(C_INT) :: k
 END 
!
 SUBROUTINE arVideoCapNext() BIND(C, NAME = 'arVideoCapNext')
 END 
!
 SUBROUTINE argSwapBuffers() BIND(C, NAME = 'argSwapBuffers')
 END 
!
 SUBROUTINE arGetTransMat(ptr, d1, d, d2) BIND(C, NAME = 'arGetTransMat')
 USE, INTRINSIC :: ISO_C_BINDING 
 TYPE(C_PTR), VALUE :: ptr
 REAL(C_DOUBLE), VALUE :: d
 REAL(C_DOUBLE):: d1(2), d2(3, 4)
 END 
!
 SUBROUTINE argDrawMode3D() BIND(C, NAME = 'argDrawMode3D')
 END SUBROUTINE argDrawMode3D
 !
 SUBROUTINE argDraw3dCamera(ix, iy) BIND(C, NAME = 'argDraw3dCamera')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: ix, iy 
 END
 !
 SUBROUTINE argConvGlpara(x, y) BIND(C, NAME = 'argConvGlpara')
 USE, INTRINSIC :: ISO_C_BINDING 
 REAL(C_DOUBLE) :: x(4, 3), y(4, 4) 
 END
! 
 SUBROUTINE arVideoCapStart() BIND(C, NAME = 'arVideoCapStart')
 END
!
 SUBROUTINE argMainLoop( i, j, k ) BIND(C, NAME = 'argMainLoop')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER, VALUE :: i
 TYPE(C_FUNPTR), VALUE :: j, k
 END
!
 INTEGER(C_INT) FUNCTION arVideoOpen(text) BIND(C, NAME = 'arVideoOpen')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(C_CHAR) :: text(*)
 END
!
 INTEGER(C_INT) FUNCTION arVideoInqSize(ix, iy) BIND(C, NAME = 'arVideoInqSize')
 USE, INTRINSIC :: ISO_C_BINDING
 INTEGER(C_INT) :: ix, iy
 END
! 
 INTEGER(C_INT) FUNCTION arParamLoad(text, i, xparam) BIND(C, NAME = 'arParamLoad')
 USE, INTRINSIC :: ISO_C_BINDING
 IMPORT :: t_ARParam
 CHARACTER(C_CHAR) :: text(*)
 TYPE(t_ARParam) :: xparam
 INTEGER(C_INT), VALUE :: i
 END
!
 SUBROUTINE arParamChangeSize(xparam, i, j, yparam) BIND(C, NAME = 'arParamChangeSize')
 USE, INTRINSIC :: ISO_C_BINDING
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam, yparam
 INTEGER(C_INT), VALUE :: i, j
 END
! 
 SUBROUTINE arInitCparam(xparam) BIND(C, NAME = 'arInitCparam' )
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
! 
 SUBROUTINE arParamDisp(xparam) BIND(C, NAME = 'arParamDisp')
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
!
 INTEGER(C_INT) FUNCTION arLoadPatt(text) BIND(C, NAME = 'arLoadPatt')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(C_CHAR) :: text(*)
 END
!
 SUBROUTINE argInit(iptr, d, i1, i2, i3, i4) BIND(C, NAME = 'argInit')
 USE, INTRINSIC :: ISO_C_BINDING
 TYPE(C_PTR), VALUE :: iptr
 REAL(C_DOUBLE)      , VALUE :: d
 INTEGER(C_INT)     , VALUE :: i1, i2, i3, i4
 END
!
 REAL(C_DOUBLE) FUNCTION arUtilTimer() BIND(C, NAME = 'arUtilTimer')
 USE, INTRINSIC :: ISO_C_BINDING
 END FUNCTION arUtilTimer
END INTERFACE
CONTAINS
!-------------------------------------------------
SUBROUTINE init() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
TYPE(t_ARParam) :: wparam
CHARACTER (LEN = 35) :: vconf = "Data/WDM_camera_flipV.xml" // c_null_char
CHARACTER (LEN = 35) :: cparam_name = "Data/camera_para.dat" // c_null_char
CHARACTER (LEN = 35) :: patt_name   = "Data/patt.hiro"       // c_null_char
!   /* open the video path */
IF( arVideoOpen( vconf ) < 0 ) STOP
!    /* find the size of the window */
IF( arVideoInqSize( ixsize, iysize ) < 0 ) STOP
PRINT *, "Image size (x,y) = ", ixsize, iysize
!    /* set the initial camera parameters */
IF ( arParamLoad( cparam_name, 1, wparam ) < 0 ) THEN
 PRINT *, "Camera parameter load error !!"
 STOP
END IF
CALL arParamChangeSize( wparam, ixsize, iysize, cparam )
CALL arInitCparam( cparam )
PRINT *, "*** Camera Parameter ***"
CALL arParamDisp( cparam )
ipatt_id = arLoadPatt( patt_name )
IF ( ipatt_id < 0 ) THEN
 PRINT *, "pattern load error !!"
 STOP
END IF
!    /* open the graphics window */
CALL argInit( C_LOC(cparam), 1.0_c_double, 0, 0, 0, 0 )
RETURN
END SUBROUTINE init
!-------------------------------------------------
SUBROUTINE mainLoop() BIND(C, NAME = 'mainLoop')
IMPLICIT NONE
INTEGER(C_INTPTR_T) :: dataPtr
TYPE(t_ARMarkerInfo), POINTER:: marker_info(:)
INTEGER(C_INT) :: marker_num
INTEGER(C_INT) :: j, k
LOGICAL, SAVE :: qfirst = .TRUE.
IF (qfirst) THEN 
 ALLOCATE(marker_info(50)) ! ?
 qfirst = .FALSE.
END IF
dataPtr = arVideoGetImage()
IF (dataPtr == 0) THEN
 CALL arUtilSleep(2)
 RETURN
END IF
IF (icount == 0) CALL arUtilTimerReset()
icount = icount + 1
CALL argDrawMode2D()
CALL argDispImage( dataPtr, 0, 0 )
IF ( arDetectMarker( dataPtr, ithresh, marker_info, marker_num ) < 0 ) THEN
 CALL cleanup()
 RETURN
END IF
CALL arVideoCapNext()
k = -1
DO j = 1, marker_num
 IF (ipatt_id == marker_info(j)%id) THEN
  IF (k == -1) THEN 
   k = j
  ELSE 
   IF (marker_info(k)%cf < marker_info(j)%cf ) k = j
  END IF
 END IF
END DO
IF ( k == -1 ) THEN 
 CALL argSwapBuffers()
 RETURN
END IF
! draw-main 
CALL arGetTransMat( C_LOC(marker_info(k)), patt_center, patt_width, patt_trans )
CALL draw()
CALL argSwapBuffers()
RETURN
END SUBROUTINE mainLoop
!-------------------------------------------------
SUBROUTINE cleanup() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE 
CALL arVideoCapStop()
CALL arVideoClose()
CALL argCleanup()
RETURN
END SUBROUTINE cleanup
!-------------------------------------------------
SUBROUTINE keyEvent( key, ix, iy ) BIND(C, NAME = 'keyEvent')
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
INTEGER(C_SIGNED_CHAR), VALUE :: key
INTEGER(C_INT)        , VALUE :: ix, iy
IF (key == 27) THEN
  PRINT *, "*** ", icount / arUtilTimer(), "(frame/sec)"
  CALL cleanup()
  STOP
END IF
RETURN
END SUBROUTINE keyEvent
!---------------------------------------------------------------
SUBROUTINE draw() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
USE m_scatter
IMPLICIT NONE
REAL(gldouble) :: gl_para(4, 4)
REAL(glfloat ) :: mat_ambient(4) = [ 0.0_glfloat, 1.0_glfloat, 0.0_glfloat, 1.0_glfloat ]
REAL(glfloat ) :: mat_flash(4)   = [ 0.0_glfloat, 0.0_glfloat, 0.0_glfloat, 1.0_glfloat ]
REAL(glfloat ) :: mat_flash_shiny(1) = [ 50.0 ]
REAL(glfloat ) :: ambi(4)            = [ 0.1_glfloat, 0.1_glfloat, 0.1_glfloat, 0.1_glfloat ]
REAL(glfloat ) :: lightZeroColor(4)  = [ 0.0_glfloat, 1.0_glfloat, 1.0_glfloat, 0.1_glfloat ]
! shadow  
REAL(gldouble) :: shadowMatrix(4, 4), plane(4)
REAL(glfloat ), SAVE :: xl = 0.0_glfloat, yl = 0.0_glfloat, zl = -1.0_glfloat, wl = 0.0_glfloat, theta = 0.0
!
INTEGER, SAVE :: icount
LOGICAL, SAVE :: qfirst = .TRUE.
REAL :: x0, y0, z0
!
INTEGER, PARAMETER :: max_pos = 2000, max_traject = 2000
REAL(glfloat), SAVE :: traject(max_pos, 3, max_traject) 
TYPE (t_particle) :: p
INTEGER :: i, ipos
! 
CALL argDrawMode3D()
CALL argDraw3dCamera( 0, 0 )
CALL glClearDepth( 1.0_gldouble )
CALL glClear(GL_DEPTH_BUFFER_BIT)
CALL glEnable(GL_DEPTH_TEST)
CALL glDepthFunc(GL_LEQUAL)
CALL argConvGlpara(patt_trans, gl_para) ! /* load the camera transformation matrix */
CALL glMatrixMode(GL_MODELVIEW)
CALL glLoadMatrixd( gl_para )
CALL glLightfv(GL_LIGHT0, GL_POSITION, [xl, yl, zl, wl]) ! light_position + light kind
!CALL glLightfv(GL_LIGHT0, GL_AMBIENT , ambi)
CALL glLightfv(GL_LIGHT0, GL_DIFFUSE , lightZeroColor)
CALL glMaterialfv(GL_FRONT, GL_SPECULAR , mat_flash)
CALL glMaterialfv(GL_FRONT, GL_SHININESS, mat_flash_shiny)	
CALL glMaterialfv(GL_FRONT, GL_AMBIENT  , mat_ambient)
CALL glMatrixMode(GL_MODELVIEW)
!
IF (qfirst) THEN
 qfirst = .FALSE.
 CALL RANDOM_SEED()
 icount = 0
END IF
icount = icount + 1 
IF (icount == max_traject) stop
!
z0 = -20.0_glfloat
CALL RANDOM_NUMBER(x0)
CALL RANDOM_NUMBER(y0)
x0 = ( x0 - 0.5_glfloat ) * 20.0_glfloat
y0 = ( y0 - 0.5_glfloat ) * 20.0_glfloat
CALL set_time_interval(0.05)
p = init_particle( t_vec(x0, y0, z0), t_vec(0.0, 0.0, 1.0) )
DO i = 1, max_traject
 CALL step(p)
 traject(i, 1, icount) = p%pos%x * 5.0_glfloat
 traject(i, 2, icount) = p%pos%y * 5.0_glfloat
 traject(i, 3, icount) = p%pos%z * 5.0_glfloat
END DO

!
CALL glTranslatef( 0.0_glfloat, 0.0_glfloat, 10.0 -z0 * 5.0_glfloat )
CALL glScalef(0.5_glfloat, 0.5_glfloat, 0.5_glfloat)
!
! light   position = [xl, yl, zl, wl]
!
CALL glPushMatrix()
CALL glPushAttrib(GL_LIGHTING_BIT)
CALL glColor3f(1.0_glfloat, 0.0_glfloat, 0.1_glfloat)
CALL glTranslatef( 0.0, 0.0, 0.0 )
CALL glutSolidSphere(5.0_gldouble, 10, 10)
CALL glPopAttrib()
CALL glPopMatrix()
!
! objects
!
CALL glPushMatrix()
CALL glEnable(GL_LIGHTING)
CALL glEnable(GL_LIGHT0)
!
CALL glLineWidth(1.0_glfloat)
CALL glColor3f(1.0_glfloat, 1.0_glfloat, 1.0_glfloat)
DO i = 1, icount 
 CALL glBegin(GL_LINE_STRIP)
 DO ipos = 1, max_traject
   CALL glVertex3f( traject(ipos, 1, i), traject(ipos, 2, i), traject(ipos, 3, i) )
 END DO
 CALL glEnd()
 CALL glFlush()
END DO 
CALL glPopMatrix()
!
CALL glDisable( GL_LIGHTING )
CALL glDisable( GL_DEPTH_TEST )
RETURN
END SUBROUTINE draw
END MODULE m_ART
!=====================================================
PROGRAM main1
USE m_art
IMPLICIT NONE
CALL glutInit()
CALL init()
CALL arVideoCapStart()
CALL argMainLoop( 0, C_FUNLOC(keyEvent), C_FUNLOC(mainLoop) )
STOP
END PROGRAM main1

三次元ライフゲーム

コンウェイのライフゲームというものがあります。昔、マイコン上でよく遊ばれていたものです。あれは2次元でのものでしたが、OpenGL があるのでせっかくだから三次元でやってみることにしました。

結論から言うと、生命発生・維持・死亡のバランスが難しく、ちょっとパラメータをいじった範囲ではあまり面白いものになりませんでした。コロニー状の定常状態は発生しても、形を保って動いてゆく図形が現れないところがつまらないです。

ふと思いついて作り始めて、どうでもいい OpenGL での描画のデバッグに時間を食ってしまい、肝心のメインの部分を全然チェックしていません。


(むしろ、独立な二次元面を重ねて三次元にしたほうが面白いかもしれないと思いました。)

実行結果

40*40*40 の立方体状の単純格子上で実行しています。
光源を無限遠の平行光にしているので、赤い点は光源ではなく、マーカーの中心から見た光源の方向になっています。


Fortran ソース

  1. 境界の扱いですが周期境界条件ではなく、壁に囲まれたようになっています。
  2. 最初の分布は乱数で決めています。
  3. 生命誕生・存続・死亡の条件は最隣接サイトの占有数で決定しています。
    1. 占有数 1〜9 :淋しくて死亡。
    2. 占有数 10〜19 :生命誕生または生きながらえる。
    3. 占有数 20〜26 :人口過密で死亡。
MODULE m_ART
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
!---- TYPE definitions ---------------------------------
TYPE, BIND(C) :: t_ARParam
 INTEGER(C_INT) :: ixsize, iysize
 REAL(C_DOUBLE) :: dmat(4, 3)
 REAL(C_DOUBLE) :: dist_factor(4)
END TYPE
!
TYPE, BIND(C) :: t_ARMarkerInfo 
 INTEGER(C_INT) :: iarea
 INTEGER(C_INT) :: id
 INTEGER(C_INT) :: dir
 REAL(C_DOUBLE)  :: cf
 REAL(C_DOUBLE)  :: pos(2)
 REAL(C_DOUBLE)  :: dline(3, 4)
 REAL(C_DOUBLE)  :: vertex(2, 4)
END TYPE
!---- global variables ---------------------------------
INTEGER(C_INT) :: ithresh = 100
INTEGER(C_INT) :: icount = 0
INTEGER(C_INT) :: ipatt_id
REAL(C_DOUBLE) :: patt_width = 80.0_c_double
REAL(C_DOUBLE) :: patt_center(2) = (/0.0_c_double, 0.0_c_double/)
REAL(C_DOUBLE) :: patt_trans(3, 4)
INTEGER(C_INT) :: ixsize, iysize
TYPE(t_ARParam) :: cparam
!---- INTERFACE BLOCK -----------------------------------
INTERFACE
 SUBROUTINE arVideoCapStop() BIND(C, NAME = 'arVideoCapStop')
 END SUBROUTINE arVideoCapStop
! 
 SUBROUTINE arVideoClose() BIND(C, NAME = 'arVideoClose')
 END SUBROUTINE arVideoClose
!
 SUBROUTINE argCleanup() BIND(C, NAME = 'argCleanup')
 END SUBROUTINE argCleanup
! 
 INTEGER(C_INT) FUNCTION arVideoGetImage() BIND(C, NAME = 'arVideoGetImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE arUtilSleep(i) BIND(C, NAME = 'arUtilSleep')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: i 
 END
!
 SUBROUTINE arUtilTimerReset() BIND(C, NAME = 'arUtilTimerReset')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE argDrawMode2D() BIND(C, NAME = 'argDrawMode2D')
 END 
!
 SUBROUTINE argDispImage( iptr, i, j ) BIND(C, NAME = 'argDispImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER, VALUE :: i, j
 END 
!
 INTEGER(C_INT) FUNCTION arDetectMarker( iptr, i, mark, k ) BIND(C, NAME = 'arDetectMarker')
 USE, INTRINSIC :: ISO_C_BINDING 
 IMPORT :: t_ARMarkerInfo
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER(C_INT)     , VALUE :: i
 TYPE(t_ARMarkerInfo), POINTER :: mark(:)
 INTEGER(C_INT) :: k
 END 
!
 SUBROUTINE arVideoCapNext() BIND(C, NAME = 'arVideoCapNext')
 END 
!
 SUBROUTINE argSwapBuffers() BIND(C, NAME = 'argSwapBuffers')
 END 
!
 SUBROUTINE arGetTransMat(ptr, d1, d, d2) BIND(C, NAME = 'arGetTransMat')
 USE, INTRINSIC :: ISO_C_BINDING 
 TYPE(C_PTR), VALUE :: ptr
 REAL(C_DOUBLE), VALUE :: d
 REAL(C_DOUBLE):: d1(2), d2(3, 4)
 END 
!
 SUBROUTINE argDrawMode3D() BIND(C, NAME = 'argDrawMode3D')
 END SUBROUTINE argDrawMode3D
 !
 SUBROUTINE argDraw3dCamera(ix, iy) BIND(C, NAME = 'argDraw3dCamera')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: ix, iy 
 END
 !
 SUBROUTINE argConvGlpara(x, y) BIND(C, NAME = 'argConvGlpara')
 USE, INTRINSIC :: ISO_C_BINDING 
 REAL(C_DOUBLE) :: x(4, 3), y(4, 4) 
 END
! 
 SUBROUTINE arVideoCapStart() BIND(C, NAME = 'arVideoCapStart')
 END
!
 SUBROUTINE argMainLoop( i, j, k ) BIND(C, NAME = 'argMainLoop')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER, VALUE :: i
 TYPE(C_FUNPTR), VALUE :: j, k
 END
!
 INTEGER(C_INT) FUNCTION arVideoOpen(text) BIND(C, NAME = 'arVideoOpen')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(C_CHAR) :: text(*)
 END
!
 INTEGER(C_INT) FUNCTION arVideoInqSize(ix, iy) BIND(C, NAME = 'arVideoInqSize')
 USE, INTRINSIC :: ISO_C_BINDING
 INTEGER(C_INT) :: ix, iy
 END
! 
 INTEGER(C_INT) FUNCTION arParamLoad(text, i, xparam) BIND(C, NAME = 'arParamLoad')
 USE, INTRINSIC :: ISO_C_BINDING
 IMPORT :: t_ARParam
 CHARACTER(C_CHAR) :: text(*)
 TYPE(t_ARParam) :: xparam
 INTEGER(C_INT), VALUE :: i
 END
!
 SUBROUTINE arParamChangeSize(xparam, i, j, yparam) BIND(C, NAME = 'arParamChangeSize')
 USE, INTRINSIC :: ISO_C_BINDING
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam, yparam
 INTEGER(C_INT), VALUE :: i, j
 END
! 
 SUBROUTINE arInitCparam(xparam) BIND(C, NAME = 'arInitCparam' )
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
! 
 SUBROUTINE arParamDisp(xparam) BIND(C, NAME = 'arParamDisp')
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
!
 INTEGER(C_INT) FUNCTION arLoadPatt(text) BIND(C, NAME = 'arLoadPatt')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(C_CHAR) :: text(*)
 END
!
 SUBROUTINE argInit(iptr, d, i1, i2, i3, i4) BIND(C, NAME = 'argInit')
 USE, INTRINSIC :: ISO_C_BINDING
 TYPE(C_PTR), VALUE :: iptr
 REAL(C_DOUBLE)      , VALUE :: d
 INTEGER(C_INT)     , VALUE :: i1, i2, i3, i4
 END
!
 REAL(C_DOUBLE) FUNCTION arUtilTimer() BIND(C, NAME = 'arUtilTimer')
 USE, INTRINSIC :: ISO_C_BINDING
 END FUNCTION arUtilTimer
END INTERFACE
CONTAINS
!-------------------------------------------------
SUBROUTINE init() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
TYPE(t_ARParam) :: wparam
CHARACTER (LEN = 35) :: vconf = "Data/WDM_camera_flipV.xml" // c_null_char
CHARACTER (LEN = 35) :: cparam_name = "Data/camera_para.dat" // c_null_char
CHARACTER (LEN = 35) :: patt_name   = "Data/patt.hiro"       // c_null_char
!   /* open the video path */
IF( arVideoOpen( vconf ) < 0 ) STOP
!    /* find the size of the window */
IF( arVideoInqSize( ixsize, iysize ) < 0 ) STOP
PRINT *, "Image size (x,y) = ", ixsize, iysize
!    /* set the initial camera parameters */
IF ( arParamLoad( cparam_name, 1, wparam ) < 0 ) THEN
 PRINT *, "Camera parameter load error !!"
 STOP
END IF
CALL arParamChangeSize( wparam, ixsize, iysize, cparam )
CALL arInitCparam( cparam )
PRINT *, "*** Camera Parameter ***"
CALL arParamDisp( cparam )
ipatt_id = arLoadPatt( patt_name )
IF ( ipatt_id < 0 ) THEN
 PRINT *, "pattern load error !!"
 STOP
END IF
!    /* open the graphics window */
CALL argInit( C_LOC(cparam), 1.0_c_double, 0, 0, 0, 0 )
RETURN
END SUBROUTINE init
!-------------------------------------------------
SUBROUTINE mainLoop() BIND(C, NAME = 'mainLoop')
IMPLICIT NONE
INTEGER(C_INTPTR_T) :: dataPtr
TYPE(t_ARMarkerInfo), POINTER:: marker_info(:)
INTEGER(C_INT) :: marker_num
INTEGER(C_INT) :: j, k
LOGICAL, SAVE :: qfirst = .TRUE.
IF (qfirst) THEN 
 ALLOCATE(marker_info(50)) ! ?
 qfirst = .FALSE.
END IF
dataPtr = arVideoGetImage()
IF (dataPtr == 0) THEN
 CALL arUtilSleep(2)
 RETURN
END IF
IF (icount == 0) CALL arUtilTimerReset()
icount = icount + 1
CALL argDrawMode2D()
CALL argDispImage( dataPtr, 0, 0 )
IF ( arDetectMarker( dataPtr, ithresh, marker_info, marker_num ) < 0 ) THEN
 CALL cleanup()
 RETURN
END IF
CALL arVideoCapNext()
k = -1
DO j = 1, marker_num
 IF (ipatt_id == marker_info(j)%id) THEN
  IF (k == -1) THEN 
   k = j
  ELSE 
   IF (marker_info(k)%cf < marker_info(j)%cf ) k = j
  END IF
 END IF
END DO
IF ( k == -1 ) THEN 
 CALL argSwapBuffers()
 RETURN
END IF
! draw-main 
CALL arGetTransMat( C_LOC(marker_info(k)), patt_center, patt_width, patt_trans )
CALL draw()
CALL argSwapBuffers()
RETURN
END SUBROUTINE mainLoop
!-------------------------------------------------
SUBROUTINE cleanup() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE 
CALL arVideoCapStop()
CALL arVideoClose()
CALL argCleanup()
RETURN
END SUBROUTINE cleanup
!-------------------------------------------------
SUBROUTINE keyEvent( key, ix, iy ) BIND(C, NAME = 'keyEvent')
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
INTEGER(C_SIGNED_CHAR), VALUE :: key
INTEGER(C_INT)        , VALUE :: ix, iy
IF (key == 27) THEN
  PRINT *, "*** ", icount / arUtilTimer(), "(frame/sec)"
  CALL cleanup()
  STOP
END IF
RETURN
END SUBROUTINE keyEvent
!---------------------------------------------------------------
SUBROUTINE draw() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
REAL(gldouble) :: gl_para(4, 4)
REAL(glfloat ) :: mat_ambient(4) = [ 0.0_glfloat, 1.0_glfloat, 0.0_glfloat, 1.0_glfloat ]
REAL(glfloat ) :: mat_flash(4)   = [ 0.0_glfloat, 1.0_glfloat, 0.0_glfloat, 1.0_glfloat ]
REAL(glfloat ) :: mat_flash_shiny(1) = [ 50.0 ]
REAL(glfloat ) :: ambi(4)            = [ 0.1_glfloat, 0.1_glfloat, 0.1_glfloat, 0.1_glfloat ]
REAL(glfloat ) :: lightZeroColor(4)  = [ 0.0_glfloat, 1.0_glfloat, 1.0_glfloat, 0.1_glfloat ]
! shadow  
REAL(gldouble) :: shadowMatrix(4, 4), plane(4)
REAL(glfloat ), SAVE :: xl = 50.0_glfloat, yl = 50.0_glfloat, zl = 200.0_glfloat, wl = 0.0_glfloat, theta = 0.0
!
INTEGER, PARAMETER :: nz = 40
INTEGER, SAVE :: icount, mat(nz, nz, nz) = 0
INTEGER :: ix, iy, iz, m(nz, nz, nz)
LOGICAL, SAVE :: qfirst = .TRUE.
REAL, ALLOCATABLE :: tmp(:, :, :)
!   
! move light position
!
xl = 100.0 * COS(theta)
yl = 100.0 * SIN(theta)
theta = theta + 0.05
! 
CALL argDrawMode3D()
CALL argDraw3dCamera( 0, 0 )
CALL glClearDepth( 1.0_gldouble )
CALL glClear(GL_DEPTH_BUFFER_BIT)
CALL glEnable(GL_DEPTH_TEST)
CALL glDepthFunc(GL_LEQUAL)
CALL argConvGlpara(patt_trans, gl_para) ! /* load the camera transformation matrix */
CALL glMatrixMode(GL_MODELVIEW)
CALL glLoadMatrixd( gl_para )
CALL glLightfv(GL_LIGHT0, GL_POSITION, [xl, yl, zl, wl]) ! light_position + light kind
CALL glLightfv(GL_LIGHT0, GL_AMBIENT , ambi)
CALL glLightfv(GL_LIGHT0, GL_DIFFUSE , lightZeroColor)
CALL glMaterialfv(GL_FRONT, GL_SPECULAR , mat_flash)
CALL glMaterialfv(GL_FRONT, GL_SHININESS, mat_flash_shiny)	
CALL glMaterialfv(GL_FRONT, GL_AMBIENT  , mat_ambient)
CALL glMatrixMode(GL_MODELVIEW)
!
IF (qfirst) THEN
 qfirst = .FALSE.
 CALL RANDOM_SEED()
 ALLOCATE(tmp(nz, nz, nz))
 CALL RANDOM_NUMBER(tmp)
 WHERE (tmp > 0.74) mat = 1
 DEALLOCATE(tmp)
 icount = 0
END IF
icount = icount + 1 
!
! Life cycle
!
IF (MOD(icount, 10) == 0) THEN ! delay
 DO ix = 1, nz
  DO iy = 1, nz
   DO iz = 1, nz
    SELECT CASE( icount_nextn(mat, ix, iy, iz) )
     CASE (0:9,20:26)
      m(ix, iy, iz) = 0
     CASE (10:19)
      m(ix, iy, iz) = 1
     CASE DEFAULT
      STOP 'error impossible!'
     END SELECT 
   END DO
  END DO
 END DO 
 mat = m
END IF
!
! light   position = [xl, yl, zl, wl]
!
CALL glPushMatrix()
CALL glPushAttrib(GL_LIGHTING_BIT)
CALL glColor3f(1.0_glfloat, 0.0_glfloat, 0.1_glfloat)
CALL glTranslatef( xl, yl, zl )
CALL glutSolidSphere(5.0_gldouble, 10, 10)
CALL glPopAttrib()
CALL glPopMatrix()
!
! Shadow
!
plane = [ 0.0_gldouble,   0.0_gldouble,  1.0_gldouble,  0.0_gldouble]     ! plane normal; 4th element = intercept
CALL glPushMatrix()
CALL projectShadowMatrix(shadowMatrix, plane, [xl, yl, zl, wl]) ! light vector; 4th element = scaling factor
CALL glMultMatrixd(shadowMatrix)
!
CALL glScalef(5.0_glfloat, 5.0_glfloat, 5.0_glfloat)
CALL glPointSize(7.0_glfloat)
CALL glTranslatef( 0.0_glfloat, 0.0_glfloat, 0.0_glfloat )
CALL glColor3f(0.0_glfloat, 0.0_glfloat, 0.0_glfloat)
!
CALL glBegin(GL_POINTS)
DO ix = 1, nz
 DO iy = 1, nz
  DO iz = 1, nz
   IF (mat(ix, iy, iz) == 1) CALL glVertex3f(REAL(ix - nz / 2, glfloat), REAL(iy - nz / 2, glfloat), REAL(iz, glfloat))
  END DO
 END DO
END DO 
CALL glEnd()
CALL glPopMatrix()
!
! object
!
CALL glPushMatrix()
CALL glEnable(GL_LIGHTING)
CALL glEnable(GL_LIGHT0)
!
CALL glScalef(5.0_glfloat, 5.0_glfloat, 5.0_glfloat)
CALL glPointSize(7.0_glfloat)
CALL glTranslatef( 0.0_glfloat, 0.0_glfloat, 0.0_glfloat )
CALL glColor3f(1.0_glfloat, 1.0_glfloat, 1.0_glfloat)
CALL glBegin(GL_POINTS)
DO ix = 1, nz
 DO iy = 1, nz
  DO iz = 1, nz
   IF (mat(ix, iy, iz) == 1) CALL glVertex3f(REAL(ix - nz / 2, glfloat), REAL(iy - nz / 2, glfloat), REAL(iz, glfloat))
  END DO
 END DO
END DO 
CALL glEnd()
CALL glPopMatrix()
!
IF ( SUM(mat) == 0 ) qfirst = .TRUE. ! restart when all dead 
!
CALL glDisable( GL_LIGHTING )
CALL glDisable( GL_DEPTH_TEST )
RETURN
END SUBROUTINE draw
!----------------------------------------------------
INTEGER FUNCTION icount_nextn(mat, ix0, iy0, iz0) ! count next neighbour 
INTEGER, INTENT(IN) :: mat(:, :, :), ix0, iy0, iz0
INTEGER :: kx, ky, kz
icount_nextn = 0
DO kx = ix0 - 1, ix0 + 1
 IF ( kx < 1 .OR. kx > SIZE(mat, 1) ) CYCLE
 DO ky = iy0 - 1, iy0 + 1 
  IF ( ky < 1 .OR. ky > SIZE(mat, 2) ) CYCLE
  DO kz = iz0 - 1, iz0 + 1
   IF ( kz < 1 .OR. kz > SIZE(mat, 3) ) CYCLE
   IF ( kx == ix0 .AND. ky == iy0 .AND. kz == iz0 ) CYCLE ! skip center
   IF ( mat(kx, ky, kz) == 1 ) icount_nextn = icount_nextn + 1
  END DO
 END DO
END DO
RETURN
END FUNCTION icount_nextn
!-------------------------------------------------------------
SUBROUTINE projectShadowMatrix(sm, plane, light)
REAL(gldouble), INTENT(OUT) :: sm(:, :)
REAL(gldouble), INTENT(IN ) :: plane(:)
REAL(glfloat ), INTENT(IN ) :: light(:)
REAL(gldouble) :: dot
INTEGER :: i
sm = 0.0
dot = DOT_PRODUCT(plane, REAL(light, gldouble))     ! <plane|light> I  dot_product * diagonal Matrix
DO i = 1, 4
 sm(:, i) = -REAL(light(:), gldouble) * plane(i)    ! |light><plane|   4x4 Matrix
END DO
DO i = 1, 4
 sm(i, i) = dot + sm(i, i)                          
END DO
RETURN
END SUBROUTINE projectShadowMatrix
END MODULE m_ART
!=====================================================
PROGRAM main1
USE m_art
IMPLICIT NONE
CALL glutInit()
CALL init()
CALL arVideoCapStart()
CALL argMainLoop( 0, C_FUNLOC(keyEvent), C_FUNLOC(mainLoop) )
STOP
END PROGRAM main1

ARToolKit で影をつける。その2

以前の結果が間違っていたので訂正します。
http://d.hatena.ne.jp/fortran66/20090214

  • glutSolidCube を用いると、座標原点が動かされるのか、影を作る射影行列が正しく働かないようです。
    • そこで glVertex3f を用いて描画するようにしました。
  • 射影面の法線ベクトルの第 4 成分には、原点から射影平面への距離を負にしたものが入ります。
  • 光源の位置ベクトルの第 4 成分には、普通の座標点と同じくスケールファクターが入ります。ゆえにそれが 0 の時は、第 1 〜 3 成分 (x、y、z) 方向の無限遠に光源が存在することになり、それが 1 の時は、(x,y,z) 位置に点光源が存在することになるようです。
    • OpenGL の 4 元位置ベクトルでは第 4 成分がスケールファクターとなっていて、(x,y,z,w) = (x/w, y/w, z/w, 1.0) というスケール関係にあります。

影を生成する射影行列の導出

まず射影面の原点を通る場合を考えるとして、以下の三次元空間のベクトルを定義する。
|l> 光源の座標ベクトル。
|n> 原点から射影面までの距離ベクトル。(法線ベクトル)
|x> 影を生じさせる物体の座標ベクトル。
|s> 物体によって生じた影の座標ベクトル。

今求めたいものは、物体の座標を影の座標に変換する射影演算子。これをMとすると、
|s>=M|x>            (1)
と書ける。

影は射影面上にあるので、影のベクトル|s>は射影面の法線ベクトル|n>と直交する。ゆえに、それらの内積は消えて
[tex:=0]              (2)
となる。

影は光源と物体を結ぶ直線状にある。その直線はパラメータtを用いると、
|s>=t|x>+(1-t)|l>        (3)
とあらわせる。

式 (2),(3) より、
[tex:=t+(1-t)=0]     (4)
となる。

これを満たすパラメータtは、比例定数\alphaを導入して、
[tex:t=\alpha ,]          (5)
[tex:1-t=-\alpha ]         (6)
とあらわせる。

ここで式 (5),(6) から\alphaを求めると、
[tex:\alpha=(-)^{-1}]     (7)
となる。

以上の結果から影の座標は、
[tex:|s>=\alpha(|x>-|l>)]   (8)
と求まる。

これを形式的に|x>をくくり出した形で書き直すと
[tex:|s>=\alpha(I-|l>]    (9)
となる。ただし、I単位行列である。式 (1) と比較すると、射影の演算子 M
[tex:M=\alpha(I-|l> とあらわせる。

ここで、OpenGL の 4 元ベクトルを用いることを考えると、比例定数[tex:\alpha]を第 4 成分に引き受けさせることができる。[tex:|l>,|x>]の第 4 成分が 1 であることに注意すると、求めるべき第 4 成分も三次元ベクトルの時と同じ形式となることが分かるので、
[tex:M=(I-|l> と書きあらわせる。

また射影面が原点を通らない場合を考えた場合も、法線ベクトル[tex:|n>]の第 4 成分を、原点から平面までの距離に三次元法線ベクトルのノルムを掛けて負にしたものとすることで、式 (11) がそのまま成立することが導ける。

参照:
http://son-son.sakura.ne.jp/programming/opengl_1.html
http://www.opengl.org/resources/code/samples/mjktips/TexShadowReflectLight.html

□追記:式 (11) が原点を通らない平面でもそのまま成立すること

原点から射影平面へ直交する三次元のベクトルは、射影面のベクトルであるが、OpenGL の 4 元ベクトルにするときに、その第 4 成分を三次元ベクトルのノルムの二乗にしておく。これは原点から射影面までの距離の二乗に等しい。

射影面の法線ベクトルを成分にわけて書いたとき
|n>=|n_x, n_y, n_z, n_w>          (A.1)
であるとするとき、あらたに
|n_0>=|n_x, n_y, n_z, 0>,          (A.2)
を定義する。

ここで、定義によって
[tex:n_w=-]           (A.3)
の関係がある。

また座標原点を射影面に移したベクトルを'を付けることにして、
|l'>=|l>-|n_0>,          (A.4)
|x'>=|x>-|n_0>,          (A.5)
|s'>=|s>-|n_0>,          (A.6)
を定義する。

さて座標原点を射影面上に持ってきているので、影の生じる点は (8) より
[tex:|s'>=|x'>-|l'>]          (A.7)
とあらわせる。

今、
[tex:==-=-n_w-=], (A.8)  
[tex:==-=-n_w-=]  (A.9)
とまとめられる。

これを用いて、
[tex:|s'>=|x'>-|l'>=|x>-|l>-(-)|n_0>] (A.10)
と書き直せる。(ただし、右辺は全体を合わせた時のみ意味がある。)

ここで、[tex:-(-)|n_0>]の項を単独で切り離す場合は、OpenGL の4 元ベクトルの第 4 成分が規格化因子であることを思い出すと、[tex:|x>-|l>]の第 4 成分で割っておかなければならないことがわかる。|l>,|x>の第 4 成分は 1 であることに注意すると、係数は打ち消しあうことがわかり結局-|n_0>となる。

したがって、
[tex:|s>=|s'>+|n_0>=(|x>-|l>-|n_0>)+|n_0>=|x>-|l>]          (A.11)
と求める結果を得る。

解き方が第 4 成分の取り扱いに関して不自然ですっきりしません。OpenGL の第 4 成分を含めた代数がよく分かりません。普通の3次元ベクトルの形から始めて、だんだんと OpenGL の 4 元ベクトルに直して行ったので複雑になったかもしれません。不備なところはチラシの裏の走り書きということで勘弁してw

実行結果


Fortran のソース

光源(赤玉)が周回運動をするようにしてあります。

MODULE m_ART
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
!---- TYPE definitions ---------------------------------
TYPE, BIND(C) :: t_ARParam
 INTEGER(C_INT) :: ixsize, iysize
 REAL(C_DOUBLE) :: dmat(4, 3)
 REAL(C_DOUBLE) :: dist_factor(4)
END TYPE
!
TYPE, BIND(C) :: t_ARMarkerInfo 
 INTEGER(C_INT) :: iarea
 INTEGER(C_INT) :: id
 INTEGER(C_INT) :: dir
 REAL(C_DOUBLE)  :: cf
 REAL(C_DOUBLE)  :: pos(2)
 REAL(C_DOUBLE)  :: dline(3, 4)
 REAL(C_DOUBLE)  :: vertex(2, 4)
END TYPE
!---- global variables ---------------------------------
INTEGER(C_INT) :: ithresh = 100
INTEGER(C_INT) :: icount = 0
INTEGER(C_INT) :: ipatt_id
REAL(C_DOUBLE) :: patt_width = 80.0_c_double
REAL(C_DOUBLE) :: patt_center(2) = (/0.0_c_double, 0.0_c_double/)
REAL(C_DOUBLE) :: patt_trans(3, 4)
INTEGER(C_INT) :: ixsize, iysize
TYPE(t_ARParam) :: cparam
!---- INTERFACE BLOCK -----------------------------------
INTERFACE
 SUBROUTINE arVideoCapStop() BIND(C, NAME = 'arVideoCapStop')
 END SUBROUTINE arVideoCapStop
! 
 SUBROUTINE arVideoClose() BIND(C, NAME = 'arVideoClose')
 END SUBROUTINE arVideoClose
!
 SUBROUTINE argCleanup() BIND(C, NAME = 'argCleanup')
 END SUBROUTINE argCleanup
! 
 INTEGER(C_INT) FUNCTION arVideoGetImage() BIND(C, NAME = 'arVideoGetImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE arUtilSleep(i) BIND(C, NAME = 'arUtilSleep')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: i 
 END
!
 SUBROUTINE arUtilTimerReset() BIND(C, NAME = 'arUtilTimerReset')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE argDrawMode2D() BIND(C, NAME = 'argDrawMode2D')
 END 
!
 SUBROUTINE argDispImage( iptr, i, j ) BIND(C, NAME = 'argDispImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER, VALUE :: i, j
 END 
!
 INTEGER(C_INT) FUNCTION arDetectMarker( iptr, i, mark, k ) BIND(C, NAME = 'arDetectMarker')
 USE, INTRINSIC :: ISO_C_BINDING 
 IMPORT :: t_ARMarkerInfo
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER(C_INT)     , VALUE :: i
 TYPE(t_ARMarkerInfo), POINTER :: mark(:)
 INTEGER(C_INT) :: k
 END 
!
 SUBROUTINE arVideoCapNext() BIND(C, NAME = 'arVideoCapNext')
 END 
!
 SUBROUTINE argSwapBuffers() BIND(C, NAME = 'argSwapBuffers')
 END 
!
 SUBROUTINE arGetTransMat(ptr, d1, d, d2) BIND(C, NAME = 'arGetTransMat')
 USE, INTRINSIC :: ISO_C_BINDING 
 TYPE(C_PTR), VALUE :: ptr
 REAL(C_DOUBLE), VALUE :: d
 REAL(C_DOUBLE):: d1(2), d2(3, 4)
 END 
!
 SUBROUTINE argDrawMode3D() BIND(C, NAME = 'argDrawMode3D')
 END SUBROUTINE argDrawMode3D
 !
 SUBROUTINE argDraw3dCamera(ix, iy) BIND(C, NAME = 'argDraw3dCamera')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: ix, iy 
 END
 !
 SUBROUTINE argConvGlpara(x, y) BIND(C, NAME = 'argConvGlpara')
 USE, INTRINSIC :: ISO_C_BINDING 
 REAL(C_DOUBLE) :: x(4, 3), y(4, 4) 
 END
! 
 SUBROUTINE arVideoCapStart() BIND(C, NAME = 'arVideoCapStart')
 END
!
 SUBROUTINE argMainLoop( i, j, k ) BIND(C, NAME = 'argMainLoop')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER, VALUE :: i
 TYPE(C_FUNPTR), VALUE :: j, k
 END
!
 INTEGER(C_INT) FUNCTION arVideoOpen(text) BIND(C, NAME = 'arVideoOpen')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(C_CHAR) :: text(*)
 END
!
 INTEGER(C_INT) FUNCTION arVideoInqSize(ix, iy) BIND(C, NAME = 'arVideoInqSize')
 USE, INTRINSIC :: ISO_C_BINDING
 INTEGER(C_INT) :: ix, iy
 END
! 
 INTEGER(C_INT) FUNCTION arParamLoad(text, i, xparam) BIND(C, NAME = 'arParamLoad')
 USE, INTRINSIC :: ISO_C_BINDING
 IMPORT :: t_ARParam
 CHARACTER(C_CHAR) :: text(*)
 TYPE(t_ARParam) :: xparam
 INTEGER(C_INT), VALUE :: i
 END
!
 SUBROUTINE arParamChangeSize(xparam, i, j, yparam) BIND(C, NAME = 'arParamChangeSize')
 USE, INTRINSIC :: ISO_C_BINDING
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam, yparam
 INTEGER(C_INT), VALUE :: i, j
 END
! 
 SUBROUTINE arInitCparam(xparam) BIND(C, NAME = 'arInitCparam' )
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
! 
 SUBROUTINE arParamDisp(xparam) BIND(C, NAME = 'arParamDisp')
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
!
 INTEGER(C_INT) FUNCTION arLoadPatt(text) BIND(C, NAME = 'arLoadPatt')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(C_CHAR) :: text(*)
 END
!
 SUBROUTINE argInit(iptr, d, i1, i2, i3, i4) BIND(C, NAME = 'argInit')
 USE, INTRINSIC :: ISO_C_BINDING
 TYPE(C_PTR), VALUE :: iptr
 REAL(C_DOUBLE)      , VALUE :: d
 INTEGER(C_INT)     , VALUE :: i1, i2, i3, i4
 END
!
 REAL(C_DOUBLE) FUNCTION arUtilTimer() BIND(C, NAME = 'arUtilTimer')
 USE, INTRINSIC :: ISO_C_BINDING
 END FUNCTION arUtilTimer
END INTERFACE
CONTAINS
!-------------------------------------------------
SUBROUTINE init() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
TYPE(t_ARParam) :: wparam
CHARACTER (LEN = 35) :: vconf = "Data/WDM_camera_flipV.xml" // c_null_char
CHARACTER (LEN = 35) :: cparam_name = "Data/camera_para.dat" // c_null_char
CHARACTER (LEN = 35) :: patt_name   = "Data/patt.hiro"       // c_null_char
!   /* open the video path */
IF( arVideoOpen( vconf ) < 0 ) STOP
!    /* find the size of the window */
IF( arVideoInqSize( ixsize, iysize ) < 0 ) STOP
PRINT *, "Image size (x,y) = ", ixsize, iysize
!    /* set the initial camera parameters */
IF ( arParamLoad( cparam_name, 1, wparam ) < 0 ) THEN
 PRINT *, "Camera parameter load error !!"
 STOP
END IF
CALL arParamChangeSize( wparam, ixsize, iysize, cparam )
CALL arInitCparam( cparam )
PRINT *, "*** Camera Parameter ***"
CALL arParamDisp( cparam )
ipatt_id = arLoadPatt( patt_name )
IF ( ipatt_id < 0 ) THEN
 PRINT *, "pattern load error !!"
 STOP
END IF
!    /* open the graphics window */
CALL argInit( C_LOC(cparam), 1.0_c_double, 0, 0, 0, 0 )
RETURN
END SUBROUTINE init
!-------------------------------------------------
SUBROUTINE mainLoop() BIND(C, NAME = 'mainLoop')
IMPLICIT NONE
INTEGER(C_INTPTR_T) :: dataPtr
TYPE(t_ARMarkerInfo), POINTER:: marker_info(:)
INTEGER(C_INT) :: marker_num
INTEGER(C_INT) :: j, k
LOGICAL, SAVE :: qfirst = .TRUE.
IF (qfirst) THEN 
 ALLOCATE(marker_info(50)) ! ?
 qfirst = .FALSE.
END IF
dataPtr = arVideoGetImage()
IF (dataPtr == 0) THEN
 CALL arUtilSleep(2)
 RETURN
END IF
IF (icount == 0) CALL arUtilTimerReset()
icount = icount + 1
CALL argDrawMode2D()
CALL argDispImage( dataPtr, 0, 0 )
IF ( arDetectMarker( dataPtr, ithresh, marker_info, marker_num ) < 0 ) THEN
 CALL cleanup()
 RETURN
END IF
CALL arVideoCapNext()
k = -1
DO j = 1, marker_num
 IF (ipatt_id == marker_info(j)%id) THEN
  IF (k == -1) THEN 
   k = j
  ELSE 
   IF (marker_info(k)%cf < marker_info(j)%cf ) k = j
  END IF
 END IF
END DO
IF ( k == -1 ) THEN 
 CALL argSwapBuffers()
 RETURN
END IF
! draw-main 
CALL arGetTransMat( C_LOC(marker_info(k)), patt_center, patt_width, patt_trans )
CALL draw()
CALL argSwapBuffers()
RETURN
END SUBROUTINE mainLoop
!-------------------------------------------------
SUBROUTINE cleanup() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE 
CALL arVideoCapStop()
CALL arVideoClose()
CALL argCleanup()
RETURN
END SUBROUTINE cleanup
!-------------------------------------------------
SUBROUTINE keyEvent( key, ix, iy ) BIND(C, NAME = 'keyEvent')
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
INTEGER(C_SIGNED_CHAR), VALUE :: key
INTEGER(C_INT)        , VALUE :: ix, iy
IF (key == 27) THEN
  PRINT *, "*** ", icount / arUtilTimer(), "(frame/sec)"
  CALL cleanup()
  STOP
END IF
RETURN
END SUBROUTINE keyEvent
!---------------------------------------------------------------
SUBROUTINE draw() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
REAL(gldouble) :: gl_para(4, 4)
REAL(glfloat ) :: mat_ambient(4) = [ 0.0_glfloat, 0.0_glfloat, 1.0_glfloat, 1.0_glfloat ]
REAL(glfloat ) :: mat_flash(4)   = [ 0.0_glfloat, 0.0_glfloat, 1.0_glfloat, 1.0_glfloat ]
REAL(glfloat ) :: mat_flash_shiny(1) = [ 50.0 ]
REAL(glfloat ) :: ambi(4)            = [ 0.1_glfloat, 0.1_glfloat, 0.1_glfloat, 0.1_glfloat ]
REAL(glfloat ) :: lightZeroColor(4)  = [ 1.0_glfloat, 0.0_glfloat, 0.0_glfloat, 0.1_glfloat ]
REAL(glfloat ), PARAMETER :: a = 10.0_glfloat, b = 28.0_glfloat, &
                              c =  8.0_glfloat / 3.0_glfloat, d = 0.01_glfloat
! shadow  
REAL(gldouble) :: shadowMatrix(4, 4), plane(4)
REAL(glfloat ), SAVE :: xl = 50.0_glfloat, yl = 50.0_glfloat, zl = 200.0_glfloat, wl = 0.0_glfloat, theta = 0.0
!
! Lorenz
REAL(glfloat ) :: dx, dy, dz
REAL(glfloat ), SAVE :: x(3000), y(3000), z(3000)
INTEGER :: i
LOGICAL, SAVE :: first = .TRUE.


xl = 100.0 * COS(theta)
yl = 100.0 * SIN(theta)
theta = theta + 0.05
!
! Lorenz attractor
IF (first) THEN
 first = .FALSE.
 x(1) = 1.0_glfloat
 y(1) = 1.0_glfloat
 z(1) = 1.0_glfloat
 DO i = 1, 3000 - 1
  dx = d * ( a    * (y(i) - x(i))     )
  dy = d * ( x(i) * (b - z(i)) - y(i) )
  dz = d * ( x(i) * y(i) - c * z(i)   )
  x(i + 1) = x(i) + dx
  y(i + 1) = y(i) + dy
  z(i + 1) = z(i) + dz
 END DO
END IF
! 
CALL argDrawMode3D()
CALL argDraw3dCamera( 0, 0 )
CALL glClearDepth( 1.0_gldouble )
CALL glClear(GL_DEPTH_BUFFER_BIT)
CALL glEnable(GL_DEPTH_TEST)
CALL glDepthFunc(GL_LEQUAL)
CALL argConvGlpara(patt_trans, gl_para) ! /* load the camera transformation matrix */
CALL glMatrixMode(GL_MODELVIEW)
CALL glLoadMatrixd( gl_para )
CALL glLightfv(GL_LIGHT0, GL_POSITION, [xl, yl, zl, wl]) ! light_position + light kind
CALL glLightfv(GL_LIGHT0, GL_AMBIENT , ambi)
CALL glLightfv(GL_LIGHT0, GL_DIFFUSE , lightZeroColor)
CALL glMaterialfv(GL_FRONT, GL_SPECULAR , mat_flash)
CALL glMaterialfv(GL_FRONT, GL_SHININESS, mat_flash_shiny)	
CALL glMaterialfv(GL_FRONT, GL_AMBIENT  , mat_ambient)
CALL glMatrixMode(GL_MODELVIEW)
!
! light   position = [xl, yl, zl, wl]
!
CALL glPushMatrix()
CALL glPushAttrib(GL_LIGHTING_BIT)
CALL glColor3f(1.0_glfloat, 0.0_glfloat, 0.1_glfloat)
CALL glTranslatef( xl, yl, zl )
CALL glutSolidSphere(5.0_gldouble, 10, 10)
CALL glPopAttrib()
CALL glPopMatrix()
!
! Shadow
!
plane = [ 0.0_gldouble,   0.0_gldouble,  1.0_gldouble,  0.0_gldouble]     ! plane normal; 4th element = intercept
CALL glPushMatrix()
CALL projectShadowMatrix(shadowMatrix, plane, [xl, yl, zl, wl]) ! light vector; 4th element = scaling factor
CALL glMultMatrixd(shadowMatrix)
CALL glColor3f(0.1_glfloat, 0.1_glfloat, 0.1_glfloat)
!
CALL glLineWidth(2.0_glfloat)
CALL glBegin(GL_LINE_STRIP)
CALL glBegin(GL_POLYGON)
DO i = 1, 3000
 IF (i > 100) THEN
  CALL glVertex3f(2 * x(i), 2 * y(i), 2 * z(i))
 END IF 
END DO
CALL glEnd()
CALL glPopMatrix()
!
! Object
!
CALL glPushMatrix()
CALL glEnable(GL_LIGHTING)
CALL glEnable(GL_LIGHT0)
!
CALL glLineWidth(2.0_glfloat)
CALL glBegin(GL_LINE_STRIP)
CALL glBegin(GL_POLYGON)
DO i = 1, 3000
 IF (i > 100) THEN
  CALL glVertex3f(2 * x(i), 2 * y(i), 2 * z(i))
 END IF 
END DO
CALL glEnd()
CALL glPopMatrix()
!
CALL glDisable( GL_LIGHTING )
CALL glDisable( GL_DEPTH_TEST )
!
RETURN
END SUBROUTINE draw
!-------------------------------------------------------------
SUBROUTINE projectShadowMatrix(sm, plane, light)
REAL(gldouble), INTENT(OUT) :: sm(:, :)
REAL(gldouble), INTENT(IN ) :: plane(:)
REAL(glfloat ), INTENT(IN ) :: light(:)
REAL(gldouble) :: dot
INTEGER :: i
sm = 0.0
dot = DOT_PRODUCT(plane, REAL(light, gldouble))     ! <plane|light> I  dot_product * diagonal Matrix
DO i = 1, 4
 sm(:, i) = -REAL(light(:), gldouble) * plane(i)    ! |light><plane|   4x4 Matrix
END DO
DO i = 1, 4
 sm(i, i) = dot + sm(i, i)                          
END DO
RETURN
END SUBROUTINE projectShadowMatrix
!-------------------------------------------------------------
END MODULE m_ART
!=====================================================
PROGRAM shadow
USE m_art
IMPLICIT NONE
CALL glutInit()
CALL init()
CALL arVideoCapStart()
CALL argMainLoop( 0, C_FUNLOC(keyEvent), C_FUNLOC(mainLoop) )
STOP
END PROGRAM shadow

ARToolKit でローレンツ・アトラクターに影をつける。

OpenGL の範囲内で ARToolKit の像に影を付けることを目標とします。

適当にググったところ、以下のページに行き当たりました。簡潔な説明で非常に参考になります。単純な幾何学での導出がされていて理解しやすいです。実際のプログラム例もあって助かりました。(ただ、ベクトルの内積や行列とベクトルの積、テンソル積などが区別なく全部『*』の記号で表されているので読み分ける必要があります。)

ここでは影は幾何光学的な射影の範囲で求めており、回折などの波動光学的な扱いは考慮していません。
式の上では射影の有効演算子が、光源ベクトルと射影面の法線ベクトルの内積テンソル積の和の形に取れるのがキモのようです。

10行ぐらい付け足せばおkかと思ったのですが、今の場合はそのままではうまくいきませんでした。影の形は適切に変化しているようなのですが、影の位置が光源の移動によらずに原点にとどまってしまいました。私のいい加減な OpenGL の知識では、それが OpenGL の座標変換の問題なのか ARToolKit を使ったために暗黙に行われている処理の為なのか、よく分かりません。

ここでは半経験的に、もっともそうな座標移動操作を加えておきました。それっぽい位置に影を付けているような気もしますが、正しいのか分かりません。理解してません、すいませんw

CALL glTranslatef(-xl * ABS(xl / zl), -yl * ABS(yl / zl), 0.0)   ! ?

追記:2009-2-19 理解しましたw この日の結果は間違いです。詳しくは 2009-2-19 の日記で。

誤解の主たる理由

  1. glutSolidCube 等の glut 系の図形描画ルーチンは座標原点を暗黙に動かしているようで、glutSolidCube 等には、ここで用いた射影演算子をそのままでは適用できないようです。
  2. 影の生成式での第4成分(w成分)の意味を誤解していました。

  • Lorenz Attractor とその影
    • テスト版


  • Lorenz Attractor とその影および光源
    • 赤い点が光源の位置に対応しています。

光源の種類としては平行光線を与えています。

ローレンツ・アトラクターの裏表での光の当たり具合の差が以前より分かりにくくなりましたが、途中色々パラメータをいじりすぎたので原因が分かりませんw

以下のソースコードでは、C言語とのBINDINGのところと配列生成子記号([,])、IMPORT 命令に Fortran2003 の機能を用いています。
以下のプログラムは定式化が間違っているので、正しい結果を与えません。

MODULE m_ART
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
!---- TYPE definitions ---------------------------------
TYPE, BIND(C) :: t_ARParam
 INTEGER(C_INT) :: ixsize, iysize
 REAL(C_DOUBLE) :: dmat(4, 3)
 REAL(C_DOUBLE) :: dist_factor(4)
END TYPE
!
TYPE, BIND(C) :: t_ARMarkerInfo 
 INTEGER(C_INT) :: iarea
 INTEGER(C_INT) :: id
 INTEGER(C_INT) :: dir
 REAL(C_DOUBLE)  :: cf
 REAL(C_DOUBLE)  :: pos(2)
 REAL(C_DOUBLE)  :: dline(3, 4)
 REAL(C_DOUBLE)  :: vertex(2, 4)
END TYPE
!---- global variables ---------------------------------
INTEGER(C_INT) :: ithresh = 100
INTEGER(C_INT) :: icount = 0
INTEGER(C_INT) :: ipatt_id
REAL(C_DOUBLE) :: patt_width = 80.0_c_double
REAL(C_DOUBLE) :: patt_center(2) = (/0.0_c_double, 0.0_c_double/)
REAL(C_DOUBLE) :: patt_trans(3, 4)
INTEGER(C_INT) :: ixsize, iysize
TYPE(t_ARParam) :: cparam
!---- INTERFACE BLOCK -----------------------------------
INTERFACE
 SUBROUTINE arVideoCapStop() BIND(C, NAME = 'arVideoCapStop')
 END SUBROUTINE arVideoCapStop
! 
 SUBROUTINE arVideoClose() BIND(C, NAME = 'arVideoClose')
 END SUBROUTINE arVideoClose
!
 SUBROUTINE argCleanup() BIND(C, NAME = 'argCleanup')
 END SUBROUTINE argCleanup
! 
 INTEGER(C_INT) FUNCTION arVideoGetImage() BIND(C, NAME = 'arVideoGetImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE arUtilSleep(i) BIND(C, NAME = 'arUtilSleep')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: i 
 END
!
 SUBROUTINE arUtilTimerReset() BIND(C, NAME = 'arUtilTimerReset')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE argDrawMode2D() BIND(C, NAME = 'argDrawMode2D')
 END 
!
 SUBROUTINE argDispImage( iptr, i, j ) BIND(C, NAME = 'argDispImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER, VALUE :: i, j
 END 
!
 INTEGER(C_INT) FUNCTION arDetectMarker( iptr, i, mark, k ) BIND(C, NAME = 'arDetectMarker')
 USE, INTRINSIC :: ISO_C_BINDING 
 IMPORT :: t_ARMarkerInfo
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER(C_INT)     , VALUE :: i
 TYPE(t_ARMarkerInfo), POINTER :: mark(:)
 INTEGER(C_INT) :: k
 END 
!
 SUBROUTINE arVideoCapNext() BIND(C, NAME = 'arVideoCapNext')
 END 
!
 SUBROUTINE argSwapBuffers() BIND(C, NAME = 'argSwapBuffers')
 END 
!
 SUBROUTINE arGetTransMat(ptr, d1, d, d2) BIND(C, NAME = 'arGetTransMat')
 USE, INTRINSIC :: ISO_C_BINDING 
 TYPE(C_PTR), VALUE :: ptr
 REAL(C_DOUBLE), VALUE :: d
 REAL(C_DOUBLE):: d1(2), d2(3, 4)
 END 
!
 SUBROUTINE argDrawMode3D() BIND(C, NAME = 'argDrawMode3D')
 END SUBROUTINE argDrawMode3D
 !
 SUBROUTINE argDraw3dCamera(ix, iy) BIND(C, NAME = 'argDraw3dCamera')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: ix, iy 
 END
 !
 SUBROUTINE argConvGlpara(x, y) BIND(C, NAME = 'argConvGlpara')
 USE, INTRINSIC :: ISO_C_BINDING 
 REAL(C_DOUBLE) :: x(4, 3), y(4, 4) 
 END
! 
 SUBROUTINE arVideoCapStart() BIND(C, NAME = 'arVideoCapStart')
 END
!
 SUBROUTINE argMainLoop( i, j, k ) BIND(C, NAME = 'argMainLoop')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER, VALUE :: i
 TYPE(C_FUNPTR), VALUE :: j, k
 END
!
 INTEGER(C_INT) FUNCTION arVideoOpen(text) BIND(C, NAME = 'arVideoOpen')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(C_CHAR) :: text(*)
 END
!
 INTEGER(C_INT) FUNCTION arVideoInqSize(ix, iy) BIND(C, NAME = 'arVideoInqSize')
 USE, INTRINSIC :: ISO_C_BINDING
 INTEGER(C_INT) :: ix, iy
 END
! 
 INTEGER(C_INT) FUNCTION arParamLoad(text, i, xparam) BIND(C, NAME = 'arParamLoad')
 USE, INTRINSIC :: ISO_C_BINDING
 IMPORT :: t_ARParam
 CHARACTER(C_CHAR) :: text(*)
 TYPE(t_ARParam) :: xparam
 INTEGER(C_INT), VALUE :: i
 END
!
 SUBROUTINE arParamChangeSize(xparam, i, j, yparam) BIND(C, NAME = 'arParamChangeSize')
 USE, INTRINSIC :: ISO_C_BINDING
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam, yparam
 INTEGER(C_INT), VALUE :: i, j
 END
! 
 SUBROUTINE arInitCparam(xparam) BIND(C, NAME = 'arInitCparam' )
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
! 
 SUBROUTINE arParamDisp(xparam) BIND(C, NAME = 'arParamDisp')
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
!
 INTEGER(C_INT) FUNCTION arLoadPatt(text) BIND(C, NAME = 'arLoadPatt')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(C_CHAR) :: text(*)
 END
!
 SUBROUTINE argInit(iptr, d, i1, i2, i3, i4) BIND(C, NAME = 'argInit')
 USE, INTRINSIC :: ISO_C_BINDING
 TYPE(C_PTR), VALUE :: iptr
 REAL(C_DOUBLE)      , VALUE :: d
 INTEGER(C_INT)     , VALUE :: i1, i2, i3, i4
 END
!
 REAL(C_DOUBLE) FUNCTION arUtilTimer() BIND(C, NAME = 'arUtilTimer')
 USE, INTRINSIC :: ISO_C_BINDING
 END FUNCTION arUtilTimer
END INTERFACE
CONTAINS
!-------------------------------------------------
SUBROUTINE init() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
TYPE(t_ARParam) :: wparam
CHARACTER (LEN = 35) :: vconf = "Data/WDM_camera_flipV.xml" // c_null_char
CHARACTER (LEN = 35) :: cparam_name = "Data/camera_para.dat" // c_null_char
CHARACTER (LEN = 35) :: patt_name   = "Data/patt.hiro"       // c_null_char
!   /* open the video path */
IF( arVideoOpen( vconf ) < 0 ) STOP
!    /* find the size of the window */
IF( arVideoInqSize( ixsize, iysize ) < 0 ) STOP
PRINT *, "Image size (x,y) = ", ixsize, iysize
!    /* set the initial camera parameters */
IF ( arParamLoad( cparam_name, 1, wparam ) < 0 ) THEN
 PRINT *, "Camera parameter load error !!"
 STOP
END IF
CALL arParamChangeSize( wparam, ixsize, iysize, cparam )
CALL arInitCparam( cparam )
PRINT *, "*** Camera Parameter ***"
CALL arParamDisp( cparam )
ipatt_id = arLoadPatt( patt_name )
IF ( ipatt_id < 0 ) THEN
 PRINT *, "pattern load error !!"
 STOP
END IF
!    /* open the graphics window */
CALL argInit( C_LOC(cparam), 1.0_c_double, 0, 0, 0, 0 )
RETURN
END SUBROUTINE init
!-------------------------------------------------
SUBROUTINE mainLoop() BIND(C, NAME = 'mainLoop')
IMPLICIT NONE
INTEGER(C_INTPTR_T) :: dataPtr
TYPE(t_ARMarkerInfo), POINTER:: marker_info(:)
INTEGER(C_INT) :: marker_num
INTEGER(C_INT) :: j, k
LOGICAL, SAVE :: qfirst = .TRUE.
IF (qfirst) THEN 
 ALLOCATE(marker_info(50)) ! ?
 qfirst = .FALSE.
END IF
dataPtr = arVideoGetImage()
IF (dataPtr == 0) THEN
 CALL arUtilSleep(2)
 RETURN
END IF
IF (icount == 0) CALL arUtilTimerReset()
icount = icount + 1
CALL argDrawMode2D()
CALL argDispImage( dataPtr, 0, 0 )
IF ( arDetectMarker( dataPtr, ithresh, marker_info, marker_num ) < 0 ) THEN
 CALL cleanup()
 RETURN
END IF
CALL arVideoCapNext()
k = -1
DO j = 1, marker_num
 IF (ipatt_id == marker_info(j)%id) THEN
  IF (k == -1) THEN 
   k = j
  ELSE 
   IF (marker_info(k)%cf < marker_info(j)%cf ) k = j
  END IF
 END IF
END DO
IF ( k == -1 ) THEN 
 CALL argSwapBuffers()
 RETURN
END IF
! draw-main 
CALL arGetTransMat( C_LOC(marker_info(k)), patt_center, patt_width, patt_trans )
CALL draw()
CALL argSwapBuffers()
RETURN
END SUBROUTINE mainLoop
!-------------------------------------------------
SUBROUTINE cleanup() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE 
CALL arVideoCapStop()
CALL arVideoClose()
CALL argCleanup()
RETURN
END SUBROUTINE cleanup
!-------------------------------------------------
SUBROUTINE keyEvent( key, ix, iy ) BIND(C, NAME = 'keyEvent')
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
INTEGER(C_SIGNED_CHAR), VALUE :: key
INTEGER(C_INT)        , VALUE :: ix, iy
IF (key == 27) THEN
  PRINT *, "*** ", icount / arUtilTimer(), "(frame/sec)"
  CALL cleanup()
  STOP
END IF
RETURN
END SUBROUTINE keyEvent
!---------------------------------------------------------------
SUBROUTINE draw() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
REAL(gldouble) :: gl_para(4, 4)
REAL(glfloat ) :: mat_ambient(4) = [ 0.0_glfloat, 0.0_glfloat, 1.0_glfloat, 1.0_glfloat ]
REAL(glfloat ) :: mat_flash(4)   = [ 0.0_glfloat, 0.0_glfloat, 1.0_glfloat, 1.0_glfloat ]
REAL(glfloat ) :: mat_flash_shiny(1) = [ 50.0 ]
REAL(glfloat ) :: ambi(4)            = [ 0.1_glfloat, 0.1_glfloat, 0.1_glfloat, 0.1_glfloat ]
REAL(glfloat ) :: lightZeroColor(4)  = [ 1.0_glfloat, 0.0_glfloat, 0.0_glfloat, 0.1_glfloat ]
REAL(glfloat ), PARAMETER :: a = 10.0_glfloat, b = 28.0_glfloat, &
                              c =  8.0_glfloat / 3.0_glfloat, d = 0.01_glfloat
! shadow  
REAL(gldouble) :: shadowMatrix(4, 4), plane(4)
REAL(glfloat ), SAVE :: xl = 40.0_glfloat, yl = -40.0_glfloat, zl = 80.0_glfloat, wl = 0.0_glfloat
!
! Lorenz
REAL(glfloat ) :: x, y, z
REAL(glfloat ), SAVE :: dx(3000), dy(3000), dz(3000)
INTEGER :: i
LOGICAL, SAVE :: first = .TRUE.
!
! Lorenz attractor
IF (first) THEN
 first = .FALSE.
 x = 1.0_glfloat
 y = 1.0_glfloat
 z = 1.0_glfloat
 DO i = 1, 3000
  dx(i) = d * ( a * (y - x)     )
  dy(i) = d * ( x * (b - z) - y )
  dz(i) = d * ( x * y - c * z   )
  x = x + dx(i)
  y = y + dy(i)
  z = z + dz(i)
 END DO
END IF
! 
CALL argDrawMode3D()
CALL argDraw3dCamera( 0, 0 )
CALL glClearDepth( 1.0_gldouble )
CALL glClear(GL_DEPTH_BUFFER_BIT)
CALL glEnable(GL_DEPTH_TEST)
CALL glDepthFunc(GL_LEQUAL)
CALL argConvGlpara(patt_trans, gl_para) ! /* load the camera transformation matrix */
CALL glMatrixMode(GL_MODELVIEW)
CALL glLoadMatrixd( gl_para )
CALL glLightfv(GL_LIGHT0, GL_POSITION, [xl, yl, zl, wl]) ! light_position
CALL glLightfv(GL_LIGHT0, GL_AMBIENT , ambi)
CALL glLightfv(GL_LIGHT0, GL_DIFFUSE , lightZeroColor)
CALL glMaterialfv(GL_FRONT, GL_SPECULAR , mat_flash)
CALL glMaterialfv(GL_FRONT, GL_SHININESS, mat_flash_shiny)	
CALL glMaterialfv(GL_FRONT, GL_AMBIENT  , mat_ambient)
CALL glMatrixMode(GL_MODELVIEW)
!
! Scaling
!
CALL glScalef( 2.0, 2.0, 2.0)
!
! light   position = [xl, yl, zl, wl]
!
CALL glPushMatrix()
CALL glPushAttrib(GL_LIGHTING_BIT)
CALL glColor3f(1.0_glfloat, 0.0_glfloat, 0.1_glfloat)
CALL glTranslatef( xl, yl, zl )
CALL glutSolidSphere(2.0_gldouble, 10, 10)
CALL glPopAttrib()
CALL glPopMatrix()
!
! Shadow
!
plane = [ 0.0_gldouble,   0.0_gldouble,  1.0_gldouble,  0.0_gldouble] ! plane normal
CALL glPushMatrix()
CALL projectShadowMatrix(shadowMatrix, plane, [xl, yl, zl, wl])
CALL glMultMatrixd(shadowMatrix)
CALL glTranslatef(-xl * ABS(xl / zl), -yl * ABS(yl / zl), 0.0)   ! ?
!
CALL glBegin(GL_LINE_STRIP)
CALL glColor3f(0.1_glfloat, 0.1_glfloat, 0.1_glfloat)
DO i = 1, 3000
 CALL glTranslatef( dx(i), dy(i), dz(i) )
 IF (i > 100) THEN
  CALL glutSolidCube(1.0_gldouble)
 END IF 
END DO
CALL glPopMatrix()
!
! Object
!
CALL glEnable(GL_LIGHTING)
CALL glEnable(GL_LIGHT0)
CALL glPushMatrix()
DO i = 1, 3000
 CALL glTranslatef( dx(i), dy(i), dz(i) )
 IF (i > 100) THEN
  CALL glutSolidCube(1.0_gldouble )
 END IF 
END DO
CALL glPopMatrix()
CALL glDisable( GL_LIGHTING )
CALL glDisable( GL_DEPTH_TEST )
!
RETURN
END SUBROUTINE draw
!-------------------------------------------------------------
SUBROUTINE projectShadowMatrix(sm, plane, light)
REAL(gldouble), INTENT(OUT) :: sm(:, :)
REAL(gldouble), INTENT(IN ) :: plane(:)
REAL(glfloat ), INTENT(IN ) :: light(:)
REAL(gldouble) :: dot
INTEGER :: i
dot = DOT_PRODUCT(plane, REAL(light, gldouble))
DO i = 1, 4
 sm(:, i) = -REAL(light(:), gldouble) * plane(i)    ! |light><plane|   4x4 Tensor
END DO
DO i = 1, 4
 sm(i, i) = dot + sm(i, i)                          ! <plane|light> I  dot_product * diagonal Matrix
END DO
RETURN
END SUBROUTINE projectShadowMatrix
!-------------------------------------------------------------
END MODULE m_ART
!=====================================================
PROGRAM shadow
USE m_art
IMPLICIT NONE
CALL glutInit()
CALL init()
CALL arVideoCapStart()
CALL argMainLoop( 0, C_FUNLOC(keyEvent), C_FUNLOC(mainLoop) )
STOP
END PROGRAM shadow

[Fortran2003]ARToolKit で Lorenz attractor

しばらくぶりに、ARToolKitをいじってみます。
ローレンツ・アトラクターは複雑な三次元図形なので、ARToolKit で任意方向から眺め回せれば愉快かなとふと思う。去年作ったプログラムの描画部分を数行ばかり変えることで実現できるはず。

アトラクター部分は『C言語による 最新 アルゴリズム事典』を参照しました。

新たな環境にインストール 
必要条件 チラシ裏

  1. F90GL が動いていること。
  2. ARToolKit 展開。インストールはいらない。展開だけでおk。
    • 実行ファイルは展開ディレクトリ下の bin に置かないと、初期ファイルが読めない。
  3. LIB Path に、展開ディレクトリ下の lib を指定。
    • 必要ライブラリ glut32.lib f90gl.lib f90glu.lib f90glut.lib libAR.lib libARgsub.lib libARvideo.lib
  4. Fortran のコンソールアプリとして製作する。
  5. ESCで終了しないとゴミプロセスが溜まる。




インターフェースでのポインター定義などを微妙に修正。

MODULE m_ART
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
!---- TYPE definitions ---------------------------------
TYPE, BIND(C) :: t_ARParam
 INTEGER(C_INT) :: ixsize, iysize
 REAL(C_DOUBLE) :: dmat(4, 3)
 REAL(C_DOUBLE) :: dist_factor(4)
END TYPE
!
TYPE, BIND(C) :: t_ARMarkerInfo 
 INTEGER(C_INT) :: iarea
 INTEGER(C_INT) :: id
 INTEGER(C_INT) :: dir
 REAL(C_DOUBLE)  :: cf
 REAL(C_DOUBLE)  :: pos(2)
 REAL(C_DOUBLE)  :: dline(3, 4)
 REAL(C_DOUBLE)  :: vertex(2, 4)
END TYPE
!---- global variables ---------------------------------
INTEGER(C_INT) :: ithresh = 100
INTEGER(C_INT) :: icount = 0
INTEGER(C_INT) :: ipatt_id
REAL(C_DOUBLE) :: patt_width = 80.0_c_double
REAL(C_DOUBLE) :: patt_center(2) = (/0.0_c_double, 0.0_c_double/)
REAL(C_DOUBLE) :: patt_trans(3, 4)
INTEGER(C_INT) :: ixsize, iysize
TYPE(t_ARParam) :: cparam
!---- INTERFACE BLOCK -----------------------------------
INTERFACE
 SUBROUTINE arVideoCapStop() BIND(C, NAME = 'arVideoCapStop')
 END SUBROUTINE arVideoCapStop
! 
 SUBROUTINE arVideoClose() BIND(C, NAME = 'arVideoClose')
 END SUBROUTINE arVideoClose
!
 SUBROUTINE argCleanup() BIND(C, NAME = 'argCleanup')
 END SUBROUTINE argCleanup
! 
 INTEGER(C_INT) FUNCTION arVideoGetImage() BIND(C, NAME = 'arVideoGetImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE arUtilSleep(i) BIND(C, NAME = 'arUtilSleep')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: i 
 END
!
 SUBROUTINE arUtilTimerReset() BIND(C, NAME = 'arUtilTimerReset')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE argDrawMode2D() BIND(C, NAME = 'argDrawMode2D')
 END 
!
 SUBROUTINE argDispImage( iptr, i, j ) BIND(C, NAME = 'argDispImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER, VALUE :: i, j
 END 
!
 INTEGER(C_INT) FUNCTION arDetectMarker( iptr, i, mark, k ) BIND(C, NAME = 'arDetectMarker')
 USE, INTRINSIC :: ISO_C_BINDING 
 IMPORT :: t_ARMarkerInfo
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER(C_INT)     , VALUE :: i
 TYPE(t_ARMarkerInfo), POINTER :: mark(:)
 INTEGER(C_INT) :: k
 END 
!
 SUBROUTINE arVideoCapNext() BIND(C, NAME = 'arVideoCapNext')
 END 
!
 SUBROUTINE argSwapBuffers() BIND(C, NAME = 'argSwapBuffers')
 END 
!
 SUBROUTINE arGetTransMat(ptr, d1, d, d2) BIND(C, NAME = 'arGetTransMat')
 USE, INTRINSIC :: ISO_C_BINDING 
 TYPE(C_PTR), VALUE :: ptr
 REAL(C_DOUBLE), VALUE :: d
 REAL(C_DOUBLE):: d1(2), d2(3, 4)
 END 
!
 SUBROUTINE argDrawMode3D() BIND(C, NAME = 'argDrawMode3D')
 END SUBROUTINE argDrawMode3D
 !
 SUBROUTINE argDraw3dCamera(ix, iy) BIND(C, NAME = 'argDraw3dCamera')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: ix, iy 
 END
 !
 SUBROUTINE argConvGlpara(x, y) BIND(C, NAME = 'argConvGlpara')
 USE, INTRINSIC :: ISO_C_BINDING 
 REAL(C_DOUBLE) :: x(4, 3), y(4, 4) 
 END
! 
 SUBROUTINE arVideoCapStart() BIND(C, NAME = 'arVideoCapStart')
 END
!
 SUBROUTINE argMainLoop( i, j, k ) BIND(C, NAME = 'argMainLoop')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER, VALUE :: i
 TYPE(C_FUNPTR), VALUE :: j, k
 END
!
 INTEGER(4) FUNCTION arVideoOpen(text) BIND(C, NAME = 'arVideoOpen')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(C_CHAR) :: text(*)
 END
!
 INTEGER(C_INT) FUNCTION arVideoInqSize(ix, iy) BIND(C, NAME = 'arVideoInqSize')
 USE, INTRINSIC :: ISO_C_BINDING
 INTEGER(C_INT) :: ix, iy
 END
! 
 INTEGER(C_INT) FUNCTION arParamLoad(text, i, xparam) BIND(C, NAME = 'arParamLoad')
 USE, INTRINSIC :: ISO_C_BINDING
 IMPORT :: t_ARParam
 CHARACTER(C_CHAR) :: text(*)
 TYPE(t_ARParam) :: xparam
 INTEGER(C_INT), VALUE :: i
 END
!
 SUBROUTINE arParamChangeSize(xparam, i, j, yparam) BIND(C, NAME = 'arParamChangeSize')
 USE, INTRINSIC :: ISO_C_BINDING
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam, yparam
 INTEGER(C_INT), VALUE :: i, j
 END
! 
 SUBROUTINE arInitCparam(xparam) BIND(C, NAME = 'arInitCparam' )
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
! 
 SUBROUTINE arParamDisp(xparam) BIND(C, NAME = 'arParamDisp')
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
!
 INTEGER(4) FUNCTION arLoadPatt(text) BIND(C, NAME = 'arLoadPatt')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(C_CHAR) :: text(*)
 END
!
 SUBROUTINE argInit(iptr, d, i1, i2, i3, i4) BIND(C, NAME = 'argInit')
 USE, INTRINSIC :: ISO_C_BINDING
 TYPE(C_PTR), VALUE :: iptr
 REAL(C_DOUBLE)      , VALUE :: d
 INTEGER(C_INT)     , VALUE :: i1, i2, i3, i4
 END
!
 REAL(8) FUNCTION arUtilTimer() BIND(C, NAME = 'arUtilTimer')
 END FUNCTION arUtilTimer
END INTERFACE
CONTAINS
!-------------------------------------------------
SUBROUTINE init() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
TYPE(t_ARParam) :: wparam
CHARACTER (LEN = 35) :: vconf = "Data\\WDM_camera_flipV.xml" // c_null_char
CHARACTER (LEN = 35) :: cparam_name = "Data\\camera_para.dat" // c_null_char
CHARACTER (LEN = 35) :: patt_name   = "Data\\patt.hiro"       // c_null_char
!   /* open the video path */
IF( arVideoOpen( vconf ) < 0 ) STOP
!    /* find the size of the window */
IF( arVideoInqSize( ixsize, iysize ) < 0 ) STOP
PRINT *, "Image size (x,y) = ", ixsize, iysize
!    /* set the initial camera parameters */
IF ( arParamLoad( cparam_name, 1, wparam ) < 0 ) THEN
 PRINT *, "Camera parameter load error !!"
 STOP
END IF
CALL arParamChangeSize( wparam, ixsize, iysize, cparam )
CALL arInitCparam( cparam )
PRINT *, "*** Camera Parameter ***"
CALL arParamDisp( cparam )
ipatt_id = arLoadPatt( patt_name )
IF ( ipatt_id < 0 ) THEN
 PRINT *, "pattern load error !!"
 STOP
END IF
!    /* open the graphics window */
CALL argInit( C_LOC(cparam), 1.0_c_double, 0, 0, 0, 0 )
RETURN
END SUBROUTINE init
!-------------------------------------------------
SUBROUTINE mainLoop() BIND(C, NAME = 'mainLoop')
IMPLICIT NONE
INTEGER(C_INTPTR_T) :: dataPtr
TYPE(t_ARMarkerInfo), POINTER:: marker_info(:)
INTEGER(C_INT) :: marker_num
INTEGER(C_INT) :: j, k
LOGICAL, SAVE :: qfirst = .TRUE.
IF (qfirst) THEN 
 ALLOCATE(marker_info(50)) ! ?
 qfirst = .FALSE.
END IF
dataPtr = arVideoGetImage()
IF (dataPtr == 0) THEN
 CALL arUtilSleep(2)
 RETURN
END IF
IF (icount == 0) CALL arUtilTimerReset()
icount = icount + 1
CALL argDrawMode2D()
CALL argDispImage( dataPtr, 0, 0 )
IF ( arDetectMarker( dataPtr, ithresh, marker_info, marker_num ) < 0 ) THEN
 CALL cleanup()
 RETURN
END IF
CALL arVideoCapNext()
k = -1
DO j = 1, marker_num
 IF (ipatt_id == marker_info(j)%id) THEN
  IF (k == -1) THEN 
   k = j
  ELSE 
   IF (marker_info(k)%cf < marker_info(j)%cf ) k = j
  END IF
 END IF
END DO
IF ( k == -1 ) THEN 
 CALL argSwapBuffers()
 RETURN
END IF
! draw-main 
CALL arGetTransMat( C_LOC(marker_info(k)), patt_center, patt_width, patt_trans )
CALL draw()
CALL argSwapBuffers()
RETURN
END SUBROUTINE mainLoop
!-------------------------------------------------
SUBROUTINE cleanup() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE 
CALL arVideoCapStop()
CALL arVideoClose()
CALL argCleanup()
RETURN
END SUBROUTINE cleanup
!-------------------------------------------------
SUBROUTINE keyEvent( key, ix, iy ) BIND(C, NAME = 'keyEvent')
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
INTEGER(C_SIGNED_CHAR), VALUE :: key
INTEGER(C_INT)        , VALUE :: ix, iy
IF (key == 27) THEN
  PRINT *, "*** ", icount / arUtilTimer(), "(frame/sec)"
  CALL cleanup()
  STOP
END IF
RETURN
END SUBROUTINE keyEvent
!---------------------------------------------------------------
SUBROUTINE draw() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
REAL(gldouble) :: gl_para(4, 4)
REAL(glfloat ) :: mat_ambient(4) = (/0.0_glfloat, 0.0_glfloat, 1.0_glfloat, 1.0_glfloat/)
REAL(glfloat ) :: mat_flash(4)   = (/0.0_glfloat, 0.0_glfloat, 1.0_glfloat, 1.0_glfloat/)
REAL(glfloat ) :: mat_flash_shiny(1) = (/50.0/)
REAL(glfloat ) :: light_position(4)  = (/100.0_glfloat, -200.0_glfloat, 200.0_glfloat, 0.0_glfloat/)
REAL(glfloat ) :: ambi(4)            = (/0.1_glfloat, 0.1_glfloat, 0.1_glfloat, 0.1_glfloat/)
REAL(glfloat ) :: lightZeroColor(4)  = (/1.0_glfloat, 0.0_glfloat, 0.0_glfloat, 0.1_glfloat/)
REAL(glfloat ), PARAMETER :: a = 10.0_glfloat, b = 28.0_glfloat, &
                              c =  8.0_glfloat / 3.0_glfloat, d = 0.01_glfloat
REAL(glfloat ) :: x, y, z
REAL(glfloat ), SAVE :: dx(3000), dy(3000), dz(3000)
INTEGER :: i
LOGICAL, SAVE :: first = .TRUE.
!   
! Lorenz attractor
IF (first) THEN
 first = .FALSE.
 x = 1.0_glfloat
 y = 1.0_glfloat
 z = 1.0_glfloat
 DO i = 1, 3000
  dx(i) = a * (y - x)
  dy(i) = x * (b - z) - y
  dz(i) = x * y - c * z
  x = x + d * dx(i)
  y = y + d * dy(i)
  z = z + d * dz(i)
 END DO
END IF
! 
CALL argDrawMode3D()
CALL argDraw3dCamera( 0, 0 )
CALL glClearDepth( 1.0_gldouble )
CALL glClear(GL_DEPTH_BUFFER_BIT)
CALL glEnable(GL_DEPTH_TEST)
CALL glDepthFunc(GL_LEQUAL)
!    /* load the camera transformation matrix */
CALL argConvGlpara(patt_trans, gl_para)
CALL glMatrixMode(GL_MODELVIEW)
CALL glLoadMatrixd( gl_para )
CALL glEnable(GL_LIGHTING)
CALL glEnable(GL_LIGHT0)
CALL glLightfv(GL_LIGHT0, GL_POSITION, light_position)
CALL glLightfv(GL_LIGHT0, GL_AMBIENT, ambi)
CALL glLightfv(GL_LIGHT0, GL_DIFFUSE, lightZeroColor)
CALL glMaterialfv(GL_FRONT, GL_SPECULAR, mat_flash)
CALL glMaterialfv(GL_FRONT, GL_SHININESS, mat_flash_shiny)	
CALL glMaterialfv(GL_FRONT, GL_AMBIENT, mat_ambient)
CALL glMatrixMode(GL_MODELVIEW)

!--- from here -------------------
CALL glTranslatef( 0.0_glfloat, 0.0_glfloat, 30.0_glfloat )
DO i = 1, 3000
 CALL glTranslatef( dx(i) / 35.0, dy(i) / 35.0, dz(i) / 35.0 )
 IF (i > 100) THEN
  CALL glutSolidCube(2.0_gldouble )
 END IF 
END DO
!--- till here ----------------------

CALL glDisable( GL_LIGHTING )
CALL glDisable( GL_DEPTH_TEST )
RETURN
END SUBROUTINE draw
!-------------------------------------------------------------
END MODULE m_ART
!=====================================================
PROGRAM main1
USE m_art
IMPLICIT NONE
CALL glutInit()
CALL init()
CALL arVideoCapStart()
CALL argMainLoop( 0, C_FUNLOC(keyEvent), C_FUNLOC(mainLoop) )
STOP
END PROGRAM main1


MODULE m_bitmap
IMPLICIT NONE
INTEGER, PARAMETER :: DWORD = 4, LONG = 4, WORD = 2, kBYTE = 1 
TYPE :: T_BITMAPINFOHEADER
 SEQUENCE
 integer(DWORD) biSize ! knowns  DWORD 
 integer(LONG) biWidth ! knowns  LONG 
 integer(LONG) biHeight ! knowns  LONG 
 integer(WORD) biPlanes ! knowns  WORD 
 integer(WORD) biBitCount ! knowns  WORD 
 integer(DWORD) biCompression ! knowns  DWORD 
 integer(DWORD) biSizeImage ! knowns  DWORD 
 integer(LONG) biXPelsPerMeter ! knowns  LONG 
 integer(LONG) biYPelsPerMeter ! knowns  LONG 
 integer(DWORD) biClrUsed ! knowns  DWORD 
 integer(DWORD) biClrImportant ! knowns  DWORD 
END TYPE
!
TYPE :: T_RGBQUAD
 SEQUENCE
 integer(kBYTE) rgbBlue ! knowns  BYTE 
 integer(kBYTE) rgbGreen ! knowns  BYTE 
 integer(kBYTE) rgbRed ! knowns  BYTE 
 integer(kBYTE) rgbReserved ! knowns  BYTE 
END TYPE
!
TYPE :: T_BITMAPFILEHEADER
 SEQUENCE
 CHARACTER(2)   :: bfType
 INTEGER(DWORD) :: bfSize
 INTEGER(WORD)  :: bfReserved1
 INTEGER(WORD)  :: bfReserved2
 INTEGER(DWORD) :: bfOffBits
END TYPE 
CONTAINS
!----------------------------------------------------
SUBROUTINE ReadBitmapFile(text, bih, bcolmap, itmp)
IMPLICIT NONE
CHARACTER(*), INTENT(IN) :: text
TYPE (T_BITMAPINFOHEADER)     , INTENT(OUT) :: bih
TYPE (T_RGBQUAD), ALLOCATABLE, INTENT(OUT) :: bcolmap(:)
INTEGER(1)     , ALLOCATABLE, INTENT(OUT) :: itmp(:)
TYPE (T_BITMAPFILEHEADER) :: bfh
INTEGER :: i, len
INTEGER :: ir = 10
OPEN(ir, file = text, FORM = 'BINARY')
READ(ir) bfh
READ(ir) bih
len = bih%biWidth * bih%biHeight
IF (bih%biBitCount == 8) THEN 
 ALLOCATE(bcolmap(256))
 READ(ir) bcolmap
 ALLOCATE( itmp(len) )
 READ(ir) itmp
ELSE IF (bih%biBitCount == 4) THEN
 ALLOCATE(bcolmap(16))
 READ(ir) bcolmap
 ALLOCATE( itmp(len) )
 itmp = 0
 READ(ir) itmp
 DO i = 1, len, 2
  itmp( (i + 1) / 2) = itmp(i + 1) * 16 + itmp(i)   
 END DO
END IF 
CLOSE(ir)
RETURN
END SUBROUTINE ReadBitmapFile
!----------------------------------------------------
SUBROUTINE ReadBitMapData(text, iwidth, iheight, ipixel)
IMPLICIT NONE
CHARACTER(*), INTENT(IN) :: text
INTEGER, INTENT(OUT) :: iwidth, iheight
INTEGER(1), ALLOCATABLE, INTENT(OUT) :: ipixel(:)
TYPE (T_BITMAPINFOHEADER) :: bih
INTEGER :: istatus, isel
INTEGER :: i, j, k, m
INTEGER :: icounter, icolor
TYPE (T_RGBQUAD), ALLOCATABLE :: bcolmap(:)
INTEGER(1)     , ALLOCATABLE :: idata(:)
!
CALL ReadBitmapFile(text, bih, bcolmap, idata)
iwidth  = bih%biWidth
iheight = bih%biHeight
ALLOCATE( ipixel( 4 * iwidth * iheight ) )
IF (bih%biBitCount == 8) THEN
 DO j = 0, iheight - 1
  DO i = 0, iwidth - 1
   k = 4 * (j * iwidth + i) + 1
   m = idata(j * iwidth + i + 1)
   IF (m < 0) m = 256 + m
   m = m + 1 
   ipixel(k    ) = bcolmap(m)%rgbRed
   ipixel(k + 1) = bcolmap(m)%rgbGreen
   ipixel(k + 2) = bcolmap(m)%rgbBlue
   ipixel(k + 3) = 255
  END DO
 END DO
ELSE IF (bih%biBitCount == 4) THEN
 icounter = 0
 DO j = 0, iheight - 1
  DO i = 0, iwidth - 1
   icounter = icounter + 1
   icolor = idata(icounter)
!   IF (icolor < 0) icolor = icolor + 256 ! 
   ipixel(k    ) = bcolmap(icolor)%rgbRed
   ipixel(k + 1) = bcolmap(icolor)%rgbGreen
   ipixel(k + 2) = bcolmap(icolor)%rgbBlue
   ipixel(k + 3) = 255
  END DO
 END DO
END IF
DEALLOCATE(bcolmap, idata)
RETURN
END SUBROUTINE ReadBitMapData
!----------------------------------------------------
END MODULE m_bitmap
!====================================================
MODULE m_interface
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
!---- TYPE definition ---------------------------------
TYPE, BIND(C) :: t_ARParam
 INTEGER(C_INT) :: ixsize, iysize
 REAL(C_DOUBLE) :: dmat(4, 3)
 REAL(C_DOUBLE) :: dist_factor(4)
END TYPE
!
TYPE, BIND(C) :: t_ARMarkerInfo 
 INTEGER(C_INT) :: iarea
 INTEGER(C_INT) :: id
 INTEGER(C_INT) :: dir
 REAL(C_DOUBLE)  :: cf
 REAL(C_DOUBLE)  :: pos(2)
 REAL(C_DOUBLE)  :: dline(3, 4)
 REAL(C_DOUBLE)  :: vertex(2, 4)
END TYPE
!---- global variable ---------------------------------
INTEGER(C_INT) :: ithresh = 100
INTEGER(C_INT) :: icount = 0
INTEGER(C_INT) :: ipatt_id
REAL(C_DOUBLE) :: patt_width = 80.0_c_double
REAL(C_DOUBLE) :: patt_center(2) = (/0.0_c_double, 0.0_c_double/)
REAL(C_DOUBLE) :: patt_trans(3, 4)
INTEGER(C_INT) :: ixsize, iysize
TYPE(t_ARParam) :: cparam
!---- INTERFACE BLOCK -----------------------------------
INTERFACE
 SUBROUTINE arVideoCapStop() BIND(C, NAME = 'arVideoCapStop')
 END SUBROUTINE arVideoCapStop
! 
 SUBROUTINE arVideoClose() BIND(C, NAME = 'arVideoClose')
 END SUBROUTINE arVideoClose
!
 SUBROUTINE argCleanup() BIND(C, NAME = 'argCleanup')
 END SUBROUTINE argCleanup
! 
 INTEGER(C_INT) FUNCTION arVideoGetImage() BIND(C, NAME = 'arVideoGetImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE arUtilSleep(i) BIND(C, NAME = 'arUtilSleep')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: i 
 END
!
 SUBROUTINE arUtilTimerReset() BIND(C, NAME = 'arUtilTimerReset')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE argDrawMode2D() BIND(C, NAME = 'argDrawMode2D')
 END 
!
 SUBROUTINE argDispImage( iptr, i, j ) BIND(C, NAME = 'argDispImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER, VALUE :: i, j
 END 
!
 INTEGER(C_INT) FUNCTION arDetectMarker( iptr, i, mark, k ) BIND(C, NAME = 'arDetectMarker')
 USE, INTRINSIC :: ISO_C_BINDING 
 import
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER(C_INT)     , VALUE :: i
 TYPE(t_ARMarkerInfo), POINTER :: mark(:)
 INTEGER(C_INT) :: k
 END 
!
 SUBROUTINE arVideoCapNext() BIND(C, NAME = 'arVideoCapNext')
 END 
!
 SUBROUTINE argSwapBuffers() BIND(C, NAME = 'argSwapBuffers')
 END 
!
 SUBROUTINE arGetTransMat(ptr, d1, d, d2) BIND(C, NAME = 'arGetTransMat')
 USE, INTRINSIC :: ISO_C_BINDING 
 TYPE(C_PTR), VALUE :: ptr
 REAL(C_DOUBLE), VALUE :: d
 REAL(C_DOUBLE):: d1(2), d2(3, 4)
 END 
!
 SUBROUTINE argDrawMode3D() BIND(C, NAME = 'argDrawMode3D')
 END SUBROUTINE argDrawMode3D
 !
 SUBROUTINE argDraw3dCamera(ix, iy) BIND(C, NAME = 'argDraw3dCamera')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: ix, iy 
 END
 !
 SUBROUTINE argConvGlpara(x, y) BIND(C, NAME = 'argConvGlpara')
 USE, INTRINSIC :: ISO_C_BINDING 
 REAL(C_DOUBLE) :: x(4, 3), y(4, 4) 
 END
! 
 SUBROUTINE arVideoCapStart() BIND(C, NAME = 'arVideoCapStart')
 END
!
 SUBROUTINE argMainLoop( i, j, k ) BIND(C, NAME = 'argMainLoop')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER, VALUE :: i
 TYPE(C_FUNPTR), VALUE :: j, k
 END
!
 INTEGER(4) FUNCTION arVideoOpen(text) BIND(C, NAME = 'arVideoOpen')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(c_char) :: text(*)
 END
!
 INTEGER(C_INT) FUNCTION arVideoInqSize(ix, iy) BIND(C, NAME = 'arVideoInqSize')
 USE, INTRINSIC :: ISO_C_BINDING
 INTEGER(C_INT) :: ix, iy
 END
! 
 INTEGER(C_INT) FUNCTION arParamLoad(text, i, xparam) BIND(C, NAME = 'arParamLoad')
 USE, INTRINSIC :: ISO_C_BINDING
 import :: t_ARParam
 CHARACTER(c_char) :: text(*)
 TYPE(t_ARParam) :: xparam
 INTEGER(C_INT), VALUE :: i
 END
!
 SUBROUTINE arParamChangeSize(xparam, i, j, yparam) BIND(C, NAME = 'arParamChangeSize')
 USE, INTRINSIC :: ISO_C_BINDING
 import :: t_ARParam
 TYPE(t_ARParam) :: xparam, yparam
 INTEGER(C_INT), VALUE :: i, j
 END
! 
 SUBROUTINE arInitCparam(xparam) BIND(C, NAME = 'arInitCparam' )
 import :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
! 
 SUBROUTINE arParamDisp(xparam) BIND(C, NAME = 'arParamDisp')
 import :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
!
 INTEGER(4) FUNCTION arLoadPatt(text) BIND(C, NAME = 'arLoadPatt')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(c_char) :: text(*)
 END
!
 SUBROUTINE argInit(iptr, d, i1, i2, i3, i4) BIND(C, NAME = 'argInit')
 USE, INTRINSIC :: ISO_C_BINDING
 INTEGER(C_INTPTR_T), VALUE :: iptr
 REAL(C_DOUBLE)      , VALUE :: d
 INTEGER(C_INT)     , VALUE :: i1, i2, i3, i4
 END
!
 REAL(8) FUNCTION arUtilTimer() BIND(C, NAME = 'arUtilTimer')
 END FUNCTION arUtilTimer
END INTERFACE
END MODULE m_interface
!=========================================================================
MODULE m_ART
USE :: m_interface
LOGICAL :: qaniki = .FALSE.
CONTAINS
!-------------------------------------------------
SUBROUTINE init() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
IMPLICIT NONE
TYPE(t_ARParam) :: wparam
CHARACTER (LEN = 35) :: vconf = "Data\\WDM_camera_flipV.xml" // c_null_char
CHARACTER (LEN = 35) :: cparam_name = "Data/camera_para.dat" // c_null_char
CHARACTER (LEN = 35) :: patt_name   = "Data/patt.hiro"       // c_null_char
!   /* open the video path */
IF( arVideoOpen( vconf ) < 0 ) STOP
!    /* find the size of the window */
IF( arVideoInqSize( ixsize, iysize ) < 0 ) STOP
PRINT *, "Image size (x,y) = ", ixsize, iysize
!    /* set the initial camera parameters */
IF ( arParamLoad( cparam_name, 1, wparam ) < 0 ) THEN
 PRINT *, "Camera parameter load error !!"
 STOP
END IF
CALL arParamChangeSize( wparam, ixsize, iysize, cparam )
CALL arInitCparam( cparam )
PRINT *, "*** Camera Parameter ***"
CALL arParamDisp( cparam )
ipatt_id = arLoadPatt( patt_name )
IF ( ipatt_id < 0 ) THEN
 PRINT *, "pattern load error !!"
 STOP
END IF
!    /* open the graphics window */
CALL argInit( LOC(cparam), 1.0_c_double, 0, 0, 0, 0 )
RETURN
END SUBROUTINE init
!-------------------------------------------------
SUBROUTINE mainLoop() BIND(C, NAME = 'mainLoop')
IMPLICIT NONE
INTEGER(C_INTPTR_T) :: dataPtr
TYPE(t_ARMarkerInfo), POINTER:: marker_info(:)
INTEGER(C_INT) :: marker_num
INTEGER(C_INT) :: j, k
LOGICAL, SAVE :: qfirst = .TRUE.
IF (qfirst) THEN 
 ALLOCATE(marker_info(50)) ! ?
 qfirst = .FALSE.
END IF
dataPtr = arVideoGetImage()
IF (dataPtr == 0) THEN
 CALL arUtilSleep(2)
 RETURN
END IF
IF (icount == 0) CALL arUtilTimerReset()
icount = icount + 1
CALL argDrawMode2D()
CALL argDispImage( dataPtr, 0, 0 )
IF ( arDetectMarker( dataPtr, ithresh, marker_info, marker_num ) < 0 ) THEN
 CALL cleanup()
 RETURN
END IF
CALL arVideoCapNext()
k = -1
DO j = 1, marker_num
 IF (ipatt_id == marker_info(j)%id) THEN
  IF (k == -1) THEN 
   k = j
  ELSE 
   IF (marker_info(k)%cf < marker_info(j)%cf ) k = j
  END IF
 END IF
END DO
IF ( k == -1 ) THEN 
 CALL argSwapBuffers()
 RETURN
END IF
! draw-main 
CALL arGetTransMat( C_LOC(marker_info(k)), patt_center, patt_width, patt_trans )
CALL draw()
CALL argSwapBuffers()
RETURN
END SUBROUTINE mainLoop
!-------------------------------------------------
SUBROUTINE cleanup() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING
IMPLICIT NONE 
CALL arVideoCapStop()
CALL arVideoClose()
CALL argCleanup()
RETURN
END SUBROUTINE cleanup
!-------------------------------------------------
SUBROUTINE keyEvent( key, ix, iy ) BIND(C, NAME = 'keyEvent')
USE, INTRINSIC :: ISO_C_BINDING 
IMPLICIT NONE
INTEGER(C_SIGNED_CHAR), VALUE :: key
INTEGER(C_INT)        , VALUE :: ix, iy
SELECT CASE (key)
 CASE (27) 
  PRINT *, "*** ", icount / arUtilTimer(), "(frame/sec)"
  CALL cleanup()
  STOP
 CASE (32)
  qaniki = .TRUE. 
 CASE DEFAULT
  CONTINUE
END SELECT
RETURN
END SUBROUTINE keyEvent
!---------------------------------------------------------------
SUBROUTINE draw() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
IMPLICIT NONE
REAL(gldouble) :: gl_para(4, 4)
REAL(glfloat ) :: mat_ambient(4)  = (/0.0_glfloat, 0.0_glfloat, 1.0_glfloat, 1.0_glfloat/)
REAL(glfloat ) :: mat_flash(4)    = (/0.0_glfloat, 0.0_glfloat, 1.0_glfloat, 1.0_glfloat/)
REAL(glfloat ) :: mat_emission(4) = (/0.1_glfloat, 0.0_glfloat, 1.0_glfloat, 1.0_glfloat/)
REAL(glfloat ) :: mat_flash_shiny(1) = (/50.0_glfloat/)
REAL(glfloat ) :: light_position(4)  = (/0.0_glfloat, 0.0_glfloat, 200.0_glfloat, 0.0_glfloat/)
REAL(glfloat ) :: ambi(4)            = (/0.1_glfloat, 0.1_glfloat, 0.1_glfloat, 0.1_glfloat/)
REAL(glfloat ) :: lightZeroColor(4)  = (/0.9_glfloat, 0.9_glfloat, 0.9_glfloat, 0.1_glfloat/)
!    
CALL argDrawMode3D()
CALL argDraw3dCamera( 0, 0 )
CALL glClearDepth( 1.0_gldouble )
CALL glClear(GL_DEPTH_BUFFER_BIT)
CALL glEnable(GL_DEPTH_TEST)
CALL glDepthFunc(GL_LEQUAL)
!    /* load the camera transformation matrix */
CALL argConvGlpara(patt_trans, gl_para)
CALL glMatrixMode(GL_MODELVIEW)
CALL glLoadMatrixd( gl_para )
CALL glEnable(GL_LIGHTING)
CALL glEnable(GL_LIGHT0)
CALL glLightfv(GL_LIGHT0, GL_POSITION, light_position)
CALL glLightfv(GL_LIGHT0, GL_AMBIENT, ambi)
CALL glLightfv(GL_LIGHT0, GL_DIFFUSE, lightZeroColor)
CALL glMaterialfv(GL_FRONT, GL_SPECULAR, mat_flash)
CALL glMaterialfv(GL_FRONT, GL_SHININESS, mat_flash_shiny)	
CALL glMaterialfv(GL_FRONT, GL_AMBIENT, mat_ambient)
CALL glMatrixMode(GL_MODELVIEW)
CALL glPushMatrix()
!
!
IF (qaniki) THEN 
 CALL glTranslatef( 0.0_glfloat, 0.0_glfloat, 130.0_glfloat )
 CALL glPushMatrix()
 CALL glScalef( 2.0_glfloat, 2.0_glfloat, 2.0_glfloat )
! core 
 lightZeroColor  = (/1.0_glfloat, 0.0_glfloat, 0.1_glfloat, 0.1_glfloat/)
 CALL glLightfv(GL_LIGHT0, GL_DIFFUSE, lightZeroColor)
 CALL glTranslatef(-10.0_glfloat, 0.0_glfloat, -19.0_glfloat )
 CALL glutSolidSphere( 3.0_gldouble, 10, 10 )
 CALL glTranslatef( 38.0_glfloat, 0.0_glfloat, -5.0_glfloat )
 CALL glutSolidSphere( 3.0_gldouble, 10, 10 )
 CALL glPopMatrix()
!
 CALL glScalef( 80.0_glfloat, 80.0_glfloat, 80.0_glfloat )
 CALL glEnable(GL_TEXTURE_2D)
 CALL glBindTexture( GL_TEXTURE_2D, 3 )
 CALL glScalef( 4.0_glfloat, 4.0_glfloat, 4.0_glfloat )
!
 CALL glBegin( GL_QUADS )
  CALL glTexCoord2f( 0.0_glfloat, 0.0_glfloat )
  CALL glVertex3f  (-0.5_glfloat, 0.0_glfloat,-0.5_glfloat )
  CALL glTexCoord2f( 0.0_glfloat, 1.0_glfloat )
  CALL glVertex3f  (-0.5_glfloat, 0.0_glfloat, 0.5_glfloat )
  CALL glTexCoord2f( 1.0_glfloat, 1.0_glfloat )
  CALL glVertex3f  ( 0.5_glfloat, 0.0_glfloat, 0.5_glfloat )
  CALL glTexCoord2f( 1.0_glfloat, 0.0_glfloat )
  CALL glVertex3f  ( 0.5_glfloat, 0.0_glfloat,-0.5_glfloat )
 CALL glEnd()
 CALL glDisable(GL_TEXTURE_2D)
ELSE

 CALL glTranslatef( 0.0_glfloat, 0.0_glfloat, 160.0_glfloat )
 CALL glScalef( 80.0_glfloat, 80.0_glfloat, 80.0_glfloat )
 CALL glutSolidOctahedron()
 CALL glPushAttrib(GL_LIGHTING_BIT)
 CALL glLightfv(GL_LIGHT0, GL_EMISSION, ambi)
 CALL glMaterialfv(GL_FRONT, GL_EMISSION, mat_emission)
 CALL glutWireOctahedron()
 CALL glPopAttrib()
END IF
!
CALL glPopMatrix()
CALL glDisable( GL_LIGHTING )
CALL glDisable( GL_DEPTH_TEST )
RETURN
END SUBROUTINE draw
!-------------------------------------------------------------
SUBROUTINE makeImage()
USE :: m_bitmap
IMPLICIT NONE
INTEGER :: iheight, iwidth
INTEGER(1), ALLOCATABLE :: image(:)
INTEGER :: iret
CALL ReadBitMapData("aniki2b.bmp", iwidth, iheight, image)
CALL glPixelStorei( GL_UNPACK_ALIGNMENT, 1 )
CALL glBindTexture( GL_TEXTURE_2D, 3 )
CALL glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_WRAP_S    , GL_REPEAT  )
CALL glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_WRAP_T    , GL_REPEAT  )
CALL glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST )
CALL glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST )
CALL glTexImage2D( GL_TEXTURE_2D, 0, 4, iwidth, iheight,0, &
                    GL_RGBA, GL_UNSIGNED_BYTE, image )
CALL glTexEnvi( GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL )
CALL glEnable ( GL_TEXTURE_2D )
RETURN
END SUBROUTINE makeImage
!---------------------------------------------------------
END MODULE m_ART
!======================================================================
PROGRAM main1
USE m_art
IMPLICIT NONE
CALL glutInit()
CALL init()
CALL makeImage()
CALL arVideoCapStart()
CALL argMainLoop( 0, C_FUNLOC(keyEvent), C_FUNLOC(mainLoop) )
STOP
END PROGRAM main1