!***********************************************************************
!	振動子と粒子の衝突
!	振動子を何としてもカオス的にする
!	MaxDeamon10a+.F90 をベースに修正、両端でエネルギー交換を行う
!	Date: 2007/12/15
!	g95 MaxDemon14_Spring1Exchg1.F90 mtfort90.o -o3 -Wno=165,112 -Wall -O -o MaxDemon14 -Wl,--subsystem,console -lGrWin -mwindows
PROGRAM MaxDemon14
USE mtmod	! MT乱数モジュール使用
IMPLICIT NONE

! GWライブラリについては以下の警告が出る
! Warning (165): Implicit interface 'gwopen' called at (1)
! Warning (112): Variable 'n' at (1) is set but never used
! 気になるのであれば、-Wno=165,112 として警告を抑制する

	! 引数
	CHARACTER*8 :: str
	INTEGER :: param1
	! 物体の数、1=等温壁、2=ボール、3=隔壁、4=ボール、5=等温壁
	INTEGER, PARAMETER :: N_BALL = 5
	! 連成振動子の数
	INTEGER, PARAMETER :: N_BODY = 16
	! 画面サイズ
	INTEGER, PARAMETER :: WIDTH = 500, HEIGHT = 300
	! 壁の位置
	REAL, PARAMETER :: W1POS = 0.0
	! 描画する物体の高さ、大きさ
	REAL, PARAMETER :: yy=200.0, rr=8.0
	! バネの自然長
	REAL, PARAMETER :: LEN0 = 50.0
	
	REAL, PARAMETER :: rx0 = 10.0, ry0=20.0	! 初期振動子の振幅
	REAL :: x(1:N_BALL)  = (/rx0, 125.0, 250.0, 375.0, real(WIDTH)+rx0/)
	REAL :: next_x(1:N_BALL)	! 1ステップ先の位置
	REAL :: dx(1:N_BALL) = (/0.0, 10.0, 0.0, 10.0, 0.0/)
	REAL :: m(1:N_BALL)  = (/1.0, 1.0, 5.0, 1.0, 1.0/)
	REAL :: t=0.0
!	REAL :: dt=0.1	! 計差精度粗くしてみる
	REAL :: dt=0.05 ! 計算精度標準
!	REAL :: dt=0.008	!高精度
	
	REAL, PARAMETER :: Kr  = 80.0	! 引力の係数
	REAL, PARAMETER :: Kin = 0.2	! 中心力の係数

	REAL, PARAMETER :: Eps = 0.05	! 引力の最小半径、0割を防ぐ
	REAL, PARAMETER :: Knr = (Kr / 3)	! 近接反発力の係数
	
	REAL, PARAMETER :: Gra = 8.0	! 重力
	
	! 連成振動子の位置と速度
	! (x,y)平面上の運動とする
	REAL :: rx_ini(1:N_BODY)
	REAL :: ry_ini(1:N_BODY)
	REAL :: rx(1:N_BODY)
	REAL :: ry(1:N_BODY)
	REAL :: rdx(1:N_BODY)
	REAL :: rdy(1:N_BODY)
	REAL :: next_rx(1:N_BODY)
	REAL :: next_ry(1:N_BODY)
	REAL :: next_rdx(1:N_BODY)
	REAL :: next_rdy(1:N_BODY)
	REAL, PARAMETER :: yy2 = 100.0
	REAL, PARAMETER :: LL0 = 50.0	! 振り子の自然長
	
	! ワープゾーンの位置
	INTEGER, PARAMETER :: DELTA_P = 50
	INTEGER :: P1FROM =  50, P1TO
	INTEGER, PARAMETER :: P2FROM = 225, P2TO = P2FROM+DELTA_P
	
	! ワープ用作業変数
	INTEGER, PARAMETER :: WARP_INTERVAL = 1500
		! 150 -- それでも気になるのでピコピコが生じないほど長くしてみた
		! 20 -- インターバル時間を変えても偏り傾向には影響ない
	INTEGER :: interval=0, warpcnt=0, totwarpcnt=0
	
	! 当たり判定フラグ
	LOGICAL :: colflag(1:N_BALL-1)
	LOGICAL :: acolflag
	
	! ヒストグラム
