I am new to programming. Can anyone tell me what is RK ? how can i run this program? and what formula is used in this program? please explain with details, Thank you
import RK
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
“”" Implements d\vec X/dt = \vec V, d\vec V/dt = - M G \vec X/|X|^3.
Takes as input a time and
a 6-d vector. The first three components are the components of \vec
X, the next three are the components of \vec V. The parameter “t” is
not used.
“”"
def Kepler(t, y):
M = 1.0
G = 1.0
position = y[0:3]
velocity = y[3:6]
# position_dot = velocity
pos_dot = velocity;
r = np.sqrt(np.dot(position,position))
r3 = r*r*r
# formula for velocity_dot
vel_dot = - position * M * G / r3;
# y_dot
return np.concatenate((pos_dot, vel_dot))
“”" Source term for a binary. It has components m1 and m2. Currently
assumes circular orbits. Can adjust binary separation and mass
ratio.
“”"
def Binary(t, y):
binary_separation = 1.0
mass_ratio = 1.0
M = 1.0
G = 1.0
m1 = M * mass_ratio / (1.+mass_ratio)
m2 = M / (1.+mass_ratio)
binary_frequency = np.sqrt(M/binary_separation**3)
binary_angle = t * binary_frequency
X1 = np.array([np.cos(binary_angle), np.sign(binary_angle),0]) / (1.0 + mass_ratio)
X2 = -X1 * mass_ratio
position = y[0:3]
velocity = y[3:6]
# position_dot = velocity
pos_dot = velocity;
p1 = position - X1
p2 = position - X2
r1 = np.sqrt(np.dot(p1,p1))
r2 = np.sqrt(np.dot(p2,p2))
# formula for velocity_dot
vel_dot = - p1 * m1 * G / r1**3 - p2 * m2 * G / r2**3
# y_dot
return np.concatenate((pos_dot, vel_dot))
Initial position
x0 = 2
y0 = 0
z0 = 0
#Initial velocity
Vx = 0
Vy = np.sqrt(1.0/2.0)
Vz = 0
tend = 10000
dt = 0.01
y = np.array((x0,y0,z0, Vx, Vy, Vz), dtype=np.float64)
t = 0