program planet_rk4
	implicit none
	real(8) :: x = 1.0
	real(8) :: y = 0.0
	real(8) :: vx = 0.0
	real(8) :: vy = 1.0
	real(8) :: h = 0.1
	real(8) :: t = 0.0
	real(8) :: A1,A2,A3,A4
	real(8) :: B1,B2,B3,B4
	real(8) :: C1,C2,C3,C4
	real(8) :: D1,D2,D3,D4
	character(256) :: argv
	
	call getarg(1,argv)
	open(9,file=argv)
	
	do while(t <= 50)
		write(9,*) x,y
		
		A1 = h * f(t,vx,x,y)
		B1 = h * g(t,vx,x)
		
		A2 = h * f(t+h/2.0,vx+A1/2.0,x+B1/2.0,y)
		B2 = h * g(t+h/2.0,vx+A1/2.0,x+B1/2.0)
		
		A3 = h * f(t+h/2.0,vx+A2/2.0,x+B2/2.0,y)
		B3 = h * g(t+h/2.0,vx+A2/2.0,x+B2/2.0)
		
		A4 = h * f(t+h,vx+A3,x+B3,y)
		B4 = h * g(t+h,vx+A3,x+B3)
		
		vx = vx + 1.0/6.0 * (A1 + 2.0*A2 + 2.0*A3 + A4)
		x = x + 1.0/6.0 * (B1 + 2.0*B2 + 2.0*B3 + B4)
		
		
		C1 = h * f(t,vy,y,x)
		D1 = h * g(t,vy,y)
		
		C2 = h * f(t+h/2.0,vy+C1/2.0,y+D1/2.0,x)
		D2 = h * g(t+h/2.0,vy+C1/2.0,y+D1/2.0)
		
		C3 = h * f(t+h/2.0,vy+C2/2.0,y+D2/2.0,x)
		D3 = h * g(t+h/2.0,vy+C2/2.0,y+D2/2.0)
				
		C4 = h * f(t+h,vy+C3,y+D3,x)
		D4 = h * g(t+h,vy+C3,y+D3)
		
		vy = vy + 1.0/6.0 * (C1 + 2.0*C2 + 2.0*C3 + C4)
		y = y + 1.0/6.0 * (D1 + 2.0*D2 + 2.0*D3 + D4)
		
		t = t + h
	end do
	
	close(9)
	stop
	
	
	contains
	real(8) function f(t,v,x,y)
		implicit none
		real(8),intent(in) :: t,v,x,y
		real(8) :: r2,r3
		
		r2 = x*x + y*y
		r3 = r2 * Sqrt(r2)
		f = -x/r3
		
	end function f
	
	real(8) function g(t,v,x)
		implicit none
		real(8),intent(in) :: t,v,x
		
		g = v
		
	end function g
end program planet_rk4