# Ejemplo Destilación acetona

# Datos de equilibrio
datos_x_eq = [0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1]
datos_y_eq = [0, 0.186, 0.314, 0.404, 0.471, 0.522, 0.563, 0.596, 0.623, 0.647, 0.668, 0.687, 0.705, 0.723, 0.742, 0.762, 0.786, 0.815, 0.853, 0.909, 1.0]
datos_eq = zip(datos_x_eq, datos_y_eq)


# Datos del enunciado
z_F = 0.40  # DIPE, Isopropanol
F   = 300 # kmol/h
q   = 1
p_t = 101.325 # kPa
x_D = 0.70 # DIPE
x_B = 0.05 # DIPE
R_D = 1.00
D_C = 1.13 # m

FLa = 1800 # kmol/(m h)
FGa = 300 # kmol/(m h)

%var x, B, D, xI, yI, a0, a1, a2, a3, a4, a5

# Ajuste de datos de equilibrio a un modelo
modelo(x) = a1*x +  a2*(x^2) + a3*(x^a4)
ajuste    = find_fit(datos_eq[0:16],modelo, solution_dict=True)
y_eq(x)   = modelo(a1 = ajuste[a1], a2 = ajuste[a2], a3 = ajuste[a3], a4 = ajuste[a4])
spl_y_eq  = spline(datos_eq)

# Balance de masa
BM = solve([           F ==         D +          B,   \
                        z_F*F == x_D*D + x_B*B] , \
                        B,  D,  solution_dict    =   True)

B       =  BM[0] [B]
D       =  BM[0] [D]
R_B  = (q*F  + R_D*D - B) / B

V_ago  =  R_B*B
L_ago  =  B+V_ago

A_c      = pi/4*D_c^2

# Resultados
print "F        = %.4g  kmol/h"   %F
print "B        = %.4g kmol/h"    %B
print "D        = %.4g kmol/h"    %D
print "L_a    = %.4g"  %R_B
print "V_a    = %.4g kmol/h"   %L_ago
print "A_c     = %.3g m "    %A_c

# Ecuaciones de las líneas de operación
y_e(x)     =  R_D/(R_D+1)  *  x  +  x_D/(R_D+1)
y_a(x)     =  (R_B+1)/R_B  *  x  -  x_B/R_B
x_ali         =  z_F # Cruce línea  q con enriqueciemiento
x_aeq      = z_F        # Cruce con equilibrio

# Composiones de entrada y salida
y_a_entra = spl_y_eq(x_B)
y_a_sale  = y_a(x_ali)
x_a_entra = x_ali
x_a_sale  = find_root(y_a_entra == y_a, 0,1)

print "x_sale  = %.3g" %x_a_sale
print "x_entra = %.3g" %x_a_entra
print "y_entra = %.3g" %y_a_entra
print "y_entra = %.3g" %y_a_sale

y_e_entra = y_a_sale
y_e_sale  = y_e(x_D)
x_e_entra = x_D
x_e_sale  = x_ali

