# -*- coding: utf-8 -*-

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


# color range
colors = [
    "#00395b", "#74adc1", "#b54036", "#ec6707", "#bfbfbf", "#999999", "#010101"
]


# %% figure 1: plot component exergy destruction

# color range for E_F and E_P bars
E_F_colors = ["#3a9dce", "#f08e2b", "#f08e2b", "#f08e2b", "#f08e2b", "#db5252"]

# read data
df_ED_NH3 = pd.read_csv("NH3_E_D.csv", index_col=0)
df_ED_R290 = pd.read_csv("R290_E_D.csv", index_col=0)

# get data from dataframes
y_NH3 = df_ED_NH3.columns.values.tolist()
y_NH3_pos = np.arange(len(y_NH3))
E_P_NH3 = df_ED_NH3.loc["E_P"]
E_D_NH3 = df_ED_NH3.loc["E_D"]
y_R290 = df_ED_R290.columns.values.tolist()
y_R290_pos = np.arange(len(y_R290))
E_P_R290 = df_ED_R290.loc["E_P"]
E_D_R290 = df_ED_R290.loc["E_D"]

# create bar diagram of absolute exergy destruction
fig, axs = plt.subplots(2, 2, constrained_layout=True,
                        sharex="row", sharey="row")
axs[0, 0].barh(y_NH3_pos, E_P_NH3, align="center", color=E_F_colors)
axs[0, 0].barh(y_NH3_pos, E_D_NH3, align="center", left=E_P_NH3, label="E_D",
               color="#6ed880")
axs[0, 0].set_xlabel("Exergy in W")
axs[0, 0].set_yticks(y_NH3_pos)
axs[0, 0].set_yticklabels(y_NH3)
axs[0, 0].invert_yaxis()
axs[0, 0].set_title("NH3")

axs[0, 1].barh(y_R290_pos, E_P_R290, align="center", color=E_F_colors)
axs[0, 1].barh(y_R290_pos, E_D_R290, align="center", left=E_P_R290,
               color="#6ed880")
axs[0, 1].set_xlabel("Exergy in W")
axs[0, 1].set_yticks(y_R290_pos)
axs[0, 1].set_yticklabels(y_R290)
axs[0, 1].set_title("R290")
axs[0, 1].set_xlim(right=1000)

# create bar diagram of percentage exergy destruction
E_F_NH3 = df_ED_NH3.loc["E_P", "E_F"]
E_F_R290 = df_ED_R290.loc["E_P", "E_F"]


axs[1, 0].barh(y_NH3_pos, E_P_NH3/E_F_NH3, align="center", color=E_F_colors)
axs[1, 0].barh(y_NH3_pos, E_D_NH3/E_F_NH3, align="center",
               left=E_P_NH3 / E_F_NH3, color="#6ed880")
axs[1, 0].set_xlabel("$\epsilon$")
axs[1, 0].set_yticks(y_NH3_pos)
axs[1, 0].set_yticklabels(y_NH3)
axs[1, 0].invert_yaxis()

axs[1, 1].barh(y_R290_pos, E_P_R290/E_F_R290, align="center",
               color=E_F_colors)
axs[1, 1].barh(y_R290_pos, E_D_R290/E_F_R290, align="center",
               left=E_P_R290 / E_F_R290, color="#6ed880")
axs[1, 1].set_xlabel("$\epsilon$")
axs[1, 1].set_yticks(y_R290_pos)
axs[1, 1].set_yticklabels(y_R290)
axs[1, 1].set_xlim(right=1)

fig.suptitle("Component Exergy Destruction", fontsize=14)
fig.legend(bbox_to_anchor=(0.08, 0, 0.9, .0), loc="center",
           ncol=3, borderaxespad=0.)
fig.savefig("diagram_E_D.svg", bbox_inches="tight")
plt.close()


# %% figure 2: epsilon depending on ambient Temperature Tamb
#              and mean geothermal temperature Tgeo

Tamb_design = 1
Tgeo_design = 10

Tamb_range = [4, 8, 12, 16, 20]
Tgeo_range = [14, 13, 12, 11, 10, 9, 8]

# read data
df_eps_Tamb_NH3 = pd.read_csv("NH3_eps_Tamb.csv", index_col=0)
df_eps_Tgeo_NH3 = pd.read_csv("NH3_eps_Tgeo.csv", index_col=0)
df_eps_Tamb_R290 = pd.read_csv("R290_eps_Tamb.csv", index_col=0)
df_eps_Tgeo_R290 = pd.read_csv("R290_eps_Tgeo.csv", index_col=0)

# create plot
fig, axs = plt.subplots(1, 2, constrained_layout=True,
                        sharex="col", sharey="all")
axs[0].plot(Tamb_range, df_eps_Tamb_NH3.loc[Tgeo_design], "x",
            label="NH3", color=colors[0], markersize=7)
axs[0].plot(Tamb_range, df_eps_Tamb_R290.loc[Tgeo_design], "x",
            label="R290", color=colors[3], markersize=7)
axs[0].set_title("ambient Temperature")
axs[0].set_ylabel("$\epsilon$")
axs[0].set_ylim([0.25, 0.6])
axs[0].set_xlabel("$T_{amb}$ in °C ($T_{geo}$ = " + str(Tgeo_design) + "°C)")
axs[0].set_ylabel("$\epsilon$")
axs[0].legend(loc="lower left")

axs[1].plot(Tgeo_range, df_eps_Tgeo_NH3.loc[Tamb_design], "x",
            label="NH3", color=colors[0], markersize=7)
axs[1].plot(Tgeo_range, df_eps_Tgeo_R290.loc[Tamb_design], "x",
            label="R290", color=colors[3], markersize=7)