!!	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, PARAMETER :: MAXITER = 100

	! 交換したエネルギーの積算値
	REAL :: leftEG=0.0, rightEG=0.0

! BEGIN
	CALL getarg(1,str) ! 第 1 引数を得る
	param1 = strtoi( str )
	P1FROM = P1FROM + param1
	P1TO = P1FROM + DELTA_P
	PRINT "('P1FROM=',I3,' P1TO=',I3)", P1FROM, P1TO
	
#ifndef NOGR
	CALL GWopen(IR,0)
	CALL GWindow(IR, -400.0, -200.0, real(WIDTH + LEN0*3), real(HEIGHT+200))
			! 窓は少し大きめにとってある
	CALL GWsetogn(IR, 0, -NOG)
#endif
	
	! 連成振動子の位置と速度を初期化する
	CALL sgrnd(4357)	! ランダムシード
	DO i=1, N_BODY
		rx_ini(i) = (i-1) * LEN0
	!!	IF (i /= 1 .and. i /= N_BODY) THEN
	!!		rx_ini(i) = rx_ini(i) + real( grnd() * 10.0 - 5.0 )
	!!	END IF
		rx(i) = rx_ini(i)
		rdx(i) = 0.0
	END DO
!!	xの初速
	DO i=1, N_BODY, 2
		rx(i) = rx(i) + rx0
	END DO
	DO i=1+1, N_BODY, 2
		rx(i) = rx(i) - rx0
	END DO
	! yについても初期化
	DO i=1, N_BODY
		ry_ini(i) = yy2
		ry(i) = ry_ini(i) - LL0 + real( grnd() * 40.0 - 20.0 )
		rdy(i) = 0.0
	END DO
	
	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
		
			! 連成振動の運動を求める
			CALL Runge(t, rx, rdx, next_rx, next_rdx, &
&							ry, rdy, next_ry, next_rdy, dt, N_BODY)
			! 両端のrdxはそのまま等温壁の運動となる
		!! => この部分の相互作用は無し
		!!	next_x(1)       = next_rx(1)
		!!	next_x(N_BALL)  = next_rx(N_BODY) - (N_BODY-1)*LEN0 + WIDTH
			
			! １ステップ先のボールの位置を求める
			DO i=2, N_BALL-1
				next_x(i) = x(i) + dx(i) * dt
			END DO
			
			! 当たり判定
			acolflag = .FALSE.
		!!	DO i=1, N_BALL-1
		!!	i=1
			IF ( next_x(2) <= x(1) ) THEN
				colflag(1) = .TRUE.
				acolflag = .TRUE.
			ELSE
				colflag(1) = .FALSE.
			END IF
			DO i=2, N_BALL-2	!! 両端の等温壁との衝突は別に扱う
				colflag(i) = isCol( x, next_x, i, i+1 )
				if ( colflag(i) ) then
