Kepler's Laws (2021.01)
===

###### tags: `Side Project`, `Python`

## Brief

## Result


## Code

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm
from matplotlib import style
import math

style.use('ggplot')

epsilon = 0.0001 #時間間隔
n = 500000 #迭代次數
position = np.zeros((n, 2)) #紀錄每個時間點的距離
velocity = np.zeros((2)) #紀錄當下的速度
acceleration = np.zeros((2)) #紀錄當下的加速度

y = [1, 1.05, 1.1, 1.15, 1.2, 1.25, 1.3]

# y = [1]

rx = 0
ry = 0
lx = 0
lx = 0
uy = 0
dy = 0
kkk = []

for i in range(len(y)):

	position[0, 0] = 1
	position[0, 1] = 0
	velocity[0] = 0
	velocity[1] = y[i]
	y0 = y[i]

	r = norm(position[0]) #計算離太陽距離
	r3 = 1 / (r**3) #計算 1/r^3
	acceleration = position[0] * -r3 #計算星球加速度
	velocity += acceleration * epsilon * 0.5 #計算星球速度	

	i = 1

	while(i < n): #重複n次
	    position[i] = position[i-1] + velocity * epsilon #計算更新後的星球位置

	    if position[i, 0] > rx:
	    	rx = position[i, 0]
	    if position[i, 0] < lx:
	    	lx = position[i, 0]
	    	ly = position[i, 1]
	    if position[i, 1] > uy:
	    	ux = position[i, 0]
	    	uy = position[i, 1]
	    if position[i, 1] < dy:
	    	dx = position[i, 0]
	    	dy = position[i, 1]

	    r = norm(position[i]) #計算恆星與星球間的距離
	    r3 = 1 / (r**3) #計算 1/r^3
	    acceleration = position[i] * -r3 #計算星球加速度
	    velocity += acceleration * epsilon #計算星球速度
	    i += 1

	name = str(y0)
	data, = plt.plot(position[:, 0], position[:, 1], label=name)
	kkk.append(data)
	plt.plot(0, 0, "ro") #sun
	# 座標數字
	plt.text(rx, ry, "("+str(round(rx, 3))+","+str(round(ry, 3))+")") 
	plt.text(lx, ly, "("+str(round(lx, 3))+","+str(round(ly, 3))+")") 
	plt.text(ux, uy, "("+str(round(ux, 3))+","+str(round(uy, 3))+")") 
	plt.text(dx, dy, "("+str(round(dx, 3))+","+str(round(dy, 3))+")") 


	# 計算a,b,離心率
	print("y = "+str(y0))
	p1 = np.array([rx,ry])
	p2 = np.array([lx,ly])
	p3 = p2-p1
	a = (math.hypot(p3[0],p3[1]))/2
	print("a = "+str(a))
	p1 = np.array([ux,uy])
	p2 = np.array([dx,dy])
	p3 = p2-p1
	b = (math.hypot(p3[0],p3[1]))/2
	print("b = "+str(b))
	if (y0) == 1:
		Eccentricity = 0
	else:
		Eccentricity = math.sqrt(1-((b*b)/(a*a)))
	# e=(根号(a^2-b^2))/a=根号(1-(b/a)^2)
	print("Eccentricity = "+str(Eccentricity))

	# 計算週期T
	E = ((y0**2)/2)-1
	p = -2*E
	p = p**3
	p = math.sqrt(p)
	p = 1/p
	PI = math.pi
	T = 2*PI*p
	print("T = "+str(T))

	#克卜勒第三
	# K = (T**2)/(((a+b)/2)**3)
	K = (T**2)/(a**3)
	print("Kepler's laws = "+str(K))
	print("*****************************************")

plt.legend(handles=kkk, loc = 'upper left')
plt.show()