Thermodynamic Evaluation of Green Energy Technologies

import numpy as np
import matplotlib.pyplot as plt

from ipywidgets import interact
from ipywidgets.widgets import FloatSlider
from scipy.spatial import ConvexHull

# Define variables and constants
R = 8.314
x = np.linspace(0, 1, endpoint=False, num=301)[1:]
# T = 365
# # mu_a1 = 1400
# mu_a1 = 2000
# # mu_b1 = 0
# mu_b1 = 1500
# mu_a2 = 0
# mu_b2 = 800
# l1 = 4400
# l2 = 4400

def plot_lower_convexhull(T, mu_a1, mu_b1, mu_a2, mu_b2, l1, l2):
    mu_1 = mu_a1 * (1 - x) + mu_b1 * x + R * T * (x * np.log(x) + (1-x) * np.log(1-x)) + l1 * x * (1-x)
    mu_2 = mu_a2 * (1 - x) + mu_b2 * x + R * T * (x * np.log(x) + (1-x) * np.log(1-x)) + l2 * x * (1-x)

    combined_xs = np.concatenate((x, x))
    combined_mus = np.concatenate((mu_1, mu_2))
    plot_buffer = 200
    max_y_plot = np.max(combined_mus) + plot_buffer
    min_y_plot = np.min(combined_mus) - plot_buffer

    upper_left_x, upper_left_y = np.min(x), max_y_plot
    upper_right_x, upper_right_y = np.max(x), max_y_plot

    convexhull_x = np.concatenate((combined_xs, [upper_left_x, upper_right_x]))
    convexhull_y = np.concatenate((combined_mus, [upper_left_y, upper_right_y]))
    convexhull_xy = np.vstack(([convexhull_x], [convexhull_y])).T
    # number of convex hull points
    num_chpts = convexhull_xy.shape[0]

    ch = ConvexHull(convexhull_xy)
    lower_ch_vertices = np.array([_ for _ in ch.vertices if _ not in [num_chpts-1, num_chpts-2]], dtype=np.int32)
    # lower_ch_vertices = np.sort(lower_ch_vertices)

    # Find all the compositions where there are a common tangent line
    index_diffs = np.ediff1d(lower_ch_vertices)
    # the common tangents occurs at index jump greater than 1
    common_tangents_diff_ids = np.where(index_diffs != 1)[0]
    common_tangents_ids = [[_, _ + 1] for _ in common_tangents_diff_ids]

    fig, ax = plt.subplots(figsize=(6, 4))
    ax.plot(x, mu_1, label="$g_1$")
    ax.plot(x, mu_2, label="$g_2$")

    # Filter out indices jump that correspond to jumping from one end of one phase to the other end of the same/other phase
#     num_lowch_pts = lower_ch_vertices.shape[0]
#     lower_ch_vertices_for_plotting = lower_ch_vertices.copy()
#     if lower_ch_vertices[:2] in [[num_lowch_pts-1, 0],
#                                  [2*num_lowch_pts-1, num_lowch_pts],
#                                  [0, num_lowch_pts],
#                                  [0, 2*num_lowch_pts-1],
#                                  [num_lowch_pts-1, num_lowch_pts],
#                                  [num_lowch_pts-1, 2*num_lowch_pts-1]]:
#          lower_ch_vertices_for_plotting = np.concatenate((lower_ch_vertices[1:], [lower_ch_vertices[0]]))
    ax.plot(convexhull_x[lower_ch_vertices], convexhull_y[lower_ch_vertices],
            c="purple", label="lower convexhull", ls="-.")

    # print(lower_ch_vertices)
    for pair_of_common_tangent in common_tangents_ids:
        # if pair_of_common_tangent not in [[num_lowch_pts-1, 0],
        #                                   [2*num_lowch_pts-1, num_lowch_pts],
        #                                   [0, num_lowch_pts],
        #                                   [0, 2*num_lowch_pts-1],
        #                                   [num_lowch_pts-1, num_lowch_pts],
        #                                           [num_lowch_pts-1, 2*num_lowch_pts-1]]:
            ax.plot(convexhull_x[lower_ch_vertices[pair_of_common_tangent]],
                    convexhull_y[lower_ch_vertices[pair_of_common_tangent]],
                    c="r", label="common tangent")
            ax.scatter(convexhull_x[lower_ch_vertices[pair_of_common_tangent]],
                       convexhull_y[lower_ch_vertices[pair_of_common_tangent]],
                       c="r")
    ax.set_xlim(0, 1)
    ax.set_ylim(-20000, 20000)
    plt.legend()
    plt.show()

interact(plot_lower_convexhull,
         T=FloatSlider(value=365, min=0, max=5000),
         mu_a1=FloatSlider(value=2000, min=-10000, max=10000),
         mu_b1=FloatSlider(value=1500, min=-10000, max=10000),
         mu_a2=FloatSlider(value=0, min=-10000, max=10000),
         mu_b2=FloatSlider(value=800, min=-10000, max=10000),
         l1=FloatSlider(value=4400, min=-30000, max=30000),
         l2=FloatSlider(value=4400, min=-30000, max=30000))