!					PRINT "('Hit! (',I2, I2,')',F6.2, F6.2, F6.2, F6.2)", &
!&						i, iter, x(i), dx(i), x(i+1), dx(i+1)
					acolflag = .TRUE.
				end if
			END DO
		!!	i=N_BALL-1
			IF ( next_x(N_BALL-1) >= x(N_BALL) ) THEN
				colflag(N_BALL-1) = .TRUE.
				acolflag = .TRUE.
			ELSE
				colflag(N_BALL-1) = .FALSE.
			END IF
			
			if ( .NOT. acolflag ) then	! 衝突していなければ
				exit	!ループ終了
			end if
			
			! ボールの衝突処理
		!!	DO i=1, N_BALL-1
		!!	i=1
			IF ( colflag(1) ) THEN
				dx(2) = - dx(2)		! ボールの向きを反転する
				! 等温壁との衝突、エネルギーを再分配する
				! 交換したエネルギーを積算する
				leftEG = leftEG + aveEG( dx(2), rdx(1) ) !!, 1 )
			!!	PRINT "('**** Letf Hit', F12.3, F12.3)", leftEG, rightEG
			END IF
			DO i=2, N_BALL-2
				if ( colflag(i) ) then
					CALL Collision( dx, m, i, i+1)
					! 等温壁との衝突処理、振動子との同期をとる
				!!	if (i .eq. 1) then
				!!		rdx(1) = dx(1)
				!!	else if (i .eq. N_BALL-1) then
				!!		rdx(N_BODY) = dx(N_BALL)
				!!	end if
				end if
			END DO
		!!	i=N_BALL-1
			IF ( colflag(N_BALL-1) ) THEN
				dx(N_BALL-1) = - dx(N_BALL-1)	! ボールの向きを反転する
				! 等温壁との衝突、エネルギーを再分配する
				! 交換したエネルギーを積算する
				rightEG = rightEG + aveEG( dx(N_BALL-1), rdx(N_BODY) ) !!, 2 )
			!!	PRINT "('**** Right Hit', F12.3, F12.3)", leftEG, rightEG
			END IF
			
		END DO
		if (iter >= MAXITER ) then
			! 両端で無限ループに陥ったときはボールを1ステップ戻して回避する
			if ( colflag(1) ) then
				PRINT *,'WARNING: 無限ループ回避(2)'
				next_x(2) = x(2) + abs(dx(2))
			end if
			if ( colflag(N_BALL-1) ) then
				PRINT *,'WARNING: 無限ループ回避(4)'
				next_x(N_BALL-1) = x(N_BALL-1) - abs(dx(N_BALL-1))
			end if
			if( colflag(2) .OR. colflag(3) ) then
				PRINT *,'ERROR: 衝突無限ループ'
				stop
			end if
		end if
		
		! 物体の移動
	!! => 両端は動かさない
	!!	DO i=1, N_BALL
		DO i=2, N_BALL-1
			x(i) = next_x(i)
		END DO
		
	!	! 順序が保たれているかどうかチェック
	!	IF ( x(2) > x(3) ) THEN
	!		PRINT *, 'ERROR: 順序違反2:3'
	!		x(2) = x(2) - dx(2)
	!	!!	x(3) = x(3) - dx(3)
	!	END IF
	!	IF ( x(3) > x(4) ) THEN
	!		PRINT *, 'ERROR: 順序違反3:4'
	!	!!	x(3) = x(3) - dx(3)
	!		x(4) = x(4) - dx(4)
	!	END IF
		
		! 振動子の位置と速度を確定する
		rx = next_rx
		ry = next_ry
		rdx = next_rdx
		rdy = next_rdy
		! 両端の等温壁の速度を確定する
	!! => 両端は動かさない
	!!	dx(1)      = next_rdx(1)
	!!	dx(N_BALL) = next_rdx(N_BODY)
		
		!ワープ処理 -- x(3)を２つのセグメント間で移動する
		IF ( interval == 0 ) THEN
			! 順方向ワープ
			IF( P2FROM < x(3) .and. x(3) < P2TO &
&				.and. (x(2) < P1FROM) &
&				.and. (P2TO < x(4)) ) THEN
				x(3) = x(3) - (P2FROM - P1FROM)
				warpcnt = warpcnt + 1
				totwarpcnt = totwarpcnt + 1
				PRINT "('+++ WARP +++', I4, I6, F12.3, F12.3)", warpcnt, totwarpcnt, leftEG, rightEG
			! 逆方向ワープ -- ここにelseを入れないと対称にならないぞ
			ELSE IF( P1FROM < x(3) .and. x(3) < P1TO &
&				.and. (x(2) < P1FROM) &
&				.and. (P2TO < x(4)) ) THEN
				x(3) = x(3) + (P2FROM - P1FROM)
				warpcnt = warpcnt - 1
				totwarpcnt = totwarpcnt + 1
				PRINT "('--- WARP ---', I4, I6, F12.3, F12.3)", warpcnt, totwarpcnt, leftEG, rightEG
			END IF
			! ワープ後しばらくの間は次のワープができない
			interval = WARP_INTERVAL
		ELSE
			interval = interval - 1
		END IF
		
