# -*- coding: utf-8 -*-


## Writer's note: I know this script is absolute, AIDS-Ridden fucking cancer 
## insofar as scriptspergs are concerned, I don't care.

import numpy as np

## Define Variables
## Cube A

##Total number of Moles
A_n=12

## moles N2 in A
A_n_N2=10

## moles of O2 in A
A_n_O2=2

## moles of Cl2 in a
A_n_Cl2=0

## mole Fraction of O2
A_y_O2=0.1667

## mole fraction of N2
A_y_N2=0.8333

## mole fraction of Cl2
A_y_Cl2=0

## Temperature, in Kelvin
A_T=293.15

## Volume, in Liters
A_V=1000

##Molar Volume
A_Vm=A_V/A_n

## Van Der Waal's constants for Nitrogen
A_N2_VDWa=1.352

A_N2_VDWb=0.0387

## Van Der Waal's constants for Oxygen
A_O2_VDWa=1.364

A_O2_VDWb=0.0319

## SOURCE : www2.ucdsb.on.ca/tiss/stretton/database/van_der_waals_constants.html

## Calculate The mixture VDW constants
A_VDWa=(A_N2_VDWa*A_y_N2)+(A_O2_VDWa*A_y_O2)

A_VDWb=(A_N2_VDWb*A_y_N2)+(A_O2_VDWb*A_y_O2)

## Heat Capacities:
A_N2_Cp=29.12379083

A_O2_Cp=29.38262912

## Cube B

## Total moles 
B_n=13

## moles N2 in A
B_n_N2=10

## moles of O2 in A
B_n_O2=2

## moles of Cl2 in B
B_n_Cl2=1

## Mole Fraction of O2
B_y_O2=0.1538

## mole fraction of N2
B_y_N2=0.7692

## mole fraction Cl2
B_y_Cl2=0.0770

## Temperature, in Kelvin
B_T=313.15

## Volume, Liters
B_V=1000

## molar Volume
B_Vm=B_V/B_n

## Van Der Waal's constants for Nitrogen
B_N2_VDWa=1.352

B_N2_VDWb=0.0387

## Van Der Waal's constants for Oxygen
B_O2_VDWa=1.364

B_O2_VDWb=0.0319

## Van Der Waal's constants for Chlorine
B_Cl2_VDWa=6.26

B_Cl2_VDWb=0.0542

## SOURCE : www2.ucdsb.on.ca/tiss/stretton/database/van_der_waals_constants.html

## Calculate The mixture VDW constants
B_VDWa=(B_N2_VDWa*B_y_N2)+(B_O2_VDWa*B_y_O2)+(B_Cl2_VDWa*B_y_Cl2)

B_VDWb=(B_N2_VDWb*B_y_N2)+(B_O2_VDWb*B_y_O2)+(B_Cl2_VDWb*B_y_Cl2)

## SOURCE: http://newmanator.net/projects/Mixing_Rules.pdf

## Heat Capacities

B_N2_Cp=29.13296396

B_O2_Cp=29.4639599

B_Cl2_Cp=34.20531862

## Property Calculations

## Section 0.1: Pressure Definitions

## Van Der Waal's Equation of State to derive pressure function
def VDWP(a,b,T,V,n) :
    return (((0.08314*T)/(V-b)))-(a/(V**2))

## Summation of Partial Pressures
## Sanity test on function, expect ~ 0.25 bar
#print(VDWP(A_N2_VDWa,A_N2_VDWb,A_T,A_Vm_N2,A_n_N2))
## PASS

## P_A Calculation:
P_A=VDWP(A_VDWa,A_VDWb,A_T,A_Vm, A_n)
print('Pressure in Cube A is:')
print(P_A)
print('bar')

print(' ')

## P_B Calculation:
P_B=VDWP(B_VDWa,B_VDWb,B_T,B_Vm, B_n)
print('Pressure in Cube B is:')
print(P_B)
print('bar')

## Equilibrium VDW constant calculations
## new mole fracitons:
y_N2=20/25
y_O2=4/25
y_Cl2=1/25

Eq_VDWa=(B_N2_VDWa*y_N2)+(B_O2_VDWa*y_O2)+(B_Cl2_VDWa*y_Cl2)

Eq_VDWb=(B_N2_VDWb*y_N2)+(B_O2_VDWb*y_O2)+(B_Cl2_VDWb*y_Cl2)

Eq_P=VDWP(Eq_VDWa,Eq_VDWb,308.9,(2000/25),25)

Eq_P_heat=VDWP(Eq_VDWa,Eq_VDWb,308.9,(2000/25),25)-VDWP(Eq_VDWa,Eq_VDWb,293.15,(2000/25),25)

Eq_P_mass=VDWP(Eq_VDWa,Eq_VDWb,293.15,(2000/25),25)-VDWP(Eq_VDWa,Eq_VDWb,293.15,(1000/12),12)
## Show P difference
print(P_B-P_A)

## Section 0.2: Heat quantification

Enthalpy_A=((A_y_N2*A_N2_Cp)+(A_y_O2*A_O2_Cp))*A_n*A_T

Enthalpy_B=((B_y_N2*B_N2_Cp)+(B_y_O2*B_O2_Cp)+(B_y_Cl2*B_Cl2_Cp))*B_n*B_T

