'''
    Solve Kepler problem at N=100 points along the orbit 
    of a Saturn-class exoplanet K2-261b discovered in 2018 around a
    sun-like star K2-261 
        see exoplanet.eu/catalog/k2-261_b/
    at uniformly distributed times numbered by 
        t_i = (i/100)*P, 
    where P = orbital period, and index i = 0...99.
    
    Find 100 eccentric anomalies E(t_i) for 
    e = eccentricity = 0.39
    
    In this preliminary version of the code, we are only solving by 
    the following non-algebraic equation.
    Kepler equation
        E - e sin(E) = 2 pi t/P
    cannot be solved on paper, but can be solved iteratively in Python. 
    Here we have 
    M  = 2 pi t/P = the right-hand side of Kepler equation, 
        a.k.a. mean anomaly (angle),
    E = eccentric anomaly (angle),
    (The reason for these names is astronomical tradition.)
    
    This allows us to start with the phase of the motion (=fraction t/P, 
    or time t in units of full orbital periods), obtain E by
    iteration, and down the line calculate theta, r(t), and the Cartesian corrdinates (x,y). But that's in another version of this code. 
        
    Summary of workflow:
    t --> M(t) --> E(t) [--> Theta(t) --> r(t) --> (x,y), v^2, dE/dt, etc.]
'''

from numpy import pi,sqrt,tan,arctan,sin,cos
import matplotlib.pyplot as plt
#   input number of points N
N = 100
# eccentricity and semi-maj. axis of the planet
e = 0.39

#   prepare plot frame and description
plt.figure(figsize=(8,6.5))
plt.xlabel('t ',fontsize=14)
plt.ylabel('E ',fontsize=14)
plt.title(" Eccentric anomaly from Kepler's equation;  e=0.39",
          fontsize=14)
plt.grid(True)

# compute and plot in x,y coordinates N points of the orbit
for i in range(N) :
    M = 2*pi * i/N   # mean anomaly growing linearly from 0 to almost 2pi
    E = M            # first guess to start with 
# Iterative solution for one t    
    for k in range(99):  # counter of iterations
        darken = 1./(1.+0.5*k)
        plt.scatter(M, E, s=4, color=(0,.7*darken,.8*darken), linewidth=2)
        E_previous = E   # save previous E for comparison
        E = M + e*sin(E) # Kepler equation iteration
        if (i%10==0):   
            print('i',i,' k',k," E =",E)  
        if (abs(E-E_previous)< 1e-9): 
            break  # required accuracy reached
#   plot one point on E(M) or E(t) curve
    plt.scatter(M,E,s=7,color=(1.,.0,.0),linewidth=2)
plt.show()