import numpy as np
import matplotlib.pyplot as plt

# Data for 2-metyl-1-butanol
Mm_MB = 88.2
rho_MB = 0.815
dH_MB = 54000.0
dH_MB_err = 0.05
T_b_MB = 401.15
T_b_MB_err = 0.005

# Data for ispopropanol
Mm_IP = 60.1
rho_IP = 0.781
dH_IP = 45000
T_b_IP = 355.45

# Kjente verdier
R = 8.3145
R_err = 0.00005
p = 1.0037
p_err = 0.00066

# Målte verdier
V_MB = np.array([0, 1, 2, 4, 5, 5, 5, 5])
V_MB_err = np.array([0.0] + [0.15] * (len(V_MB) - 1))
V_IP = np.array([5, 5, 5, 5, 4, 2, 1, 0])
V_IP_err = np.array([0.15] * (len(V_IP) - 1) + [0.0])

n_MB = rho_MB * V_MB / Mm_MB
n_MB_err = rho_MB * V_MB_err / Mm_MB
n_IP = rho_IP * V_IP / Mm_IP
n_IP_err = rho_IP * V_IP_err / Mm_IP

x_MB = n_MB / (n_MB + n_IP)
x_MB_err = np.sqrt(
    (n_IP * n_MB_err / (n_MB + n_IP) ** 2) ** 2 +
    (- n_IP * n_IP_err / (n_MB + n_IP) ** 2) ** 2
)

eta = np.array([1.3762, 1.3762, 1.3809, 1.3851, 1.3900, 1.3967, 1.4008, 1.4094])    # Åsnekok
eta_err = np.ones(len(eta)) * 0.00005

# Kalibreringskurve
pol, matrix = np.polyfit(x_MB, eta, 1, cov=True, full=False)
a = pol[0]
b = pol[1]
a_err = np.sqrt(matrix[0][0])
b_err = np.sqrt(matrix[1][1])


def get_mol_fraction(refraction_index):
    return (refraction_index - b) / a


def plot_calibration_curve():
    plt.scatter(x_MB, eta, color='black')
    plt.plot(get_mol_fraction(eta), eta, color='black')
    plt.errorbar(x_MB, eta, xerr=x_MB_err*2, yerr=eta_err*2, fmt='o', color='black')

    plt.xlabel(r'Molfraksjon 2-metyl-1-butanol i væskefase $x_{MB}$ [-]')
    plt.ylabel("Brytningsindeks \u03B7 [-]")
    ax = plt.gca()
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
   # plt.show()
    #plt.savefig("kalibreringskurve")


# Målte verdier
T_k = np.array([118.0, 99.8, 95.0, 88.2, 78.0, 80.2, 83.6, 93]) + 273.15
T_err = np.ones(len(T_k)) * 0.1
eta_v = np.array([1.40835, 1.4056, 1.40515, 1.4027, 1.3752, 1.38185, 1.39075, 1.40095])
eta_v_err = np.ones(len(eta_v)) * 0.00005
eta_k = np.array([1.4083, 1.38945, 1.38495, 1.38035, 1.37495, 1.37595, 1.3771, 1.3842])
eta_k_err = np.ones(len(eta_v)) * 0.00005


# Finne molfraksjon MB i væske (x) og gass (y)
x = get_mol_fraction(eta_v)
x_err = np.sqrt(
    (eta_v_err / a)**2 +
    ((b - eta_v) * a_err / a**2)**2 +
    (b_err / b)**2
)

y = get_mol_fraction(eta_k)
y_err = np.sqrt(
    (eta_k_err / a)**2 +
    ((b - eta_k) * a_err / a**2)**2 +
    (b_err / b)**2
)


# Aktivitetskoeffisienter
ln_gamma = (dH_MB/R) * ((1/T_k) - (1/T_b_MB)) + np.log(p * y / x)
ln_gamma_err = np.sqrt(
    (((1 / (R*T_k)) - (1/(R*T_b_MB))) * dH_MB_err) ** 2 +
    (((1 / (R**2 * T_k)) - (1 / (R**2 * T_b_MB))) * R_err) ** 2 +
    (- dH_MB * T_err / (R * T_k**2)) ** 2 +
    (-dH_MB * T_b_MB_err / (R * T_b_MB**2)) ** 2 +
    (p_err / p) ** 2 +
    (- a_err / a) ** 2 +
    (- b_err / b) ** 2

)
gamma = np.exp(ln_gamma)
gamma_err = np.exp(ln_gamma) * ln_gamma_err