print(' ')
print('The total system enthalpy of Cube A is:')
print(Enthalpy_A)
print(' ')
print('the total system enthalpy of Cube B is:')
print(Enthalpy_B)
print(' ')

## From Solving: (738.06816+(-31.63228)*t+26.45148*t^2+191.06777t^3+((-0.18665)/(t^2))=734.5)
## Real Heat capacities are p gross.
T_Equilibrium=308.9 ##Kelvin

## heat transfer Coefficient calculation
## Calculation is treating the interface between the hot, chlorinated "air" and
## the cold, clean "air" - the N2, O2 mix
##plate length
L_bar=(4*(1/4)) ## 4 by 1 m^2 divided by... 4 meter perimeter

h_bar=1.32*((20/L_bar)**0.25)

## Section 0.3: Equilibrioum conditions
eq_y_Cl2=(1/25)
eq_y_N2=(20/25)
eq_y_O2=(4/25)
eq_y_A=(12/25)

## Section 0.4: Chemical Potential Calculations

## function defined to approximate the change in chemical potential with temperature
## NOTA BENE: this does not include pressure because it is at a low point, 
## Such an assumption will likely result in mild error with the potential 
## starting at a slightly positive value
def Mu(mu0, alpha, Delta_T):
    return mu0+(alpha*Delta_T)

Delta_T_A=A_T-273.15

Delta_T_B=B_T-273.15

Delta_T_Eq=T_Equilibrium-273.15

mu0_N2=0
alpha_N2=-191.5

mu0_O2=0
alpha_O2=-205.02

mu0_Cl2=0
alpha_Cl2=-222.96

## Initial potentials 

A_N2_Chemical_Potential=Mu(mu0_N2, alpha_N2, Delta_T_A)

A_O2_Chemical_Potential=Mu(mu0_O2,alpha_O2, Delta_T_A)

B_N2_Chemical_Potential=Mu(mu0_N2, alpha_N2, Delta_T_B)

B_O2_Chemical_Potential=Mu(mu0_O2,alpha_O2, Delta_T_B)

B_Cl2_Chemical_Potential=Mu(mu0_Cl2,alpha_Cl2, Delta_T_B)

## Equilibrium potentials 

A_N2_Chemical_Potential_Eq=Mu(mu0_N2, alpha_N2, Delta_T_Eq)

A_O2_Chemical_Potential_Eq=Mu(mu0_O2,alpha_O2, Delta_T_Eq)

B_N2_Chemical_Potential_Eq=Mu(mu0_N2, alpha_N2, Delta_T_Eq)

B_O2_Chemical_Potential_Eq=Mu(mu0_O2,alpha_O2, Delta_T_Eq)

B_Cl2_Chemical_Potential_Eq=Mu(mu0_Cl2,alpha_Cl2, Delta_T_Eq)

## Calculations focus only on the Chlorine, for time and illustration purposes
## this is, in part because Cl2 is the only unique species undergoing mass 
## diffusion

## the part where PDE's would come in, are herein. for my first iteration,
## I will show an approximation for the del(mu)/del(x) term in the soret effect
## the primary justification for this is theat the system is low pressure, and
## High Volume in a medium temp range. These conditions mean that the Van Der
## Waal's EOS used initially can reasonably be assumed to account for most of
## the deviation from ideality of the real gas

## Approximation 1: linear approximation/step change 
del_mu_Cl2=B_Cl2_Chemical_Potential-B_Cl2_Chemical_Potential_Eq
del_x_Cl2=eq_y_Cl2-A_y_Cl2

## Approximation 2: Exponential relation to time and derivative.

## Approximation 3: erf function 

##prepare for use in soret
Soret_potential_coefficient=del_mu_Cl2/del_x_Cl2

## Section 1.0: Soret effect Calculation

## define full system molar volume
V_m_tot=2000/25

## define Ideal gas molar Enthalpy 

Cp_ideal_Cl2=33.91

Cp_ideal_N2=29.125

Cp_ideal_O2=29.355

hm_ideal_Cl2=Cp_ideal_Cl2*B_T

hm_ideal_A=((A_y_N2*Cp_ideal_N2)+(A_y_O2*Cp_ideal_O2))*A_T

## Real molar Enthalpies

hm_Cl2=B_Cl2_Cp*B_T

hm_A_real=((A_y_N2*A_N2_Cp)+(A_y_O2*A_O2_Cp))*A_T

Vm_A=1000/12
Vm_Cl2=1000/13

Thermal_Diffusion_Coefficient_Soret=((Vm_A*Vm_Cl2)/((Vm_Cl2*B_y_Cl2)+(Vm_A*eq_y_A)))*(((hm_A_real/Vm_A)-(hm_Cl2/Vm_Cl2))/(eq_y_Cl2*Soret_potential_coefficient))

Soret_Calc=(-y_Cl2*(1-y_Cl2))*Thermal_Diffusion_Coefficient_Soret*((T_Equilibrium-B_T)/T_Equilibrium)

## Indication: the Soret effect dominates for ~10% of the total diffusion 
## this is compared tot he quilibrium endpoint of the Cl concentration

## Newton's Law of cooling answer
C_mixture=(y_Cl2*B_Cl2_Cp)+(y_O2*B_O2_Cp)+(y_N2*B_N2_Cp)

t_heating=(C_mixture/(1*h_bar))*np.log((B_T-A_T)/(B_T-T_Equilibrium))