fortran66のブログ

fortran について書きます。

Rutherford Scattering 2

ARToolKit

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

■計算結果

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

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

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

裏から見た図。

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

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


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

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

■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