!***********************************************************************
!	両側から粒子に押されている隔壁の分布を調べる
!		反射[ *左粒子 |隔壁 *右粒子 ]反射
!	両端に乱数を導入して分布の様子を調べる
! Date: 2007/02/14
! Compile:
! g95 MaxDemon04.F90 mtfort90.o -Wno=165 -Wall -O -o MaxDemon04 -Wl,--subsystem,console -lGrWin -mwindows
! PreProcess:
! -D=LIMIT シミュレーション上限指定
! -D=NOGR  グラフィック無し

! 乱数パターンの選択
! #define	RANDPAT_A	1
! #define	RANDPAT_B	1
! #define	RANDPAT_C	1
#define	RANDPAT_D	1

PROGRAM MaxDemon04
USE mtmod	! MT乱数モジュール使用
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, PARAMETER :: AVE = 1.0
	
	! 物体の初期位置
	REAL :: x(1:N_BALL)  = (/125.0, 250.0, 375.0/)
	! 物体の初速
	REAL :: dx(1:N_BALL) = (/AVE, 0.0, AVE/)
	! 物体の質量
	REAL :: m(1:N_BALL)  = (/1.0, 5.0, 1.0/)
	
	! ワープゾーンの位置
	INTEGER, PARAMETER :: P1FROM =  50, P1TO = 100
	INTEGER, PARAMETER :: P2FROM = 225, P2TO = 275
	
	! ワープ用作業変数
	INTEGER, PARAMETER :: WARP_INTERVAL = 1500
		! 150 -- それでも気になるのでピコピコが生じないほど長くしてみた
		! 20 -- インターバル時間を変えても偏り傾向には影響ない
	INTEGER :: interval=0, warpcnt=0, totwarpcnt=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
	CALL sgrnd(4357)	! ランダムシード
#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
#ifndef NOGR
		! ワープゾーンを描く
		CALL GWsrect(IR, real(P1FROM), 0.0, real(P1TO),real(HEIGHT), 12)
		CALL GWsrect(IR, real(P2FROM), 0.0, real(P2TO),real(HEIGHT), 12)
#endif
		! 衝突ループ
		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	! 方向フラグ反転
		
		! 順序が保たれているかどうかチェック
		IF ( x(1) > x(2) ) THEN
		!!	PRINT *, 'ERROR: 順序違反1:2'
			x(1) = x(1) - dx(1)
		!!	x(2) = x(2) - dx(2)
		END IF
		IF ( x(2) > x(3) ) THEN
		!!	PRINT *, 'ERROR: 順序違反2:3'
			x(2) = x(2) - dx(2)
		!!	x(3) = x(3) - dx(3)
		END IF
		
		! 両端の壁への衝突判定
		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)		! 反射
				! 反射時にできるだけ平均速度に近づける
				CALL AveDX( i, dx )
			END IF
		END DO
		
		!ワープ処理 -- x(2)を２つのセグメント間で移動する
		IF ( interval == 0 ) THEN
			! 順方向ワープ
			IF( P2FROM < x(2) .and. x(2) < P2TO &
&				.AND. (x(1) < P1FROM) &
&				.AND. (P2TO < x(3)) ) THEN
				x(2) = x(2) - (P2FROM - P1FROM)
				warpcnt = warpcnt + 1
				totwarpcnt = totwarpcnt + 1
				PRINT "('+++ WARP +++', I4, I6)", warpcnt, totwarpcnt
			! 逆方向ワープ -- ここにelseを入れないと対称にならないぞ
			ELSE IF( P1FROM < x(2) .and. x(2) < P1TO &
&				.AND. (x(1) < P1FROM) &
&				.AND. (P2TO < x(3)) ) THEN
				x(2) = x(2) + (P2FROM - P1FROM)
				warpcnt = warpcnt - 1
				totwarpcnt = totwarpcnt + 1
				PRINT "('--- WARP ---', I4, I6)", warpcnt, totwarpcnt
			END IF
			! ワープ後しばらくの間は次のワープができない
			interval = WARP_INTERVAL
		ELSE
			interval = interval - 1
		END IF
		
#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 >= 5000 ) THEN ! 5000回に１回サンプリング
	!!		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
		! 1000回でシミュレーションを打ち切る
	!	IF (tothistcnt > 3000) THEN
		IF (totwarpcnt > 1000) 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

!***********************************************************************
! 反射時に速度を平均に近づける
! i : ボールの番号 in
! dx() : 衝突前後の速度 in-out
! パラメータ AVE : 平均速度
SUBROUTINE AveDX( i, dx )
IMPLICIT NONE
	INTEGER, INTENT(IN) :: i
	REAL, INTENT(INOUT) :: dx(:)
	
	! エネルギーの流出入を積算で保持する
	REAL, SAVE :: transeg(1:2) = (/0.0, 0.0/)
	REAL :: diff = 0.0, rndave
!	REAL :: dx0, eg

#if defined(RANDPAT_A)
	rndave = real(AVE) + real( (grnd() - grnd()) / 2.0)	! 乱数で上下に振る
#elif defined(RANDPAT_B)
	rndave = real(AVE) + real( (grnd() - grnd()) )	! 分散を大きく
#elif defined(RANDPAT_C)
	rndave = real(AVE) + real( grnd() * 2.0 - 1.0 )	! 分散を思い切り大きく
#elif defined(RANDPAT_D)
	rndave = GenVelo() / 1.125	! Maxwellの速度分布、速度の平均値は実測
#else
	rndave = real(AVE)	! 乱数無し
#endif
	
!	PRINT "('rndave=', F7.3)", rndave
	
	IF ( dx(i) > 0.0 ) THEN
	!	dx0 = dx(i)
		diff = (rndave - dx(i)) / 2.0	! 平均に近づける
		dx(i) = dx(i) + diff
		
	!	PRINT "('From=', F7.3 ,' To=', F7.3)", dx0, dx(i)
		
		! 反射時に移動したエネルギー量積算
	!	eg = (dx(i) ** 2 - dx0 ** 2) / 2.0
	!	transeg(1) = transeg(1) + eg
		
	ELSE IF ( dx(i) < 0.0 ) THEN
	!	dx0 = dx(i)
		diff = (rndave - (-dx(i))) / 2.0	! 平均に近づける
		dx(i) = dx(i) - diff
		
	!	PRINT "('From=', F7.3 ,' To=', F7.3)", dx0, dx(i)
		
		! 反射時に移動したエネルギー量積算
	!	eg = (dx(i) ** 2 - dx0 ** 2) / 2.0
	!	transeg(2) = transeg(2) + eg
	END IF
	
	! 累積したエネルギーを表示
!2000	FORMAT( '(', I1 ,')  ', F7.3, '  <', F7.3, ' : ', F7.3, '>')
!		WRITE(*,2000) i, diff, transeg(1), transeg(2)
	
	RETURN
END SUBROUTINE AveDX

!***********************************************************************
! 擬似的な気体速度分布の生成
! Maxwellの速度分布を破棄法にて求める
REAL FUNCTION GenVelo() RESULT(rx)
IMPLICIT NONE
	REAL, PARAMETER :: X_MAX = 4.0
	REAL, PARAMETER :: Y_MAX = 0.4
	REAL :: ry
	DO
		rx = real(grnd() * X_MAX)
		ry = real(grnd() * Y_MAX)
		IF ( ry < rx * rx * exp(- rx * rx) ) THEN
			RETURN
		END IF
	END DO
END FUNCTION GenVelo

END PROGRAM MaxDemon04
!*** END of FILE *******************************************************