from math import exp
import matplotlib.pyplot as plt


l1 = 3 * 10**(-3)
l2 = 200 * 10**(-6)
n = 10
m = 5
seg_of_time = 100001
T0 = 300
C0 = 5.54 * 10**4
Cg0 = 20.59 / 18 * 0.86
Tw = 313
Tb = 310
Q = 100
Dg = 2.88 * 10**(-5)
cp = 2000
lambda_l = 0.618
lambda_g = 0.02645
sigma_l = 1.5 * 10**(-7)
condens = 0
dx1 = (l1/(l1+l2)) / n
dx2 = (l2/(l1+l2)) / m
dt = 100000 / seg_of_time
Tl = [[Tw/Tw if i==0 else T0/Tw for i in range(n)] for _ in range(seg_of_time)]
Tg = [[Tb/Tw if i==m-1 else T0/Tw for i in range(m)] for _ in range(seg_of_time)]
Cg = [[Cg0/C0] * m for _ in range(seg_of_time)]

def vaper_pressure(T):
	a = -7.7645
	b = 1.45838
	c = -2.7758
	d = -1.2303
	x = 1 - T*Tw/647.3
	return 22120000 * exp((a*x+b*x**1.5+c*x**3+d*x**6)/(1-x))

def vaper_c(T):
	R = 8.314
	return vaper_pressure(T)/(R*T*Tw) / C0

def calc_rho(C):
	return 18/1000*C*C0

Cg[0][0] = vaper_c(T0/Tw)

for i in range(1, seg_of_time):
	for j in range(1, n-1):
		Tl[i][j] = Tl[i-1][j] + sigma_l/10 * ((Tl[i-1][j-1]+Tl[i-1][j+1]-2*Tl[i-1][j])/dx1**2) * dt
	Tl[i][n-1] = Tl[i][n-2] + lambda_g/lambda_l * dx1/dx2 * (Tg[i-1][1]-Tg[i-1][0]) + Tw/lambda_l * Dg/C0 * dx1/dx2 * (Cg[i-1][1]-Cg[i-1][0]) * Q
	Tg[i][0] = Tl[i][n-1]
	for j in range(1, m-1):
		rho = calc_rho(Cg[i-1][j])
		sigma_g = lambda_g/(rho*cp)
		Tg[i][j] = Tg[i-1][j] + sigma_g/10 * ((Tg[i-1][j-1]+Tg[i-1][j+1]-2*Tg[i-1][j])/dx2**2) * dt

	Cg[i][0] = vaper_c(Tg[i][0])
	for j in range(1, m-1):
		C_calc = Cg[i-1][j] + Dg/10 * ((Cg[i-1][j-1]+Cg[i-1][j+1]-2*Cg[i-1][j])/dx2**2) * dt
		if C_calc < vaper_c(Tg[i][j]):
			Cg[i][j] = C_calc
		else:
			Cg[i][j] = vaper_c(Tg[i][j])
			condens += (C_calc-vaper_c(Tg[i][j])) * dx2
	#condens += Cg[i-1][m-1] + Dg/10 * ((Cg[i-1][m-2]-Cg[i-1][m-1])/dx2**2) * dt - Cg[i][m-2]
	Cg[i][m-1] = Cg[i][m-2]
for i in range(6):
	plt.plot(range(n), Tl[10**i])
print(condens)
plt.show()