import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt('star.txt')
HJD = data[:, 0]
Y = data[:, 1]
Sigma = np.array([])
Period = np.array([])
# Potrebujeme velmi jemne delenie pre periodu, aby sme ju nahodou neminuli.
d = 0.00002
for P in np.arange(0.9, 0.95, d):
    fi = np.modf(HJD/P)[0]
    X = np.ones((len(fi), 5))
    X[:, 1] = np.cos(2*np.pi*fi)
    X[:, 2] = np.sin(2*np.pi*fi)
    X[:, 3] = np.cos(4*np.pi*fi)
    X[:, 4] = np.sin(4*np.pi*fi)
    V = np.dot(X.T, X)
    U = np.dot(X.T, Y)
    b = np.dot(np.linalg.inv(V), U)
    Yp = np.dot(X, b)
    DY = Y - Yp
    s = np.sqrt(np.dot(DY.T, DY)/(len(X) - len(X[0])))
    Sigma = np.append(Sigma, s)
    Period = np.append(Period, P)
    # Pre nazornost si vykreslime svetelnu krivku pre periody v okoli spravnej
    # hodnoty a sledujeme ako sa zvacsuje presnost fitu. Tu trochu podvadzame, 
    # pretoze dopredu vieme, ze spravna perioda je 0.9366 dna.
    if abs(P - 0.9366) < d*50:
        plt.clf()
        plt.plot(fi, Y, 'bo', fi, Yp, '.r')
        plt.gca().invert_yaxis()
        plt.xlabel('faza')
        plt.ylabel('zdanliva hviezdna velkost [mag]')
        plt.title('P = {} d'.format(np.round(P, 5)))
        plt.show()
        if abs(P - 0.9366) < d/2:
            plt.pause(2)
        else:
            plt.pause(0.2)
    
P = Period[np.where(Sigma == min(Sigma))]
plt.figure()
plt.plot(Period, Sigma)
plt.xlabel('Period [d]')
plt.ylabel('standard deviation')
plt.title('P = {} d'.format(np.round(P[0], 4)))