import numpy as np
import matplotlib.pyplot as plt

SQ = np.loadtxt('SQ.dat')
MJD = SQ[:, 0]
Int = SQ[:, 1]
plt.plot(MJD, Int, '.')
plt.gca().invert_yaxis()
#plt.xlim([52320,52350])   #detail konkretneho zakrytu
plt.show()
plt.figure()
t0 = [0]
for i in range(len(MJD)):
    if Int[i] > 8.4 and abs(MJD[i] - t0[-1]) > 100:
        t0.append(MJD[i])
#pociatocne odhady pozicii minim
del t0[0]
b0 = [8.15, 0.3, 0.7]
b0.extend(t0)   #pociatocne odhady vstupnych parametrov
t_fit = np.arange(min(MJD), max(MJD), 0.1)
for i in range(5):
    f0 = b0[0]   #nulova hladina
    A = b0[1]   #amplituda
    s = b0[2]   #polosirka
    t0 = b0[3:]   #vektor okamihov minim
    X = np.zeros((len(MJD), len(b0)))
    X[:, 0] = 1
    f = f0
    f_fit = f0
    for j in range(len(t0)):
        dt2 = np.power(MJD - t0[j], 2)
        X[:, 1] += np.exp(-dt2/2/s**2)
        X[:, 2] += A/s**3*np.multiply(np.exp(-dt2/2/s**2), dt2)
        X[:, 3+j] = A/s**2*np.multiply(np.exp(-dt2/2/s**2), MJD - t0[j])
        f += A*np.exp(-dt2/2/s**2)
        f_fit += A*np.exp(-np.power(t_fit - t0[j], 2)/2/s**2)
    plt.clf()
    plt.plot(MJD, Int, '.', t_fit, f_fit, 'r')
    plt.gca().invert_yaxis()
    #plt.xlim([52320,52350])
    plt.show()
    plt.pause(1)
    dY = Int - f;
    V = np.dot(np.transpose(X), X)
    U = np.dot(np.transpose(X), dY)
    db=np.dot(np.linalg.inv(V), U)
    b0 = db + b0