import numpy as np
import matplotlib.pyplot as plt

matrix = np.loadtxt('E:\cont0109a.txt')
x = matrix[:,0]   #vlnova dlzka
y = matrix[:,1]   #intenzita
n = len(x)

b0 = [1.0,-0.6,4340.0,5.0]
for i in range(5):
    G0 = b0[0]   #intenzita kontinua
    A = b0[1]   #amplituda
    x0 = b0[2]   #stred ciary
    sigma = b0[3]   #polosirka ciary
    G = G0+A*np.exp(-np.power(x-x0,2)/2./sigma**2)   #prekladana Gaussova funkcia
    X = np.ones((n,4))
    X[:,1] = np.exp(-np.power(x-x0,2)/2./sigma**2)   #derivacia podla A
    X[:,2] = A*np.exp(-np.power(x-x0,2)/2./sigma**2)*(x-x0)/sigma**2   #derivacia podla x0
    X[:,3] = A*np.exp(-np.power(x-x0,2)/2./sigma**2)*np.power(x-x0,2)/sigma**3   #derivacia podla sigma
    dY=y-G
    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
    plt.plot(x,y,'b')
    plt.plot(x,G,'r')
    plt.show()
    chi2_G = np.dot(dY.T,dY)

b0 = [1.0,-3.0,4340.0,5.0]
for i in range(5):
    L0 = b0[0]     #intenzita kontinua
    A = b0[1]      #amplituda
    x0 = b0[2]     #stred ciary
    gamma = b0[3];
    L = L0+A*gamma/(np.power(x-x0,2)+gamma**2)    #prekladana Lorentzova funkcia
    X = np.ones((n,4))
    X[:,1] = gamma/(np.power(x-x0,2)+gamma**2)   #derivacia podla A
    X[:,2] = 2*A*gamma*(x-x0)/np.power(np.power(x-x0,2)+gamma**2,2)   #derivacia podla x0
    X[:,3] = A*(np.power(x-x0,2)-gamma**2)/np.power(np.power(x-x0,2)+gamma**2,2)   #derivacia podla gamma
    dY=y-L
    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
    plt.plot(x,y,'b')
    plt.plot(x,L,'g')
    plt.show()
    chi2_L = np.dot(dY.T,dY)

plt.plot(x,y,'b')
plt.plot(x,L,'g')
plt.plot(x,G,'r')
plt.show()
print('chi2 pre Gaussovu funkciu je ',chi2_G)
print('chi2 pre Lorentyovu funkciu je ',chi2_L)