# your code here

n = 0
a, b, t, p = 1, 1/(2**0.5), 1/4, 1
pi = (a+b)**2/4/t
err = 1000
history = []

while err >= 1e-15:
  print (f'n: {n}    pi: {pi:20.14f}')
  history.append([a,b,t,p,pi])  
  pi_old = pi
  a, b, t, p = (a+b)/2, (a*b)**0.5, t - p*(a-b)**2/4, 2*p
  pi = (a+b)**2/4/t
  err = abs(pi-pi_old)
  n += 1

# your code here

import numpy as np
import matplotlib.pyplot as plt

history = np.array(history)

plt.figure(figsize=(10,5), dpi=100)
plt.plot(history[:,4], label=r'$\pi_n$')
plt.legend()
plt.grid()
plt.xlabel('iteration')
plt.ylabel(r'$\pi_n$')
plt.show()

plt.figure(figsize=(10,5), dpi=100)
plt.plot(history[:,0], label=r'$a_n$')
plt.plot(history[:,1], label=r'$b_n$')
plt.legend()
plt.grid()
plt.xlabel('iteration')
plt.ylabel(r'$a_n$ and $b_n$')
plt.show()

# your code here

import numpy as np
import matplotlib.pyplot as plt

mu = 3.986e14
RE = 6400000

h0 = 800000
r0 = RE + h0
beta0 = np.pi/18

v0 = 9000

p = r0*v0*np.cos(beta0)
q = v0**2/2 - mu/r0
epsilon = np.sqrt( 1 + 2*p**2*q/mu**2 )

C = -np.arccos( p**2/r0/mu/epsilon - 1/epsilon )

print(f'p: {p}')
print(f'epsilon: {epsilon}')
print(f'C: {C}')

# your code here

theta = np.linspace(0, 2*np.pi, 1000)
r = p**2/mu/(1+epsilon*np.cos(theta-C))

plt.figure(figsize=(10,10), dpi=100)
plt.plot(0, 0, 'x', markersize=10, label='Earth\'s center')
plt.plot(6400*np.cos(theta), 6400*np.sin(theta), label='Earth\'s surface')
plt.plot(r*np.cos(theta)/1000, r*np.sin(theta)/1000, label='Satellite\'s trajectory')
plt.xlabel(r'$x$ position (km)')
plt.ylabel(r'$y$ position (km)')
plt.axis('equal')
plt.legend()
plt.grid()
plt.show()

# your code here

hmax = np.max(r) - RE
hmin = np.min(r) - RE

print(f'maximum h: {hmax/1000}(km)')
print(f'minimum h: {hmin/1000}(km)')

# your code here

import numpy as np
import matplotlib.pyplot as plt

v0_test = np.arange(9000, 0, -10)

for v0 in v0_test:

  p = r0*v0*np.cos(beta0)
  q = v0**2/2 - mu/r0
  epsilon = np.sqrt( 1 + 2*p**2*q/mu**2 )
  C = -np.arccos( p**2/r0/mu/epsilon - 1/epsilon )
  theta = np.linspace(0, 2*np.pi, 1000)
  r = p**2/mu/(1+epsilon*np.cos(theta-C))

  if np.min(r) <= 6400e3:
    break

print(f'minimum v_0: {v0} (m/s)')
print(f'minimum r: {np.min(r)/1000} (km)')

plt.figure(figsize=(10,10), dpi=100)
plt.plot(0, 0, 'x', markersize=10, label='Earth\'s center')
plt.plot(6400*np.cos(theta), 6400*np.sin(theta), label='Earth\'s surface')
plt.plot(r*np.cos(theta)/1000, r*np.sin(theta)/1000, label='Satellite\'s trajectory')
plt.xlabel(r'$x$ position (km)')
plt.ylabel(r'$y$ position (km)')
plt.axis('equal')
plt.legend()
plt.grid()
plt.show()