!*********************************************************************** ! 両側から粒子に押されている隔壁の分布を調べる ! 反射[ *左粒子 |隔壁 *右粒子 ]反射 ! 隔壁の位置のヒストグラムをとったところ、やはり中央付近にある確率が高い。 ! Date: 2007/02/13 ! Compile: ! g95 MaxDemon01.F90 -Wno=165 -Wall -O -o MaxDemon01 -Wl,--subsystem,console -lGrWin -mwindows ! PreProcess: ! -D=LIMIT シミュレーション上限指定 ! -D=NOGR グラフィック無し PROGRAM MaxDemon01 IMPLICIT NONE ! GWライブラリについては以下の警告が出る ! Warning (165): Implicit interface 'gwopen' called at (1) ! 気になるのであれば、-Wno=165 として警告を抑制する ! 物体の数、1=ボール、2=隔壁、3=ボール INTEGER, PARAMETER :: N_BALL = 3 ! 画面サイズ INTEGER, PARAMETER :: WIDTH = 500, HEIGHT = 300 ! 描画する物体の高さ、大きさ REAL, PARAMETER :: yy = 150.0, rr = 8.0 ! 物体の初期位置 REAL :: x(1:N_BALL) = (/125.0, 250.0, 375.0/) ! 物体の初速 REAL :: dx(1:N_BALL) = (/1.0, 0.0, 1.0/) ! REAL :: dx(1:N_BALL) = (/5.0, 0.0, 5.0/) ! 試しに初速を変えてみる ! 物体の質量 REAL :: m(1:N_BALL) = (/1.0, 5.0, 1.0/) ! REAL :: m(1:N_BALL) = (/1.0, 1.23456789, 1.0/) ! 試しに質量を変えてみる ! 当たり判定フラグ LOGICAL :: colflag(1:N_BALL-1) LOGICAL :: acolflag !1回でも衝突があればtrue REAL :: next_x ! 最大繰り返し計算回数 INTEGER, PARAMETER :: MAXITER = 100 ! ヒストグラム INTEGER, PARAMETER :: N_HIST = 20 INTEGER :: hist(1:N_HIST) = (N_HIST * 0) INTEGER :: histcnt = 0, pos #ifdef LIMIT INTEGER :: tothistcnt = 0 #endif INTEGER :: i, iter #ifndef NOGR INTEGER :: IR, NOG = 1 !GWライブラリ関連 #endif INTEGER :: direct = 1 ! 方向フラグ ! BEGIN #ifndef NOGR CALL GWopen(IR,0) CALL GWindow(IR, -20.0, -20.0, real(WIDTH)+20.0, real(HEIGHT)+20.0 ) CALL GWsetogn(IR, 0, -NOG) ! Object Group #endif DO ! 衝突ループ DO iter=1, MAXITER ! 衝突直後の衝突がなくなるまで、つまり同時衝突がなくなるまで ! このめんどくさい処理がないと「すり抜け」が生じてしまう ! 当たり判定 acolflag = .FALSE. DO i=1, N_BALL-1 colflag(i) = isCol( x, dx, i, i+1 ) IF ( colflag(i) ) THEN ! PRINT "('Hit! (',I2,')',F6.2, F6.2, F6.2, F6.2)", & !& iter, x(i), dx(i), x(i+1), dx(i+1) acolflag = .TRUE. END IF END DO IF ( .NOT. acolflag ) THEN ! 衝突していなければ EXIT !ループ終了 END IF ! 衝突処理 ! * 1:2, 2:3 の衝突 IF ( direct > 0 ) THEN ! 一方向に偏らないように、1回ごとに走査する向きを変える DO i=1, N_BALL-1 IF ( colflag(i) ) THEN CALL Collision( dx, m, i, i+1 ) END IF END DO ELSE DO i=N_BALL-1, 1, -1 IF ( colflag(i) ) THEN CALL Collision( dx, m, i, i+1 ) END IF END DO END IF END DO IF (iter >= MAXITER ) THEN PRINT *, 'ERROR: 衝突無限ループ' stop END IF ! 物体の移動 DO i=1, N_BALL x(i) = x(i) + dx(i) END DO direct = - direct ! 方向フラグ反転 ! 両端の壁への衝突判定 DO i=1, N_BALL next_x = x(i) + dx(i) ! 上限/下限を超えたら IF ( next_x > real(WIDTH) .OR. next_x < 0.0 ) THEN dx(i) = - dx(i) ! 反射 END IF END DO #ifndef NOGR ! 描画 ! ボールを描く CALL GWsetpen(IR, 2, -1, -1, -1) ! 2=暗緑色 CALL GWsetbrs(IR, 2, 1, -1) !塗りつぶし DO i=1, N_BALL, 2 CALL GWellipse(IR, x(i)-rr, yy-rr, x(i)+rr, yy+rr ) END DO ! 壁を描く DO i=2, N_BALL, 2 CALL GWsrect(IR, x(i)-rr, yy-30.0, x(i)+rr, yy+30.0, 13 ) !13=赤 END DO ! 外枠を描く CALL GWsetpen(IR, 4, -1, -1, -1) ! 4=濃紺 CALL GWsetbrs(IR, -1, 0, -1) ! 中空 CALL GWrect(IR, 0.0, 50.0, real(WIDTH), real(HEIGHT)-50.0 ) CALL GWflush(IR, -NOG) #endif ! ヒストグラムをとる histcnt = histcnt + 1 IF ( histcnt >= 100 ) THEN ! 100回に1回サンプリング histcnt = 0 ! x(2)の位置を把握する pos = int( x(2) / (WIDTH / N_HIST) ) + 1 hist(pos) = hist(pos) + 1 PRINT *, hist(1:N_HIST) #ifdef LIMIT tothistcnt = tothistcnt + 1 #endif END IF #ifdef LIMIT ! 3000回でシミュレーションを打ち切る IF (tothistcnt > 3000) THEN EXIT END IF #endif END DO #ifndef NOGR CALL GWSETBK(IR, 1) CALL GWREFRESH(IR) CALL GWQUIT(IR) #endif CONTAINS !*********************************************************************** ! 衝突判定 ! 返り値: .FALSE.=はずれ、.TRUE.=あたり LOGICAL FUNCTION isCol( x, dx, i, j ) IMPLICIT NONE INTEGER, INTENT(IN) :: i, j REAL, INTENT(IN) :: x(:), dx(:) REAL :: now_x1, now_x2, next_x1, next_x2 now_x1 = x(i) now_x2 = x(j) next_x1 = x(i) + dx(i) next_x2 = x(j) + dx(j) ! now_は等しくても衝突しない、こうしないと繰り返し衝突が起こる ! next_は等しければ衝突する IF ( (now_x1 > now_x2 .AND. next_x1 <= next_x2) & & .OR. (now_x1 < now_x2 .AND. next_x1 >= next_x2) ) THEN isCol = .TRUE. ELSE isCol = .FALSE. END IF RETURN END FUNCTION isCol !*********************************************************************** ! 衝突処理 ! dx() : 衝突前後の速度 in-out ! m() : 質量 in ! i, j : 対象となる質点の番号 in SUBROUTINE Collision( dx, m, i, j) IMPLICIT NONE REAL, INTENT(INOUT) :: dx(:) REAL, INTENT(IN) :: m(:) INTEGER, INTENT(IN) :: i, j REAL :: dwork ! PRINT "('Collision1 i=',I2,':',F7.2,' j=',I2,':',F7.2)", & !& i, dx(i), j, dx(j) ! 衝突前の運動量、エネルギー ! PRINT "('P1,E1 = ',F7.2,', ',F7.2)", & !& m(i) * dx(i) + m(j) * dx(j), & !& m(i) ** dx(i) + m(j) ** dx(j) dwork = ( dx(i) * m(i) - dx(i) * m(j) + 2.0 * dx(j) * m(j) ) & & / ( m(i) + m(j) ) dx(j) = ( 2.0 * dx(i) * m(i) - dx(j) * m(i) + dx(j) * m(j) ) & & / ( m(i) + m(j) ) dx(i) = dwork ! PRINT "('Collision2 i=',I2,':',F7.2,' j=',I2,':',F7.2)", & !& i, dx(i), j, dx(j) ! 衝突後の運動量、エネルギー ! PRINT "('P2,E2 = ',F7.2,', ',F7.2)", & !& m(i) * dx(i) + m(j) * dx(j), & !& m(i) ** dx(i) + m(j) ** dx(j) ! PRINT "('COL_velo[',I1,'][',I1,'], ',F7.2,', ',F7.2)", i, j, dx(i), dx(j) RETURN END SUBROUTINE Collision END PROGRAM MaxDemon01 !*** END of FILE *******************************************************