import numpy as np
import matplotlib.pyplot as plt

SQ = np.loadtxt('E:/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()
t0 = [0]
for i in range(len(MJD)):
    if Int[i]>8.4 and abs(MJD[i]-t0[-1])>10:
        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_vector = 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_vector = f0
    for j in range(len(t0)):
        X[:,1] += np.exp(-np.power(MJD-t0[j],2)/2/s**2)
        X[:,2] += A/s**3*np.multiply(np.exp(-np.power(MJD-t0[j],2)/2/s**2),np.power(MJD-t0[j],2))
        X[:,3+j] = A/s**2*np.multiply(np.exp(-np.power(MJD-t0[j],2)/2/s**2),MJD-t0[j])
        f += A*np.exp(-np.power(MJD-t0[j],2)/2/s**2)
        f_vector += A*np.exp(-np.power(t_vector-t0[j],2)/2/s**2)
    plt.plot(MJD,Int,'.',t_vector,f_vector,'r')
    plt.gca().invert_yaxis()
    #plt.xlim([52320,52350])
    plt.show()
    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