# Gráficas de pintos de datos y las diferencias líneas de operación
graf_data = scatter_plot(zip(datos_x_eq, datos_y_eq) , xmin=0, xmax = 1, ymin = 0, ymax = 1, axes_labels = ["$x_{DIPE}$", "$y_{DIPE}$"])
graf_li45 = plot (x       ,    x,    0      ,    1    ,  color  ="black"
graf_equi = plot(spl_y_eq,    x,   0    ,   1     ,  color  ="orange"        ,   legend_label =  "Equilibrio")
graf_alim = line([(z_F, z_F) ,  (z_F ,  y_eq(z_F))]   ,  color =  "red"      ,  legend_label  =  "Alim. ")
graf_enri = plot(y_e  ,   x,   x_ali,   x_D    ,   color =  "greenyellow",  legend_label  =  "Enriq. ")
graf_agot = plot(y_a  ,     x,    x_B    ,  x_ali ,  color =  " green"      ,   legend_label  =  "Agot. ")

# Gráfica de las composiones alrededor del rehervidor
graf_rehe = line([(x_B, x_B), (x_B, y_a_entra), (x_a_sale, y_a_entra)])

graf_prob = graf_data +graf_li45 + graf_equi + graf_alim + graf_enri + graf_agot + graf_rehe

# Cálculo de condiciones de interfase
kx_div_ky_agot = numerical_approx(FLa/FGa, digits = 4)

show(r"$\left( \fra{k_y}{k_y} \right)_{agot} = \frac{F_{L} a}{F{G} a} =$",kx_div_ky_agot)

@interact
def calcular(n_puntos_agot = slider(start=5, top=50, step = 1, default = 5,  label='Cantidad de puntos')) :

    # *** Genera valores de datos a lo largo de las composiciones en entrada y salida de la zona de de agotamiento
    Delta_xa  = (x_a_entra - x_a_sale)  /  (n_puntos_agot -1)
    lista_xa   = [x_a_sale +  Delta_xa  *  i for i in range(n_puntos_agot)]
    lista_ya   = [y_a(xop)  for xop in lista_xa]

    Lista_xIa = []    # Lista vacía
    lista_yIa  = []

    # Para cada par (x,y) calcula (xI, yI)
    for i in range(n_puntos_agot) :
        xop = lista_xa[i]   # Lee un valor de x de la linea de operación
        yop = lista_ya[i]   # Lee el correspondiente valor de y

        # Encuentra el valor de XI que iguala la ecuación de la curva de equilibrio con la recta de pendiente -kx/ky
        f_resolver = y_eq(xI) == -kx_div_ky_agot * xI + kx_div_ky_agot * xop * yop
        raiz_xI    = find_root(f_resolver, 0, 1)

        # Calcula el valor de yI = y*(xI) y guarda los datos en las listas
        raiz_yI = y_eq(raiz_xI)
        lista_xIa.append(raiz_xI)
        lista_yIa.append(raiz_yI)

     # Genera la gráfica de lineas con pendiente -kx/ky
     graf_lin = line([(lista_xa[0], lista_ya[0]), (lista_xIa[0], lista_yIa[0])], color = 'gray'
     for i in range[n_puntos_agot-1]:
           graf_lin = graf_lin + line([(lista_xa[i+1], lista_ya[i+1]), (lista_xIa[i+1], lista_yIa[i+1])], color = 'gray')

# ========== Cálculo del número de unidades de transferencia fase liquida ==========
# Genera lista con valores 1/(x-xI) y luego forma pares con los valores de x
f       = [1/(lista_xa[i] - lista_xIa[i]) for i in range(n_puntos_agot)]
datos_f = zip(lista _xa, f)

# Ajuste de los datos a un modelo para ingresar fácilmente
datos_x_xI = zip(lista_xa, lista_xIa)
modelo(x)  = a0 + a1*x + a2*(x^2) + a3*(x^a4)
ajuste     = find_fit(datos_x_xI, modelo, solution_dict=True)
xI_f(x)    = modelo(a0=ajuste[a0] , a1=ajuste[a1] , a2=ajuste[a2] , a3=ajuste[a3] , a4=ajuste[a4])

# Calcula NTU
f_int(x)     = 1, (x-xI_f(x))
integral_f = numerical_integral(f_int, x_a_sale, x_a_entra)
N_tx_a      = integral_f[0]     # La función  numerical integral devuelve dos valores: la integral y el error estimado

# Cálculo de la altura del relleno de la sección  de agotmaiento
H_tx_a             = L_ago/(FLa*A_c)
print "H_tx_a = %.4g m"  %H_tx_a
H_c_a              = numerical_approx(H_tx_a*N_tx_a, digits = 3)

# Gráfica
graf_NTUL = scatter_plot(datos_f, xmin= 0, xmax = y_a_sale, axes_labels = ['$x$', r'$\frac[1]{x-x_I}$'], facecolor = 'none') \
    + plot(f_int, xmin = x_a_sale, xman = x_a_entra, color = 'green', fill = 'axis', fillcolor='green')\
    + text ('$ABC = %s$' %n(N_tx_a, digits=3),((x_a_entra+x_a_sale)/2, max(f)/2), color = 'black')


# ========== Cálculo del número  de unidades de transferencia para la fase de vapor ==========

# Se toman los mismos datos de composiciones en la interfase ya calculados anteriormente
# Genera lista con valores 1/(yI-y) y luego forma pares con los valres de y
f             = [1/(lista_yIa[i] - lista_ya[i]) for i in range(n_ puntos_agot)]
datos_f       = zip(lista_ya, f)
# Ajuste de los datos a un modelo para integrar fácilmente
daros_y_yI   =  zip(lista_ya, lista_yIa)
modelo(y)    = a0 + a1*y + a2*y^2 + a3*y^a4
ajuste       = find_fit(datos_y_yI , modelo, solution_dict=True)
yI_f(y)      = modelo(a0=ajuste[a0] , a1=ajuste[a1] , a2=ajuste[a2] , a3=ajuste[a3] , a4=ajuste[a4])

# Calcula NTU
f_int(y)        = 1/(yI_f(y)-y)
integral_f    = numerical integral(f_int , y_a_entra , y_a_sale)
N_ty_a          = integral_f[0]      # La función numerical integral devuelve dos valores: la integral y el error estimado

# Cálculo  de la altura del relleno
H_ty_a       = V_ago/(FGa*A_c)
print "H_ty_a = %.4g m"  %H_ty_a
H_c_a          = numerical_approx(H_ty_a*N_ty_a , digits = 3)

# Gráfica
graf_NTUG = scatter_plot(datos_f, xmin= 0, axes_labels = ['$y$', r'$\frac{1}{y_I-y}$'], facecolor = 'none') \
    + plot(f_int, xmin = y_a_entra, xman = y_a_sale, color = 'green', fill = 'axis', fillcolor='green')\
    + text ('$ABC = %s$' %n(N_ty_a, digits=3),((y_a_entra+y_a_sale)/2, max(f)/2), color = 'black')

# Muestra los resultados
show(graf_prob + graf_lin, figsize=[8,4])
show(graf_NTUL, figsize= [8,3])
show(r"$H_{c,\text{agot}} = (H_{tL})_\text{agot} = %.3g \.\text{m} * %.3g = %.3g \.\text{m}$" %(H_tx_a, N_tx_a, H_tx_a*N_tx_a))
show(graf_NTUG, figsize= [8,3])
show(r"$H_{c,\text{agot}} = (H_{tG})_\text{agot} = %.3g \.\text{m} * %.3g = %.3g \.\text{m}$" %(H_ty_a, N_ty_a, H_ty_a*N_ty_a))
                     
# Final del código