import numpy as np
import matplotlib.pyplot as plt

D = np.arange(366)
M = 0.9856*D - 3.14   # stredna anomalia Zeme
# Ekliptikalna dlzka Zeme.
lambda_ = 283.14 + M + 1.91*np.sin(M*np.pi/180) + 0.02*np.sin(2*M*np.pi/180)
# Ekliptikalna sirka Zeme.
delta = np.arcsin(np.sin(23.433*np.pi/180)*np.sin(lambda_*np.pi/180))
t = D/365.25*2*np.pi   # stredny slnecny cas
# Rozdiel medzi pravym a strednym slnecnym casom.
dt = 0.46*np.cos(t)-7.34*np.sin(t)-3.36*np.cos(2*t)-9.33*np.sin(2*t);

plt.plot(dt, delta*180/np.pi, '.')
plt.title('Analema')
plt.xlabel('dt [min]')
plt.ylabel('delta [°]')