The planet method by python(physical computing lecture )

Use the odeint routine in the scipy.integrate library to integrate the motion of a planet around the sun.
the motion is planar and in the radial coordinate x= rcos$\theta$ , y=rsin$\theta$ , the equation reduce to a one - dimensional problem driven bt a potential U(r)

with the second Newton's law for r reading:

solve the three coupled differential equations for r, v_r and $\theta$ writing the rountine func(var,t)
that calcualates the three derivatives dr/dt , dv_r/dt and d$\theta$/dt that can be directly obtained from the equation above.

for instance, setting for simplicity L=1 , m1=1 , m2 = 1 , G= 1 and the boundary condition as : r0=1 m vr_0 = 0.9 , $\theta$_0 = 0.5

evloving the equation for a sufficiently long period of time you should get:

r0 = 1
vr0=0.9
theta0=0.5
var0=np.array([r0,vr0,theta0])
t=np.linspace(0,200,1000)

r,vr,theta=odeint(func,var0,t).T

pl.plot(rnp.cos(theta),rnp.sin(theta))

so what is the problem?

Changing the initial conditions one can get an ellipse or an open trajectory.
초기 조건을 바꿔가면서 타원형의 궤적인지 직선인지를 확인해보자
Can you parameterize the boundary conditions in order to predict the outcome?
경계 조건을 통해 다음의 결과를 예측할 수 있을까?

Hint: express the initial condition on vr in terms of the parameter x=v/v_esc, with
v=\sqrt{v_x^2+v_y^2} the total initial speed and v{esc} the escape velocity at r.
Verify that for x<1 the output is an ellipse, and for x>=1 it is an open trajectory.

우리가 물리학적으로 배운 탈출속력은 다음과 같다.
$ v_esc = \sqrt(2GM/r)$

만약 속력이 다음의 탈출속력보다 빠르다면 중력의 영향을 뚫고 앞으로 나아갈 수 있을 것이다.
때문에 이번 스푸트니크 1호가 발사함에 있어 다음의 escape 속력보다 빨라야지만 나아갈 수 있을 것이다.

그래서 이번 페이지에서는 우리가 원하는 속력보다 빠를 때를 예측하여 다음의 궤적을 예상해보자.

import numpy as np
from matplotlib import pyplot as pl
from scipy.integrate import odeint

L=1
G=1
m1=1
m2=1

def func(var,t):
r,vr,theta = var

dr_dt = vr
dvr_dt = L**2/(m1*m2*r**3) -G*m2/r**2
dtheta_dt = L/(m1*(r**2))

return [dr_dt, dvr_dt, dtheta_dt]

r0 = 1
x0 = 0.9
tmax = 100

v_esc = np.sqrt(2Gm2/r0)

#vr02 = 2Gm2/r0x02-L2/(m1r0**2)
vr0 = 0.9 # last point

theta0 = 0.5
var0 = np.array([r0,vr0,theta0])
t= np.linspace(0,200,1000)

" differential equation "
r,vr,theta = odeint(func,var0,t).T

pl.plot(rnp.cos(theta) , rnp.sin(theta))
pl.show()

x=vr0/v_esc
print(x)

if x<1 :
print('ellipse')
else :
print('line')