!***********************************************************************
!	両側から粒子に押されている隔壁の分布を調べる
!		反射[ *左粒子 |隔壁 *右粒子 ]反射
!	隔壁の位置のヒストグラムをとったところ、やはり中央付近にある確率が高い。
! 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		!１回でも衝突があれば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
				! 一方向に偏らないように、１回ごとに走査する向きを変える
				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回に１回サンプリング
			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 *******************************************************