import numpy as np, matplotlib.pyplot as plt, regning as regn, \
    kalibreringskurve as kalib, konstanter as k, scipy.optimize as sco

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

###### Henter ut data ########
data = read_data('tuva_sin_data.txt')
V_mb_points = data[0]   #[ml]
V_iso_points = data[1]  #[ml]
T_points = data[2] + 273.15      #[K]
p_points = data[3]      #[cmHg]
brytning_k1_points = data[4]
brytning_k2_points = data[5]
brytning_k3_points = data[6]

brytning_v1_points = data[7]
brytning_v2_points = data[8]
brytning_v3_points = data[9]

######## regner ut gjennomsnitt av brytningene #############
brytning_v_points = (brytning_v1_points + brytning_v2_points + brytning_v3_points)/3
brytning_k_points = (brytning_k1_points + brytning_k2_points + brytning_v3_points)/3

########### finner molfraksjoner av iso med standardavvik fra kalibreringskurve ########
x_iso_v, s_x_v = kalib.x_iso_fra_brytning(brytning_v_points)
x_iso_k, s_x_k = kalib.x_iso_fra_brytning(brytning_k_points)

###### molfraksjoner av mb ##########
x_mb_v = 1-x_iso_v
x_mb_k = 1-x_iso_k

######## Regner ut aktivitetskoeffisienter #######
gamma_iso, s_gamma_iso = regn.aktivitetskoeff(y = x_iso_k,s_y=s_x_k, x = x_iso_v, s_x = s_x_v,
                                              T = T_points, p = p_points, komponent='iso')
gamma_mb, s_gamma_mb = regn.aktivitetskoeff(y = x_mb_k,s_y=s_x_k, x = x_mb_v, s_x = s_x_v,
                                            T = T_points, p = p_points, komponent='mb')

aktivitet_mb = gamma_mb*x_mb_v
s_a_mb = np.sqrt((gamma_mb*s_x_v)**2 + (x_mb_v*s_gamma_mb)**2)

aktivitet_iso = gamma_iso*x_iso_v
s_a_iso = np.sqrt((gamma_iso*s_x_v)**2 + (x_iso_v*s_gamma_iso)**2)


jn = input('plot noe?')
if jn == 'j':
    ##### plotter usikkerhet ####
    jn = input('plot usikkerhet?')
    if jn == 'j':
        plt.plot(x_mb_k, s_gamma_mb, color='black', label = 'y', marker='.', markersize='10')
        plt.plot(x_mb_v, s_gamma_mb, color= 'black', linestyle = '--', label = 'x', marker='.', markersize='10')
        plt.legend()
        plt.ylabel(r'$s_{\gamma_{mb}}$ / [-]')
        plt.xlabel('x, y / [-]')
        #plt.savefig('usikkerhet_aktkoeff')
        plt.show()

        plt.bar([x+1 for x in range(len(s_gamma_mb))], s_gamma_mb, color= 'black', width = 0.2)
        plt.xlabel('Måling nr.')
        plt.ylabel(r'$s_{\gamma_{mb}}$ / [-]')
        plt.savefig('usikkerhet_aktkoeff_maalinger')
        plt.show()

    ###### plotter aktivitetskoeffisienter
    jn = input('plot aktkoeff?')
    if jn == 'j':
        plt.errorbar(x_mb_v, gamma_mb, yerr=s_gamma_mb, xerr=s_x_v, color = 'black', linestyle='', capsize = 3)

        plt.xlabel(r'$x_{mb}$')
        plt.ylabel(r'$\gamma_{mb}$')
        plt.savefig('aktivitetskoeffisient_mb_mot_xmb')
        plt.show()

        plt.errorbar(x_iso_v, gamma_iso, yerr=s_gamma_iso, xerr=s_x_v, color='black', linestyle='', capsize=3)

        plt.xlabel(r'$x_{iso}$')
        plt.ylabel(r'$\gamma_{iso}$')
        # plt.savefig('aktivitetskoeffisient_iso_mot_xiso')
        plt.show()

    jn = input('plot aktivitet?')
    if jn == 'j':
        plt.errorbar(x_mb_v, aktivitet_mb, yerr=s_a_mb, xerr=s_x_v, color='black', linestyle='', capsize=3)

        plt.xlabel(r'$x_{mb}$')
        plt.ylabel(r'$a_{mb}$')
        plt.savefig('aktivitet_mb_mot_xmb')
        plt.show()

        plt.errorbar(x_iso_v, aktivitet_iso, yerr=s_a_iso, xerr=s_x_v, color='black', linestyle='', capsize=3)

        plt.xlabel(r'$x_{iso}$')
        plt.ylabel(r'$a_{iso}$')
        plt.savefig('aktivitet_iso_mot_xiso')
        plt.show()

    ####### lager T-x diagram #########
    jn = input('plot T-x?')
    if jn == 'j':
        coeff_k = np.polyfit(x_iso_k,T_points,2)
        coeff_v = np.polyfit(x_iso_v,T_points,2)

        x_akse = np.linspace(0,1,100)
        y_fit_k = np.zeros(len(x_akse))
        y_fit_v = np.zeros(len(x_akse))
        for i in range(len(coeff_k)):
            y_fit_k += coeff_k[i] * (x_akse ** (len(coeff_k) - i-1))
            y_fit_v += coeff_v[i] * (x_akse ** (len(coeff_v) - i-1))

        plt.plot(x_akse, y_fit_k, color='black')
        plt.plot(x_akse, y_fit_v, color='black')

        # plt.plot(x_iso_k, T_points, color='r')
        plt.errorbar(x_iso_k,T_points, xerr=s_x_k, yerr=k.s_T, color = 'black', capsize = 2 , linestyle = '')

        #plt.plot(x_iso_v, T_points, color='b')
        plt.errorbar(x_iso_v,T_points, xerr=s_x_v, yerr=k.s_T, color = 'black', capsize=2,linestyle = '')
        plt.xlabel(r'$x_{iso}$ / [-]')
        plt.ylabel('T / [K]')
        plt.savefig('T_x_med_polyfit')
        plt.show()

    jn = input('x_k som funk av x_v')
    if jn == 'j':
        #plt.scatter(x_iso_v, x_iso_k, color='black')

        x_akse = np.linspace(min(x_iso_v),max(x_iso_v),100)

        coeff = np.polyfit(x_iso_v,x_iso_k,2)
        y_fit = np.polyval(coeff,x_akse)

        #plt.plot(x_akse,y_fit, color = 'black')
        plt.errorbar(x_iso_v, x_iso_k, xerr=s_x_v, yerr=s_x_k, color='black', capsize=2, linestyle='')
        plt.xlabel(r'$x_{iso}$ i væskefase / [-]')
        plt.ylabel(r'$x_{iso}$ i kondensat / [-]')
        plt.savefig('xk_mot_xv_errorbar')
        plt.show()