axs[1].set_title("geothermal Temperature")
axs[1].set_xlabel("$T_{geo}$ in °C ($T_{amb}$ = 2.8°C)")
axs[1].legend(loc="lower left")

fig.suptitle("Exergetic Efficency depending on Temperature", fontsize=14)
fig.savefig("diagram_eps_Tamb_Tgeo.svg", bbox_inches="tight")
plt.close()


# %% figure 3: epsilon and COP depending on mean geothermal temperature Tgeo
#              and mean heating system temperature Ths

Tgeo_range = [14, 12, 10]
Ths_range = [45, 40, 35]

# read data
df_cop_Tgeo_Ths_NH3 = pd.read_csv("NH3_cop_Tgeo_Ths.csv", index_col=0)
df_eps_Tgeo_Ths_NH3 = pd.read_csv("NH3_eps_Tgeo_Ths.csv", index_col=0)
df_cop_Tgeo_Ths_R290 = pd.read_csv("R290_cop_Tgeo_Ths.csv", index_col=0)
df_eps_Tgeo_Ths_R290 = pd.read_csv("R290_eps_Tgeo_Ths.csv", index_col=0)

# create plot
fig, axs = plt.subplots(2, 2, constrained_layout=True,
                        sharex="all", sharey="row")
i = 0
for Tgeo in Tgeo_range:
    axs[0, 0].plot(Ths_range, df_eps_Tgeo_Ths_NH3.loc[Tgeo], "x",
                   color=colors[i], label="$T_{geo}$ = " + str(Tgeo) +
                   " °C", markersize=7, linewidth=2)
    axs[1, 0].plot(Ths_range, df_cop_Tgeo_Ths_NH3.loc[Tgeo], "x",
                   color=colors[i], markersize=7, linewidth=2)
    axs[0, 1].plot(Ths_range, df_eps_Tgeo_Ths_R290.loc[Tgeo], "x",
                   color=colors[i], markersize=7, linewidth=2)
    axs[1, 1].plot(Ths_range, df_cop_Tgeo_Ths_R290.loc[Tgeo], "x",
                   color=colors[i], markersize=7, linewidth=2)
    i += 1

axs[0, 0].set_ylabel("$\epsilon$")
axs[0, 0].set_ylim([0.4, 0.6])
axs[0, 0].set_title("NH3")

axs[1, 0].set_ylabel("COP")
axs[1, 0].set_xlabel("$T_{heating system}$ in °C")
axs[1, 0].set_ylim([3.5, 6.5])

axs[0, 1].set_title("R290")

axs[1, 1].set_xlabel("$T_{heating system}$ in °C")

fig.legend(bbox_to_anchor=(0.08, -0.1, 0.9, .0), loc="lower left",
           ncol=3, mode="expand", borderaxespad=0.)
fig.suptitle(
    "COP and $\epsilon$ depending heating and geothermal mean temperature",
    fontsize=12)
fig.savefig("diagram_cop_eps_Tgeo_Ths.svg", bbox_inches="tight")
plt.close()


# %% figure 4: epsilon and COP depending on mean geothermal temperature Tgeo
#              and heating load Q_cond

Q_range = [4.3e3, 4e3, 3.7e3, 3.4e3, 3.1e3, 2.8e3]
Tgeo_range = [14, 12, 10]

# read data
df_cop_Tgeo_Q_NH3 = pd.read_csv("NH3_cop_Tgeo_Q.csv", index_col=0)
df_eps_Tgeo_Q_NH3 = pd.read_csv("NH3_eps_Tgeo_Q.csv", index_col=0)
df_cop_Tgeo_Q_R290 = pd.read_csv("R290_cop_Tgeo_Q.csv", index_col=0)
df_eps_Tgeo_Q_R290 = pd.read_csv("R290_eps_Tgeo_Q.csv", index_col=0)

# create plot
fig, axs = plt.subplots(2, 2, constrained_layout=True,
                        sharex="all", sharey="row")

i = 0
for Tgeo in Tgeo_range:
    axs[0, 0].plot(Q_range, df_eps_Tgeo_Q_NH3.loc[Tgeo], "x",
                   color=colors[i], label="$T_{geo}$ = " + str(Tgeo) +
                   " °C", markersize=7, linewidth=2)
    axs[1, 0].plot(Q_range, df_cop_Tgeo_Q_NH3.loc[Tgeo], "x",
                   color=colors[i], markersize=7, linewidth=2)
    axs[0, 1].plot(Q_range, df_eps_Tgeo_Q_R290.loc[Tgeo], "x",
                   color=colors[i], markersize=7, linewidth=2)
    axs[1, 1].plot(Q_range, df_cop_Tgeo_Q_R290.loc[Tgeo], "x",
                   color=colors[i], markersize=7, linewidth=2)
    i += 1

axs[0, 0].set_ylabel("$\epsilon$")
axs[0, 0].set_ylim([0.4, 0.6])
axs[0, 0].set_title("NH3")

axs[1, 0].set_ylabel("COP")
axs[1, 0].set_xlabel(r"$Q_{cond}$ in kW")
axs[1, 0].set_ylim([3, 6])
axs[1, 0].set_xlim([2500, 4500])

axs[0, 1].set_title("R290")

axs[1, 1].set_xlabel(r"$Q_{cond}$ in kW")

plt.xticks(np.arange(2500, 5000, step=1000), np.arange(2.5, 5, step=1))
fig.legend(bbox_to_anchor=(0.08, -0.1, 0.9, .0), loc="lower left",
           ncol=3, mode="expand", borderaxespad=0.)
fig.suptitle(
    "COP and $\epsilon$ depending on load and geothermal mean temperature",
    fontsize=12)
fig.savefig("diagram_cop_eps_Tgeo_Q.svg", bbox_inches="tight")
plt.close()