!***********************************************************************
!	熱源モジュール	Thermal Pool
!	多数の分子間でランダムにエネルギーを交換し、仮想的な熱源を実現する
!	分子を２個だけにしてみた
!	質の悪い乱数（線形合同法）によって周期性を見る
!	Compile:
!		g95 -c thpool_N2_lrand.f90
!		g95 thpool_N2_lrand.o lrand.o -o thpool_N2_lrand
!	Date: 2007/07/15

MODULE thpool_N2_lrand
USE lrand	! 線形合同乱数モジュール使用

!INTEGER, PARAMETER :: N_MOL = 6	! 分子の数
INTEGER, PARAMETER :: N_MOL = 2	! 分子の数
REAL*8, PARAMETER :: MAX_EXCHANGE = 2.5d0	! １回に交換する最大値

REAL*8 :: Mol(1:N_MOL)

 CONTAINS
!**********************************************************************
! 初期化
!	初期平均エネルギー
!	初期の全エネルギーを返す
REAL*8 FUNCTION pool_init(e0)
IMPLICIT NONE
	REAL*8, INTENT(IN) :: e0
	
	CALL l_srand(8)	! ランダムシード
	Mol = e0	! 配列の全要素を初期化
	CALL pool_run( 1000 )	! 十分に安定するまで回す

	pool_init = N_MOL * e0	! 最初に与えた全エネルギーを返す
	return
END FUNCTION pool_init

!**********************************************************************
! 熱浴内でのエネルギー交換を進める
! 1stepで単純に N_MOLの交換を行う
SUBROUTINE pool_run(step)
IMPLICIT NONE
	INTEGER, INTENT(IN) :: step
	INTEGER :: st, i, j
	REAL*8 :: chg
	
	DO st = 0, step * N_MOL
		! 異なる２つの分子 i, j を選び出す
		i = mod(l_rand(), N_MOL) + 1
		DO
			j = mod(l_rand(), N_MOL) + 1
			IF (i /= j) THEN
				EXIT
			END IF
		END DO
		! i -> j にエネルギーを渡す
		! 交換する大きさを決める
		chg = l_rand_max(MAX_EXCHANGE)
		IF (chg > Mol(i)) THEN	! その大きさを交換できなければ
			Mol(j) = Mol(j) + Mol(i)	! iを全てjに渡す
			Mol(i) = 0.0
		ELSE
			Mol(j) = Mol(j) + chg
			Mol(i) = Mol(i) - chg
		END IF
	END DO
	
END SUBROUTINE pool_run

!**********************************************************************
! エネルギーの受け渡しを行う
REAL*8 FUNCTION pool_get(i)
IMPLICIT NONE
	INTEGER, INTENT(IN) :: i	! 対象となる分子
	pool_get = Mol(i)
END FUNCTION pool_get

SUBROUTINE pool_set(i,e0)
IMPLICIT NONE
	INTEGER, INTENT(IN) :: i	! 対象となる分子
	REAL*8, INTENT(IN) :: e0
	Mol(i) = e0
END SUBROUTINE pool_set

!**********************************************************************
! 熱浴の全エネルギーを返す
REAL*8 FUNCTION pool_total()
IMPLICIT NONE
	INTEGER :: i
	REAL*8 :: sum
	
	sum = 0.0d0
	DO i=0, N_MOL
		sum = sum + Mol(i)
	END DO
	pool_total = sum
	return
END FUNCTION pool_total

!**********************************************************************
! 熱浴内の様子を出力する
SUBROUTINE pool_list
IMPLICIT NONE
	INTEGER, PARAMETER :: N_HIST = 20
	INTEGER :: hist(1:N_HIST)
	INTEGER :: i, pos
	REAL*8 :: sum
	REAL*8, PARAMETER :: DELTA = 0.5d0	! ヒストグラムの刻み幅
	
	hist = 0	! 配列の0クリアー
	sum = 0.0
	DO i=0, N_MOL
		pos = int( Mol(i) / DELTA ) + 1
		IF (pos > N_HIST) THEN
			pos = N_HIST
		END IF
		hist(pos) = hist(pos) + 1
		sum = sum + Mol(i)
	END DO
	PRINT *, hist(1:N_HIST), ':', sum

END SUBROUTINE pool_list

END MODULE thpool_N2_lrand

!**********************************************************************
! テスト用メイン

!!PROGRAM main
!!USE	thpool_N2_lrand
!!IMPLICIT NONE
!!	INTEGER :: i
!!	REAL*8 :: dummy
!!
!!	dummy = pool_init(10.0d0)
!!	CALL pool_list
!!	
!!	DO i=0,100
!!		CALL pool_run(100)
!!		CALL pool_list
!!	END DO
!!
!!END PROGRAM main

!*** END of FILE *******************************************************