!***********************************************************************
!	振動子と粒子の衝突
!	振動子を何としてもカオス的にする
!	MaxDeqamon14_Spring1.F90 を高精度にし、バネ粒子系を軽くした。
!	Date: 2007/11/24
!	g95 MaxDemon14_Spring1Light.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*8, PARAMETER :: W1POS = 0.0d0
	! 描画する物体の高さ、大きさ
	REAL*8, PARAMETER :: yy=200.0d0, rr=8.0d0
	! バネの自然長
	REAL, PARAMETER :: LEN0 = 50.0d0
	
	REAL*8, PARAMETER :: rx0 = 10.0d0, ry0=20.0d0	! 初期振動子の振幅
	REAL*8 :: x(1:N_BALL)  = (/rx0, 125.0d0, 250.0d0, 375.0d0, real(WIDTH)+rx0/)
	REAL*8 :: next_x(1:N_BALL)	! 1ステップ先の位置
	REAL*8 :: dx(1:N_BALL) = (/0.0d0, 10.0d0, 0.0d0, 10.0d0, 0.0d0/)
	REAL*8 :: m(1:N_BALL)  = (/1.0d0, 1.0d0, 5.0d0, 1.0d0, 1.0d0/)
	REAL*8 :: t=0.0d0
!	REAL*8 :: dt=0.1d0	! 計差精度粗くしてみる
	REAL*8 :: dt=0.05d0 ! 計算精度標準
!	REAL*8 :: dt=0.008d0	!高精度
	
	REAL*8, PARAMETER :: Kr  = 8.0d0	! 引力の係数
	REAL*8, PARAMETER :: Kin = 0.02d0	! 中心力の係数

	REAL*8, PARAMETER :: Eps = 0.05d0	! 引力の最小半径、0割を防ぐ
	REAL*8, PARAMETER :: Knr = (Kr / 3)	! 近接反発力の係数
	
	REAL*8, PARAMETER :: Gra = 0.8d0	! 重力
	
	! 連成振動子の位置と速度
	! (x,y)平面上の運動とする
	REAL*8 :: rx_ini(1:N_BODY)
	REAL*8 :: ry_ini(1:N_BODY)
	REAL*8 :: rx(1:N_BODY)
	REAL*8 :: ry(1:N_BODY)
	REAL*8 :: rdx(1:N_BODY)
	REAL*8 :: rdy(1:N_BODY)
	REAL*8 :: next_rx(1:N_BODY)
	REAL*8 :: next_ry(1:N_BODY)
	REAL*8 :: next_rdx(1:N_BODY)
	REAL*8 :: next_rdy(1:N_BODY)
	REAL*8, PARAMETER :: yy2 = 100.0d0
	REAL*8, PARAMETER :: LL0 = 50.0d0	! 振り子の自然長
	
	! ワープゾーンの位置
	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

! 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, 10
!!		rx(i) = rx(i) + rx0
!!	END DO
!!	DO i=5+1, N_BODY, 10
!!		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
				colflag(i) = isCol( x, next_x, 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
			
			! ボールの衝突処理
			DO i=1, N_BALL-1
				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
			
		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
			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)", warpcnt, totwarpcnt
			! 逆方向ワープ -- ここに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)", warpcnt, totwarpcnt
			END IF
			! ワープ後しばらくの間は次のワープができない
			interval = WARP_INTERVAL
		ELSE
			interval = interval - 1
		END IF
		