#ifndef NOGR
		! 描画
		DO i=1, N_BALL, 2
			CALL GWsrect(IR, x(i)-rr, yy-30.0, x(i)+rr, yy+30.0, 5 ) ! 5=紫
		END DO
		
		! ボールを描く
		CALL GWsetpen(IR, 2, -1, -1, -1)	! 2=暗緑色
		CALL GWsetbrs(IR, 2, 1, -1)	!塗りつぶし
		DO i=2, N_BALL, 2
			CALL GWellipse(IR, x(i)-rr, yy-rr, x(i)+rr, yy+rr )
		END DO
		! 連成振動子を描く
		CALL GWsetpen(IR, 3, -1, -1, -1)	! 3=オリーブ色
		CALL GWsetbrs(IR, 3, 1, -1)	!塗りつぶし
		DO i=1, N_BODY
			CALL GWellipse(IR, rx(i)-rr, ry(i)-rr, rx(i)+rr, ry(i)+rr )
		END DO

		! 振り子の中心を描く
		CALL GWsetpen(IR, 4, -1, -1, -1)	! 4=紺色
		CALL GWsetbrs(IR, 4, 1, -1)	!塗りつぶし
		DO i=1, N_BODY
			CALL GWellipse(IR, rx_ini(i)-rr/2, ry_ini(i)-rr/2, rx_ini(i)+rr/2, ry_ini(i)+rr/2 )
		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
	!!	CALL printEG(rx, rdx, dx, m, dt)	! エネルギーの表示
		
		! ヒストグラムをとる
	!!	histcnt = histcnt + 1
	!!	IF ( histcnt >= 5000 ) THEN ! 5000回に１回サンプリング
	!!		histcnt = 0
	!!		! x(3)の位置を把握する
	!!		pos = int( x(3) / (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
		IF (totwarpcnt > 3000) THEN
			EXIT
		END IF
#endif
	END DO
	
#ifndef NOGR
	CALL GWSETBK(IR, 1)
	CALL GWREFRESH(IR)
	CALL GWQUIT(IR)
#endif

CONTAINS

!***********************************************************************
! 左右の等温壁でボールの反射とエネルギーの交換を行う
! 引数：ボールの速度、等温壁の速度
!! ３つめの引数は単なるデバッグ表示用
REAL FUNCTION aveEG( dx1, rdx1 ) !! , idx )
IMPLICIT NONE
	REAL, INTENT(INOUT) :: dx1, rdx1
!!	INTEGER, INTENT(IN) :: idx
	REAL :: ave, delta
	REAL :: ball_eg, ball_v
	REAL :: wall_eg, wall_v
	
	! 等温壁とエネルギーを再分配する。質量は共に等しく１
	ave = ( dx1**2 + rdx1**2 ) / 2
	delta = dx1**2 - ave
	
!!PRINT "('aveEG Zero   =', F12.5, F12.5, F12.5, F12.5)", dx1, rdx1, ave, delta
	
	delta = delta * 0.5		!交換係数
	
	ball_eg = dx1**2 - delta		! 負の値にはならないはず
	IF ( ball_eg < 0 ) THEN
		PRINT *, "ERROR: ball_eg < 0"
		aveEG = 0.0
		RETURN
	END IF
	ball_v = sqrt( ball_eg )	! |ボールの速さ| = Root(エネルギー)
	
	wall_eg = rdx1**2 + delta		! 負の値にはならないはず
	IF ( wall_eg < 0 ) THEN
		PRINT *, "ERROR: wall_eg < 0"
		aveEG = 0.0
		RETURN
	END IF
	wall_v = sqrt( wall_eg )	! |等温壁の速さ| = Root(エネルギー)
	
	IF ( dx1 > 0 ) THEN		! 方向に合わせてボールの速度を設定
		dx1 = ball_v
	ELSE
		dx1 = - ball_v
	END IF
	
	IF ( rdx1 > 0 ) THEN	! 方向に合わせて等温壁の速度を設定
		rdx1 = wall_v
	ELSE
		rdx1 = - wall_v
	END IF
	
!!PRINT "('aveEG After  =', I2, F12.5, F12.5, F12.5, F12.5)", idx, dx1, rdx1, ave, delta
	
	aveEG = delta	! 交換したエネルギーの大きさを返す
END FUNCTION aveEG

!***********************************************************************
! 衝突判定
! 返り値: 0=はずれ、1=あたり
LOGICAL FUNCTION isCol( x, next_x, i, j )
IMPLICIT NONE
	INTEGER, INTENT(IN) :: i, j
	REAL, INTENT(IN) :: x(:), next_x(:)
	REAL :: now_x1, now_x2, next_x1, next_x2
	
	now_x1 = x(i)
	now_x2 = x(j)
	next_x1 = next_x(i)
	next_x2 = next_x(j)
	
	! 通常判定
	! now_は等しくても衝突しない、こうしないと繰り返し衝突が起こる
	! next_は等しければ衝突する
	if ( (now_x1 .gt. now_x2 .and. next_x1 .le. next_x2) &
&	.or. (now_x1 .lt. now_x2 .and. next_x1 .ge. 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
	INTEGER, INTENT(IN) :: i, j
	REAL, INTENT(INOUT) :: dx(:)
	REAL, INTENT(IN) :: m(:)
	REAL :: dwork
	
!	PRINT "('Collision1 i=',I2,':',F7.2,' j=',I2,':',F7.2)", &
!&				i, dx(i), 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)
	return
END SUBROUTINE Collision

!**********************************************************************
! xの導関数
REAL FUNCTION fx(t, x, vx, y, vy, i, N)
IMPLICIT NONE
	INTEGER, INTENT(IN) :: i, N
	REAL, INTENT(IN) :: t, x(N), vx(N), y(N), vy(N)
	fx = vx(i)	! 定数=1.0
	return
END FUNCTION fx

!**********************************************************************
! yの導関数
REAL FUNCTION fy(t, x, vx, y, vy, i, N)
IMPLICIT NONE
	INTEGER, INTENT(IN) :: i, N
	REAL, INTENT(IN) :: t, x(N), vx(N), y(N), vy(N)
	fy = vy(i)	! 定数=1.0
	return
END FUNCTION fy

!**********************************************************************
! vxの導関数
REAL FUNCTION fvx(t, x, vx, y, vy, rxy, r0xy, i, N)
IMPLICIT NONE
	INTEGER, INTENT(IN) :: i, N
	REAL, INTENT(IN) ::  t, x(N), vx(N), y(N), vy(N), rxy(N,N), r0xy(N)
	REAL :: fx
	INTEGER :: j
!!	REAL, PARAMETER :: Kside = 1.0	! 両端のバネ定数
	
	! 粒子間の相互作用
	fx = 0.0
	DO j = MAX(1, i-1), MIN(i+1, N)	! 前後１個までと相互作用する
		IF ( j /= i .AND. rxy(i,j) > Eps ) THEN	! 0割を防ぐ
			fx = fx - (x(i)-x(j)) / rxy(i,j) + Knr * (x(i)-x(j)) / (rxy(i,j)**2)
		END IF
	END DO
	
	fvx = Kr * fx	! 引力の係数
	
	! 初期位置への中心力
!!	IF (i == 1 .or. i >= N) THEN	! 両端
!!		fvx = fvx + Kside * ( rx_ini(i) - x(i) )	! 係数値が異なる
!!	ELSE
		IF ( r0xy(i) > Eps ) THEN	! 0割を防ぐ
			fvx = fvx + Kin * (r0xy(i) - LL0) * (rx_ini(i) - x(i)) / r0xy(i)
		END IF
!!	END IF
	
	return
END FUNCTION fvx

!**********************************************************************
! vyの導関数
REAL FUNCTION fvy(t, x, vx, y, vy, rxy, r0xy, i, N)
IMPLICIT NONE
	INTEGER, INTENT(IN) :: i, N
	REAL, INTENT(IN) ::  t, x(N), vx(N), y(N), vy(N), rxy(N,N), r0xy(N)
	REAL :: fy
	INTEGER :: j
	
	! 粒子間の相互作用
	fy = 0.0
	DO j = MAX(1, i-1), MIN(i+1, N)	! 前後１個までと相互作用する
		IF ( j /= i .AND. rxy(i,j) > Eps ) THEN	! 0割を防ぐ
			fy = fy - (y(i)-y(j)) / rxy(i,j) + Knr * (y(i)-y(j)) / (rxy(i,j)**2)
		END IF
	END DO
	
	fvy = Kr * fy	! 引力の係数
	
	! 初期位置への求心力
	IF ( r0xy(i) > Eps ) THEN	! 0割を防ぐ
		fvy = fvy + Kin * (r0xy(i) - LL0) * (ry_ini(i) - y(i)) / r0xy(i)
	END IF
	
	! 重力で引っ張ってみる
	fvy = fvy - Gra
	
	return
END FUNCTION fvy

!**********************************************************************
! ルンゲ・クッタ法
SUBROUTINE Runge(t, x, vx, xout, vxout, &
					y, vy, yout, vyout, dt, N)
IMPLICIT NONE
	INTEGER, INTENT(IN) :: N
	REAL, INTENT(INOUT) :: t
	REAL, INTENT(IN) :: x(:), vx(:)
	REAL, INTENT(OUT) :: xout(:), vxout(:)
	REAL, INTENT(IN) :: y(:), vy(:)
	REAL, INTENT(OUT) :: yout(:), vyout(:)
	REAL, INTENT(IN) :: dt
	real :: h1
	real :: kx(4,N), kvx(4,N)
	real :: kx1(N), x1(N), x2(N), x3(N)
	real :: kvx1(N), vx1(N), vx2(N), vx3(N)
	real :: ky(4,N), kvy(4,N)
	real :: ky1(N), y1(N), y2(N), y3(N)
	real :: kvy1(N), vy1(N), vy2(N), vy3(N)
	real :: rxy(N,N)	! 粒子間相互の距離
	real :: r0xy(N)		! 振り子の基準点からの距離
	integer :: i, j
	
	h1 = dt/2.0
	
	! Step-1
	DO i=1,N
		DO j=i+1,N
			rxy(i,j) = ( (x(i)-x(j))**2 + (y(i)-y(j))**2 )
		!!	rxy(i,j) = rxy(i,j) * SQRT(rxy(i,j))	! 距離^3を算出
			rxy(j,i) = rxy(i,j)
		END DO
	END DO
	DO i=1,N
		r0xy(i) = ( (x(i)-rx_ini(i))**2 + (y(i)-ry_ini(i))**2 )
	END DO
	DO i=1,N
		kx(1,i)  =  fx(t, x, vx, y, vy, i, N) * dt
		kvx(1,i) = fvx(t, x, vx, y, vy, rxy, r0xy, i, N) * dt
		x1(i)  =  x(i) +  kx(1,i)/2.0
		vx1(i) = vx(i) + kvx(1,i)/2.0
	END DO
	DO i=1,N
		ky(1,i)  =  fy(t, x, vx, y, vy, i, N) * dt
		kvy(1,i) = fvy(t, x, vx, y, vy, rxy, r0xy, i, N) * dt
		y1(i)  =  y(i) +  ky(1,i)/2.0
		vy1(i) = vy(i) + kvy(1,i)/2.0
	END DO
	
	! Step-2
	DO i=1,N
		DO j=i+1,N
			rxy(i,j) = ( (x1(i)-x1(j))**2 + (y1(i)-y1(j))**2 )
		!!	rxy(i,j) = rxy(i,j) * SQRT(rxy(i,j))	! 距離^3を算出
			rxy(j,i) = rxy(i,j)
		END DO
	END DO
	DO i=1,N
		r0xy(i) = ( (x1(i)-rx_ini(i))**2 + (y1(i)-ry_ini(i))**2 )
	END DO
	DO i=1,N
		kx(2,i)  =  fx(t+h1, x1, vx1, y1, vy1, i, N) * dt
		kvx(2,i) = fvx(t+h1, x1, vx1, y1, vy1, rxy, r0xy, i, N) * dt
		x2(i)  =  x(i) +  kx(2,i)/2.0
		vx2(i) = vx(i) + kvx(2,i)/2.0
	END DO
	DO i=1,N
		ky(2,i)  =  fy(t+h1, x1, vx1, y1, vy1, i, N) * dt
		kvy(2,i) = fvy(t+h1, x1, vx1, y1, vy1, rxy, r0xy, i, N) * dt
		y2(i)  =  y(i) +  ky(2,i)/2.0
		vy2(i) = vy(i) + kvy(2,i)/2.0
	END DO
	
	! Step-3
	DO i=1,N
		DO j=i+1,N
			rxy(i,j) = ( (x2(i)-x2(j))**2 + (y2(i)-y2(j))**2 )
		!!	rxy(i,j) = rxy(i,j) * SQRT(rxy(i,j))	! 距離^3を算出
			rxy(j,i) = rxy(i,j)
		END DO
	END DO
	DO i=1,N
		r0xy(i) = ( (x2(i)-rx_ini(i))**2 + (y2(i)-ry_ini(i))**2 )
	END DO
	DO i=1,N
		kx(3,i)  =  fx(t+h1, x2, vx2, y2, vy2, i, N) * dt
		kvx(3,i) = fvx(t+h1, x2, vx2, y2, vy2, rxy, r0xy, i, N) * dt
		x3(i)  =  x(i) +  kx(3,i)
		vx3(i) = vx(i) + kvx(3,i)
	END DO
	DO i=1,N
		ky(3,i)  =  fy(t+h1, x2, vx2, y2, vy2, i, N) * dt
		kvy(3,i) = fvy(t+h1, x2, vx2, y2, vy2, rxy, r0xy, i, N) * dt
		y3(i)  =  y(i) +  ky(3,i)
		vy3(i) = vy(i) + kvy(3,i)
	END DO
	
	! Step-4
	DO i=1,N
		DO j=i+1,N
			rxy(i,j) = ( (x3(i)-x3(j))**2 + (y3(i)-y3(j))**2 )
		!!	rxy(i,j) = rxy(i,j) * SQRT(rxy(i,j))	! 距離^3を算出
			rxy(j,i) = rxy(i,j)
		END DO
	END DO
	DO i=1,N
		r0xy(i) = ( (x3(i)-rx_ini(i))**2 + (y3(i)-ry_ini(i))**2 )
	END DO
	DO i=1,N
		kx(4,i)  =  fx(t+dt, x3, vx3, y3, vy3, i, N) * dt
		kvx(4,i) = fvx(t+dt, x3, vx3, y3, vy3, rxy, r0xy, i, N) * dt
	END DO
	DO i=1,N
		ky(4,i)  =  fy(t+dt, x3, vx3, y3, vy3, i, N) * dt
		kvy(4,i) = fvy(t+dt, x3, vx3, y3, vy3, rxy, r0xy, i, N) * dt
	END DO
	
	t = t+dt
	
	DO i=1,N
		kx1(i)  =  (kx(1,i) + 2.0* kx(2,i) + 2.0* kx(3,i) +  kx(4,i)) / 6.0
		kvx1(i) = (kvx(1,i) + 2.0*kvx(2,i) + 2.0*kvx(3,i) + kvx(4,i)) / 6.0
	END DO
	DO i=1,N
		ky1(i)  =  (ky(1,i) + 2.0* ky(2,i) + 2.0* ky(3,i) +  ky(4,i)) / 6.0
		kvy1(i) = (kvy(1,i) + 2.0*kvy(2,i) + 2.0*kvy(3,i) + kvy(4,i)) / 6.0
	END DO

!	DO i=1,N
!		xout(i)  =  x(i)+ kx1(i)
!		vxout(i) = vx(i)+kvx1(i)
!	END DO
!	DO i=1,N
!		yout(i)  =  y(i)+ ky1(i)
!		vyout(i) = vy(i)+kvy1(i)
!	END DO
		xout  =  x +  kx1
		vxout = vx + kvx1
		yout  =  y +  ky1
		vyout = vy + kvy1

	RETURN
END SUBROUTINE Runge

!**********************************************************************
! 文字列(8)=>整数変換
INTEGER FUNCTION strtoi( num ) RESULT(i)
	character*8, intent(IN) :: num
	integer :: n, n1, k
	LOGICAL :: sgnflag = .false.
	
	k = len_trim(num)
	
	!先頭がマイナスだったら
	if (num(1:1) .eq. '-') then
		sgnflag = .true.
		n1 = 2
    else
		n1 = 1
	end if
	
	! 数字を１つずつ取り出して変換
	i = 0
	do n = n1, k
		i = i * 10 + ichar(num(n:n)) -48	! 48=ascii('0')
	end do
	
	!符号を整える
	if ( sgnflag ) then
		i = -i
	end if
END FUNCTION

END PROGRAM MaxDemon14
!*** END of FILE *******************************************************