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

def data_read(filename):
    with open(filename, 'r') as file:
        lines = [line.split()[1:] for line in file.readlines()]
    data = [[float(x) for x in line] for line in lines]
    return np.array(data)

def gen_plot():
    plt.errorbar(x_points, y_points, xerr=s_x_points, yerr=s_y_points,
                 color='black', linestyle='', linewidth=0.5, capsize=2)
    plt.plot(x_line, y_line, color='black', linewidth = 0.5)
    plt.fill_between(x_line, upper, lower, color='black', alpha=0.2)

    plt.xlabel('-1/RT / [mol/J]')
    plt.ylabel('ln(p/p0) / [-]')

def gen_residual_plot():
    residual_y_points = y_points - y_line
    upper_std = avvik
    lower_std = -avvik

    plt.errorbar(x_points,residual_y_points, xerr=s_x_points, yerr=s_y_points,
                 color='black', linestyle='', linewidth=1, capsize=4)

    plt.fill_between(x_points,upper_std,lower_std, color = 'black', alpha=0.15)
    plt.plot(x_points,[0 for x in x_points], color='black')
    plt.xlabel('-1/RT /[mol/J]')
    plt.ylabel('ln(p) - ( a(-1/RT) + b) / [-]')

#### Leser inn data ####
data = data_read('data.txt')
T_points = data[0]
p_points = data[1]

######## Konverterer til SI (p i bar) ######
T_points = T_points + 273.15
p_points = p_points/1000

##### regner ut x- og y-verdier #####
x_points = -1/(k.R *T_points)
s_x_points = regn.usikkerhet_x(T_points)

y_points = np.log(p_points)
s_y_points = regn.usikkerhet_y(p_points)

######### lineærregresjon #########
reg_data = np.polyfit(x_points,y_points, 1, cov=True)
coeff = reg_data[0]
cov = reg_data[1]

x_line = x_points
y_line = coeff[0]*x_line + coeff[1]

s_stigningstall = np.sqrt(cov[0][0])
s_konstantledd = np.sqrt(cov[1][1])

upper = (coeff[0] + s_stigningstall)*x_line + (coeff[1] + s_konstantledd)
lower = (coeff[0] - s_stigningstall)*x_line + (coeff[1] - s_konstantledd)

####### regner ut varians og standardavvik fra linja ####
avvik, dev = regn.varians(y_points=y_points, y_line=y_line)

######## Viser data #######
print('Standardavvik :', avvik)
print()
print('∆vapH :', coeff[0], '±', s_stigningstall)
print('Konstantledd :', coeff[1], '±', s_konstantledd)

jn = input('plot data?')
if jn == 'j':
    gen_plot()
    plt.show()

    jn = input('lagre plot?')
    if jn == 'j':
        gen_plot()
        plt.savefig('ln(p)_mot_T')

jn = input('plot residual?')
if jn == 'j':
    gen_residual_plot()
    plt.show()

    jn = input('lagre residual?')
    if jn == 'j':
        gen_residual_plot()
        plt.savefig('residual_plot')