#ifndef NOGR
		! 描画
		DO i=1, N_BALL, 2
			CALL GWsrect(IR, real(x(i)-rr), real(yy-30.0), real(x(i)+rr), real(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, real(x(i)-rr), real(yy-rr), real(x(i)+rr), real(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, real(rx(i)-rr), real(ry(i)-rr), real(rx(i)+rr), real(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, real(rx_ini(i)-rr/2), real(ry_ini(i)-rr/2), real(rx_ini(i)+rr/2), real(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

!***********************************************************************
! エネルギーの表示
! 両端は -LEN0, WIDTH+LEN0 で固定している
SUBROUTINE printEG(rx, rdx, dx, m, dt)
IMPLICIT NONE
	REAL*8, INTENT(IN) :: rx(:), rdx(:), dx(:), m(:), dt
	REAL*8 :: eg0v, eg1v, eg0p, eg1p, eg2, eg3, eg4
	INTEGER :: i, prev_i
	REAL*8 :: prev_x, next_x

	! 連成振動子の運動エネルギー
	eg1v = 0.0
	DO i=1,N_BODY
		eg0v = dt * rdx(i) ** 2 / 2.0		!質量は1.0
!		PRINT "('v',I2,': ',F11.3,': ',F11.3)", i, rdx(i), eg0v
		eg1v = eg1v + eg0v
	END DO
	
	! 連成振動子のバネに蓄積されているエネルギー
	eg1p = 0.0
	DO i=1,N_BODY+1
		prev_i = i-1
		if (prev_i .le. 0) then		! 左端
			prev_x = - LEN0
		else
			prev_x = rx(prev_i)
		end if
		if (i .gt. N_BODY) then		! 右端
			next_x = WIDTH + LEN0
		else
			next_x = rx(i)
		end if
		! 自然長からのずれを数える
		eg0p = (LEN0 - (next_x - prev_x)) ** 2 / 2.0 !質量,バネ定数は1.0
!		PRINT "('p',I2,': ',F8.3,': ',F8.3,': ',F11.3)",&
!&			i, prev_x, next_x, eg0p
		eg1p = eg1p + eg0p
	END DO
	
	eg2  = dt * dx(2) ** 2 * m(2) / 2.0
	eg3  = dt * dx(3) ** 2 * m(3) / 2.0
	eg4  = dt * dx(4) ** 2 * m(4) / 2.0
	
	PRINT "(F6.3,': ',F13.3,': ',F13.3,': ',F8.3,': ',F6.3,': ',F6.3,': ',F13.3)",&
			eg1v, eg1p, eg1v+eg1p, eg2, eg3, eg4, &
			eg1v+eg1p+eg2+eg3+eg4

	return
END SUBROUTINE printEG

!***********************************************************************
! 衝突判定
! 返り値: 0=はずれ、1=あたり
LOGICAL FUNCTION isCol( x, next_x, i, j )
IMPLICIT NONE
	INTEGER, INTENT(IN) :: i, j
	REAL*8, INTENT(IN) :: x(:), next_x(:)
	REAL*8 :: 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*8, INTENT(INOUT) :: dx(:)
	REAL*8, INTENT(IN) :: m(:)
	REAL*8 :: 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*8 FUNCTION fx(t, x, vx, y, vy, i, N)
IMPLICIT NONE
	INTEGER, INTENT(IN) :: i, N
	REAL*8, INTENT(IN) :: t, x(N), vx(N), y(N), vy(N)
	fx = vx(i)	! 定数=1.0
	return
END FUNCTION fx

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

!**********************************************************************
! vxの導関数
REAL*8 FUNCTION fvx(t, x, vx, y, vy, rxy, r0xy, i, N)
IMPLICIT NONE
	INTEGER, INTENT(IN) :: i, N
	REAL*8, INTENT(IN) ::  t, x(N), vx(N), y(N), vy(N), rxy(N,N), r0xy(N)
	REAL*8 :: fx
	INTEGER :: j
!!	REAL*8, PARAMETER :: Kside = 1.0	! 両端のバネ定数
	
	! 粒子間の相互作用
	fx = 0.0d0
	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*8 FUNCTION fvy(t, x, vx, y, vy, rxy, r0xy, i, N)
IMPLICIT NONE
	INTEGER, INTENT(IN) :: i, N
	REAL*8, INTENT(IN) ::  t, x(N), vx(N), y(N), vy(N), rxy(N,N), r0xy(N)
	REAL*8 :: fy
	INTEGER :: j
	
	! 粒子間の相互作用
	fy = 0.0d0
	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*8, INTENT(INOUT) :: t
	REAL*8, INTENT(IN) :: x(:), vx(:)
	REAL*8, INTENT(OUT) :: xout(:), vxout(:)
	REAL*8, INTENT(IN) :: y(:), vy(:)
	REAL*8, INTENT(OUT) :: yout(:), vyout(:)
	REAL*8, INTENT(IN) :: dt
	real*8 :: h1
	real*8 :: kx(4,N), kvx(4,N)
	real*8 :: kx1(N), x1(N), x2(N), x3(N)
	real*8 :: kvx1(N), vx1(N), vx2(N), vx3(N)
	real*8 :: ky(4,N), kvy(4,N)
	real*8 :: ky1(N), y1(N), y2(N), y3(N)
	real*8 :: kvy1(N), vy1(N), vy2(N), vy3(N)
	real*8 :: rxy(N,N)	! 粒子間相互の距離
	real*8 :: 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 *******************************************************