import numpy as np
import matplotlib.pyplot as plt

# Entropy equation
def s(N, V, U):
return (N * V *U) ** (1/3)

# Verify a) as sanity check
s(50, 13, 80)

# Vary U_a from 0 to 80 in intervals of 0.1
U_a = np.linspace(0, 80, 801)

# PART B
# Find S_a and S_b as functions of U_a and U - U_a respectively
S_a = s(30, 9, U_a)
S_b = s(20, 4, 80 - U_a)
# Add them up to find total entropy of system
S_tot = S_a + S_b
S_tot

# Plot S_tot vs U_a
plt.plot(U_a, S_tot)
plt.xlabel("U_a")
plt.ylabel("Total Entropy")
plt.show()

# Find U_a for which total entropy was maximized
U_a[np.argmax(S_tot)]

# Max value of S_tot
S_tot[np.argmax(S_tot)]

# PART E
# Change volumes
# Find S_a and S_b as functions of U_a and U - U_a respectively
S_a = s(30, 7.8, U_a)
S_b = s(20, 5.2, 80 - U_a)
# Add them up to find total entropy of system
S_tot = S_a + S_b
S_tot

# Max value of S_tot
S_tot[np.argmax(S_tot)]

# Find U_a for which total entropy was maximized
U_a[np.argmax(S_tot)]

def T(N, V, U):
return 3 * U ** (2/3) / ( (N * V) ** (1/3))

# Temp in part a)
T(50, 13, 80)

# Temp in part e)
T_a = T(30, 7.8, 48.0)
T_b = T(20, 5.2, 32.0)
T_a, T_b