# Aktiviteter
activity = gamma * x
activity_err = np.sqrt(
    (x * gamma_err) ** 2 +
    (gamma * x_err) ** 2
)


def plot_activity_coefficient():
    plt.scatter(x, gamma, color='black')
    plt.errorbar(x, gamma, xerr=x_err*2, yerr=gamma_err*2, fmt='o', color='black')

    plt.xlabel(r'Molfraksjon 2-metyl-1-butanol i væskefase $x_{MB}$ [-]')
    plt.ylabel("Aktivitetskoeffisient \u03B3 [-]")
    ax = plt.gca()
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    plt.show()
    #plt.savefig("aktivitetskoeffisient")


def plot_activity():
    #plt.scatter(x, activity, color='black')
    plt.errorbar(x, activity, xerr=x_err*2, yerr=activity_err*2, fmt='o', color='black')

    plt.xlabel(r'Molfraksjon 2-metyl-1-butanol i væskefase $x_{MB}$ [-]')
    plt.ylabel(r'Aktivitet $a_{MB}$ [-]')
    ax = plt.gca()
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    plt.show()
    # plt.savefig("aktivitet")


def plot_makabe_tiele():
    plt.scatter(x, y, color='black')
    plt.errorbar(x, y, xerr=x_err*2, yerr=y_err*2, fmt='o', color='black')
    plt.plot([0, 1], [0, 1], color='black')

    plt.xlabel(r'Molfraksjon 2-metyl-1-butanol i væskefase $x_{MB}$ [-]')
    plt.ylabel(r'Molfraksjon 2-metyl-1-butanol i gassfase $y_{MB}$ [-]')
    ax = plt.gca()
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    #plt.show()
    plt.savefig("makabe_tiele")


def plot_phase_diagram():
    plt.errorbar(x, T_k, xerr=x_err*2, yerr=T_err*2, fmt='o', color='black')
    plt.errorbar(y, T_k, xerr=y_err*2, yerr=T_err*2, fmt='o', color='black')

    p1 = np.polyfit(x, T_k, 2)
    p2 = np.polyfit(y, T_k, 2)

    x_blaaah = np.linspace(0, 1, 50)
    plt.plot(x_blaaah, np.polyval(p1, x_blaaah), color='black')
    plt.plot(x_blaaah, np.polyval(p2, x_blaaah), color='black')

    plt.xlabel(r'Molfraksjon 2-metyl-1-butanol $x_{MB}$ [-]')
    plt.ylabel(r'Koketemperatur for blandingen $T_b$ [K]')
    ax = plt.gca()
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    #plt.show()
    plt.savefig("fasediagram")


def print_errors():
    print("a:", a_err)
    print("b:", b_err)
    print("eta:", eta_err)
    print("eta_k", eta_k_err)
    print("eta_v", eta_v_err)
    print("p:", p_err)
    print("x:", x_err)
    print("y:", y_err)


def print_activity_coefficients():
    print("\n\n")
    print("Molfrak MB".center(50) + "|".center(3) + "Aktivitetskoeffisient".center(25) +
          "|".center(3) + "Usikkerhet (2 std.avvik)".center(30))
    print("-" * 78)
    for i in (range(len(gamma))):
        print((str(x[i]) + "  +-  " + str(x_err[i] * 2)).center(50) + "|".center(3) + str(gamma[i]).center(25) +
              "|".center(3) + str(gamma_err[i] * 2).center(30))
    print("\n\n")


def print_pol_coeffs():
    print("a:", a, "+-", a_err)
    print("b:", b, "+-", b_err)


def main():
    print_activity_coefficients()
    print_pol_coeffs()
    print_errors()
    plot_activity_coefficient()
    plt.clf()
    plot_activity()
    plt.clf()
    plot_calibration_curve()
    plt.clf()
    plot_makabe_tiele()
    plt.clf()
    plot_phase_diagram()


main()