from scipy.optimize import fsolve
from scipy import constants as k
from math import log, exp

# SI Units
# R = 8.31447
# Temperature = Kelvin
# Pressure = kPaA

def EOS(T,P):

	### Propane Properties ###

	MW = 44.09700012 # kg/kg.mol
	Tc = 369.8980103 # Kelvin
	Pc = 4256.660156 # kPaA
	w = 0.152400002

	#Imperial Units BWRS Coefficients:
	Bo, Ao, Co, Do, Eo, b, a, d, alpha, c, gamma = 0.964762008,18634.69982,7961779083,4.53707E+11,2.56053E+13,5.462687934,40066.39826,15051999.65,2.014020806,27446095153,4.56182043

	#SI Units BWRS Coefficients:
	Bo, Ao, Co, Do, Eo, b, a, d, alpha, c, gamma = 0.060228125,500.7254656,66030166.22,2090428137,65541613091,0.021289462,67.21043849,14027.43031,0.000490006,14209931.35,0.017778556

	### BWRS-EOS ###

	def func(rho):
		return P - ( rho*k.R*T + (Bo*k.R*T - Ao - Co/T**2 + Do/T**3 - Eo/T**4)*rho**2 + (b*k.R*T - a - d/T)*rho**3 + alpha*(a + d/T)*rho**6 + c*rho**3/T**2*(1 + gamma*rho**2)*exp(-gamma*rho**2) )

	rho_molar = fsolve(func, 50)[0]  # m3/kg.mol
	molar_volme = 1/rho_molar*1000   # cm3/mol
	rho_mass = rho_molar*MW          # kg/m3
	Z = P/rho_molar/k.R/T

	# Departure Functions per K.E Starling
	# REPORT ORO-5249-2 DEVELOPMENT OF WORKING FLUID THERMODYNAMIC PROPERTIES INFORMATION FOR GEOTHERMAL CYCLES - PHASE I
	# Enthalpy Departure - Equation 6
	# Entropy Departure - Equation 10

	enthalpy_departure = (Bo*k.R*T - 2*Ao - 4*Co/T**2 + 5*Do/T**3 - 6*Eo/T**4)*rho_molar + 0.5*(2*b*k.R*T - 3*a - 4*d/T)*rho_molar**2 + 1/5*alpha*(6*a + 7*d/T)*rho_molar + c/gamma/T**2 * (3 - (3 + 0.5*gamma*rho_molar**2 - gamma**2 * rho_molar**4)*exp(-gamma*rho_molar**2))
	entropy_departure = -k.R*log(rho_molar*k.R*T) - (Bo*k.R + 2*Co/T**3 - 3*Do/T**4 + 4*Eo/T**5)*rho_molar - 0.5*(b*k.R + d/T**2)*rho_molar**2 + alpha*d*rho_molar**5/(5*T**2) + 2*c/(gamma*T**3)*(1-(1+0.5*gamma*rho_molar**2)*exp(-gamma*rho_molar**2))

	# Cp Ideal from Aspen HYSYS

	CpID = 8.43833E-31*T**6 + -8.61031E-27*T**5 + 3.70274E-11*T**4 + -1.17682E-07*T**3 + 5.24515E-05*T**2 + 0.186450054*T + 17.41831505
	CvID = CpID - 8.31447

	print(entropy_departure)

	return molar_volme, enthalpy_departure, entropy_departure, CpID, CvID

state_1 = EOS(378.15, 500)
state_2 = EOS(463.15, 2500)

enthalpy_delta = 0.5*(state_1[3]+state_2[3])*(463.15-378.15) + (state_2[1]-state_1[1])
entropy_delta = 0.5*(state_1[3]+state_2[3])*log(463.15/378.15) + (state_2[2]-state_1[2])


print(f"State 1 V (cm3/mol): {state_1[0]}")
print(f"State 2 V (cm3/mol): {state_2[0]}")
print(f"Enthalpy Delta kJ/kg.mol: {enthalpy_delta}")
print(f"Entropy Delta kJ/kg.mol: {entropy_delta}")