import numpy as np, matplotlib.pyplot as plt, konstanter as k, regning as regn, sympy

def read_data(filename):
    with open(filename,'r') as file:
        lines = np.array([line.split() for line in file.readlines()])
        lines = np.transpose(lines)

    data = np.array([None for line in lines])
    for i in range(len(lines)):
        data[i] = np.array([float(x) for x in lines[i][1:]])

    return data

##### leser ut data ######
data_matrix = read_data('data.txt')
V0_points = data_matrix[0]
T_points = data_matrix[1]
E_points = data_matrix[2]
p_points = data_matrix[3]

##### konverterer til SI (p i bar)####
V0_points = V0_points/1000
p_points = p_points*1333e-5
T_points = T_points + 273.15

########### regner ut E' og molalitet######
E_merket_points, s_E_merket_points, \
m_points, s_m_points, \
x_points, s_x_points = [[None for x in E_points] for i in range(6)]

for i in range(len(E_points)):
    E_merket_points[i], s_E_merket_points[i] = regn.E_merket(E=E_points[i], T = T_points[i], V0 = V0_points[i], p = p_points[i], usikkerhet=True)
    m_points[i], s_m_points[i] = regn.molalitet(V0=V0_points[i], usikkerhet = True)
    x_points[i], s_x_points[i] = regn.x_points(m=m_points[i], s_m = s_m_points[i], usikkerhet=True)

######### lineærregresjon #########
fit_x = [x for x in x_points]
fit_x.extend(fit_x)

fit_y = [y for y in E_merket_points]
fit_y.extend(fit_y)

fit_data = np.polyfit(fit_x,fit_y,1, cov=True)
coeff = fit_data[0]
cov = fit_data[1]

s_a = np.sqrt(cov[0][0])
s_b = np.sqrt(cov[1][1])

x_line = np.linspace(0,max(x_points),2)

y_line = coeff[0]*x_line + coeff[1]
upper = (coeff[0] + s_a)*x_line + (coeff[1] + s_b)
lower = (coeff[0] - s_a)*x_line + (coeff[1] - s_b)

########### finner E0, A og gamma #########
A, s_A = regn.A(T_points=T_points, a = coeff[0], s_a = s_a, usikkerhet=True)

gamma = [regn.midlere_aktkoeff(V0, usikkerhet=True) for V0 in V0_points]

for g,m in zip(gamma,m_points):
    print('g :', g[0],'±',g[1], '           ',  'm :', m)

print(s_m_points)
print()
print('E0 :', coeff[1], '±', s_b )
print('A :', A, '±', s_A )

########## plotting ############
jn = input('plot e_merket?')
if jn == 'j':
    plt.plot(x_line,y_line, color = 'black')
    plt.errorbar(x_points,E_merket_points, xerr=s_x_points, yerr=s_E_merket_points,
                 linestyle='', color='black', linewidth=2)
    plt.fill_between(x_line, lower, upper, color='black', alpha=0.3)

    plt.xlabel('m_HCl^(1/2) / [mol/kg]')
    plt.ylabel("E' / [V]")
    plt.savefig('fra_0_med_xy_feil_og_sky')
    plt.show()

######### plotter aktivitetskoeffisient ########
jn = input('plot aktkoeff?')
if jn == 'j':
    plt.errorbar(m_points, [g[0] for g in gamma], xerr=s_m_points, yerr=[g[1] for g in gamma],
                 linestyle='', color ='black', linewidth=2)
    plt.xlabel('m_HCl [mol/kg]')
    plt.ylabel('Midlere aktivitetskoeffisient [-]')

    plt.savefig('aktivitetskoeffisient')
    plt.show()