import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt('hd22049.dat')
conv = 10**(-10)
lam = data[:, 0]*conv   # vlnova dlzka v m
F = data[:, 1]   # ziarivy tok
plt.plot(lam/conv, F)   # vlnovu dlzku prevedieme na Angstromy
plt.xlabel('vlnova dlzka [A]')
plt.ylabel('Intenzita')

# a)
bins = np.arange(min(lam), max(lam), 100*conv)
x = np.array([])
y = np.array([])
for i in range(1, len(bins)):
    interval = np.where((lam > bins[i - 1]) & (lam <= bins[i]))
    if interval:
        y_max = max(F[interval])
        y = np.append(y, y_max)
        x_max = lam[np.where(F == y_max)][0]
        x = np.append(x, x_max)

plt.plot(x/conv, y,'*')

# b)
# Z Wienovho zakona lam_max*T = b urcime priblizne teplotu T = 5000 K.
# Dosadenim ziskame pribliznu hodnotu pre B = 3E-6. Hodnotu pre A dostaneme
# dopocitanim z funkcnej hodnoty: A = F_max*lam_max**5*(exp(B/lam_max) - 1).
B = 3*10**(-6)
A = 10**(-25)
for i in range(5):
    f = A/x**5/(np.exp(B/x) - 1)
    dy = y - f
    X = np.zeros((len(x), 2))
    X[:, 0] = f/A   # derivacia podla A
    X[:, 1] = -A*np.exp(B/x)/x**6/(np.exp(B/x) - 1)**2   # derivacia podla B
    V = np.dot(X.T, X)
    U = np.dot(X.T, dy)
    db = np.dot(np.linalg.inv(V), U)
    A += db[0]
    B += db[1]
    plt.pause(1)
    plt.clf()
    plt.plot(lam/conv, F)
    plt.plot(x/conv, y,'*')
    plt.plot(x/conv, f)
    plt.xlabel('vlnova dlzka [A]')
    plt.ylabel('Intenzita')
    plt.show()

# c)
h = 6.636*10**(-34)
c = 299792458
k_B = 1.3806*10**(-23)
print('Hodnota parametra B je {} a z nej dostaneme teplotu T = {} K.'
      .format(np.round(B, 9), np.round(h*c/k_B/B)))