import matplotlib.pyplot as plt
import scipy.optimize as opt

# VdW réduite (π+3/β²)(3β-1)=8τ

def pressionRedVdW (tRed,vRed) :
    return 8*tRed/(3*vRed-1)-3/vRed**2

N=100
VredMax=4
h=(VredMax-1/3)/N
volumes=[1/3+(i+1)*h for i in range(N)]

temperatures=[0.8,0.9,1.0,1.1,1.3,1.5,1.8]

plt.ylim(ymax=2)
for tRed in temperatures :
    pressions = [pressionRedVdW(tRed,vRed) for vRed in volumes]
    plt.plot(volumes,pressions,label=str(tRed))
plt.legend()
plt.show()

# donne l'éq vérifiée par le min et le max de la courbe
# il y a 3 solutions : une avant 1/3, impossible, les deux autres correctes
# le fsolve avec 1 donne toujours la plus petite, pour tout t<1
# avec 2/tRed, on a toujours la plus grande (même pour tRed=0.001)
def VdWminMax(v,tRed):
    return 4*tRed*v**3-(3*v-1)**2

def minMaxVdW (tRed):
    if tRed>=1 :
        return -1
    vMin=opt.fsolve(VdWminMax,1,args=(tRed))[0]
    vMax=opt.fsolve(VdWminMax,2/tRed,args=(tRed))[0]
    return (vMin,vMax)