---
# System prepended metadata

title: Untitled

---

File: 
```python
# -*- coding: utf-8 -*-
"""
Created on Sun, April 17 2022
This script is written to analytically and semi-analytical-emperically fitting to solve overall kGAv for spray system that absorbs CO2 from the surrounding air/flue gas 
@author: Awan
"""
import math
import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import *
import pandas as pd
from Input_Params import *
from Side_Functions import *


plt.close('all')
WritetoFolder='Results'



# def KGAv_Theo(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading):
#     hm, U, d = hm_calc(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading)
#     A_1d = 4*3.14*d*d/4
#     Vol_1d = (4/3)*3.14*(d**3)/8
#     Q_L = Q_l_flux*A_D/3600         # Liquid flow rate [m3/s] when Q_l_flux is in m3/m2-h
#     Q_G =Q_g_flux*A_D/3600          # Gas Flow Rate [m3/s] (For aroonwilas paper, the flow rate is assumed total gas flow rate for now as info is not clerarly mentioned in paper)
#     Q_G_mole = 3600*(Q_G*1000/(22.4*1000))*((T/298)/(P_tot/101.325))/(A_D)    #Gas Flow Rate in kmol/m2-hr

#     X2_theo = X1*math.exp(-hm*Q_L*A_1d*H/(Q_G*Vol_1d*U))
#     Y2_theo = X2_theo/(1-X2_theo)
#     KGAv_theo = (Q_G_mole*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_theo) + Y1 - Y2_theo)
#     cap_eff = (X1-X2_theo)*100/X1                ### Change the definition as per paper that it is being validated to, if needed
#     return KGAv_theo, cap_eff

X = np.zeros(2)         #X[0] is loading and X[1] is X_CO2

# def loading_func(z, CO2_loading):
#     hm, U, d = hm_calc(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading)
#     AB = (P_tot/(C_init*Rbar*T))*(6/(d*U))*hm
#     return AB*P_CO2_in/P_tot

# def X_CO2_func(z,X_CO2):
#     hm, U, d = hm_calc(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading)
#     B = (Q_l_flux/Q_g_flux)*(6/(d*U))*hm
#     return B*X_CO2


def X_func(z, X):
    hm, d = hm_calc(Q_l_flux, Q_g_flux, X[1]*P_tot, T, X[0]/(1-X[0]), U)
    AB = (P_tot/(C_init*Rbar*T))*(6/(d*U))*hm
    B = (Q_l_flux/Q_g_flux)*(6/(d*U))*hm
    print(AB*X[1], B*X[1])
    return AB*X[1], B*X[1]


npoints = len(x_exp)
Inputs = np.zeros([len(x_exp),4])
KGAv_Theo_all = np.zeros(npoints)
cap_eff_all = np.zeros(npoints)
KGAv_Fit_all = np.zeros(npoints)
cap_eff_Fit_all = np.zeros(npoints)


# def X2_calc_integrand():#, Q_l_flux, Q_g_flux, P_CO2_in, T):
#     hm, U, d = hm_calc(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading_bottomUp)
#     A_1d = 4*3.14*d*d/4
#     Vol_1d = (4/3)*3.14*(d**3)/8
#     Q_L = Q_l_flux*A_D/3600         # Liquid flow rate [m3/s] when Q_l_flux is in m3/m2-h
#     Q_G =Q_g_flux*A_D/3600          # Gas Flow Rate [m3/s] (For aroonwilas paper, the flow rate is assumed total gas flow rate for now as info is not clerarly mentioned in paper)

#     return -hm*Q_L*A_1d/(Q_G*Vol_1d*U)

for i in range(npoints):
    if Var == 'Liquid Flow Rate':
        Q_l_flux = x_exp[i]
        Inputs[i] = [Q_l_flux, Q_g_flux, P_CO2_in, T]
    if Var == 'Gas Flow Rate':
        Q_g_flux = x_exp[i]
        Inputs[i] = [Q_l_flux, Q_g_flux, P_CO2_in, T]
    
    U = U_calc2(Q_l_flux, Q_g_flux, T)
    
    print('Starting X_sol')
    # loading_sol = scipy.integrate.solve_ivp(loading_func, [0, H], [0])
    # X_CO2_sol = scipy.integrate.solve_ivp(X_CO2_func, [0, H], [0.01])
    
    X_sol = scipy.integrate.solve_ivp(X_func, (0, H), [0, 0.1])
    
    # CO2_loading_bottomUp = loading_sol.y[::-1]
    # z_array_bottomUp = loading_sol.t[::-1]
    # z_array_bottomUp[:] = H-z_array_bottomUp[:] 
    
   # hm, U, d = hm_calc(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading_bottomUp)
    # ln_X2_array = X2_calc_integrand()
    # Integrand = scipy.integrate.simps(ln_X2_array, z_array_bottomUp)
    # print (Integrand)
    # X2 = X1*math.exp(Integrand)
    #lnX2_by_X1 = scipy.integrate.quad(X2_calc_integrand, 0, H)
    
    #print(X2)
    
    # for j in range(len(loading_sol.y)):
    #     hm_array = 
    # KGAv_theo, cap_eff = KGAv_Theo(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading)
    # KGAv_Theo_all[i] = KGAv_theo
    # cap_eff_all[i] = cap_eff


# def KGAv_Fit(Inputs, K, c):
    
#     for i in range(npoints):
#         [Q_l_flux, Q_g_flux, P_CO2_in, T] = Inputs[i]
#         hm, U, d = hm_calc(Q_l_flux, Q_g_flux, P_CO2_in, T)
#         A_1d = 4*3.14*d*d/4
#         Vol_1d = (4/3)*3.14*(d**3)/8
#         Q_L = Q_l_flux*A_D/3600         # Liquid flow rate [m3/s] when Q_l_flux is in m3/m2-h
#         Q_G =Q_g_flux*A_D/3600          # Gas Flow Rate [m3/s] (For aroonwilas paper, the flow rate is assumed total gas flow rate for now as info is not clerarly mentioned in paper)
#         Q_G_mole = 3600*(Q_G*1000/(22.4*1000))*((T/298)/(P_tot/101.325))/(A_D)    #Gas Flow Rate in kmol/m3-hr
        
#         X2_Fit = X1*math.exp(-hm*Q_L*A_1d*H*K/(Q_G*Vol_1d*U))
#         Y2_Fit = X2_Fit/(1-X2_Fit)
#         KGAv_Fit = (Q_G_mole*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_Fit) + Y1 - Y2_Fit) + c
#         KGAv_Fit_all[i] = KGAv_Fit
#         cap_eff_Fit = (X1-X2_Fit)*100/X1                ### Change the definition as per paper that it is being validated to, if needed
#         cap_eff_Fit_all[i] = cap_eff_Fit

#     return KGAv_Fit_all


# params, extras = curve_fit(KGAv_Fit, Inputs, y_exp)
# K = params[0]
# c = params[1]
# print(K,c)
# KGAv_Fit_all = KGAv_Fit(Inputs, K, c)

plt.figure()
plt.plot(x_exp, KGAv_Theo_all, 'go-', label = 'Analytical Model')
# plt.plot(x_exp, KGAv_Fit_all, 'ro-', label = 'Model with Fitting')
plt.plot(x_exp, y_exp, 'bo-', label = 'Expt Data')
plt.legend()










# if Var == 'Liquid Flow Rate':
#     Q_l_flux_array = np.zeros(len(x_exp))           
#     Q_l_flux_array[:] = x_exp[:]
    
#     Q_g_flux = 382      #m3/m2-h
#     Q_G =Q_g_flux*A_D/3600   #Gas Flow Rate [m3/s] (Assumed total gas flow rate for now as info is not clerarly mentioned in paper)
#     Q_G_mole = 3600*(Q_G*1000/(22.4*1000))*(((273+T)/298)/(P_tot/101.325))/(A_D)    #Gas Flow Rate in kmol/m3-hr
    
#     hm_array, U_array, d_array = hm_calc_QLarray(Q_l_flux_array, Q_g_flux, A_D, K_factor, d_orifice, P_CO2_in, P_tot, T, nozzle_type)
    
    
#     def KGAv_Theo_for_Ql(Q_l_flux_array, hm_array, U_array, d_array):
#         Q_L_array = np.zeros(len(Q_l_flux_array))
#         X2_theo = np.zeros(len(Q_l_flux_array))
#         Y2_theo = np.zeros(len(Q_l_flux_array))
#         KGAv_theo = np.zeros(len(Q_l_flux_array))
        
#         for i in range(len(Q_l_flux_array)):
#             Q_L_array[i] = Q_l_flux_array[i]*A_D/3600
            
#             d = d_array[i]
#             A_1d = 4*3.14*d*d/4
#             Vol_1d = (4/3)*3.14*(d**3)/8
            
#             X2_theo[i] = X1*math.exp(-hm_array[i]*Q_L_array[i]*A_1d*H/(Q_G*Vol_1d*U_array[i]))
#             #X2_calc[i] = X1*math.exp(-c*H*((Q_L_array[i]*d_orifice)**.5)/Q_G)
#             Y2_theo[i] = X2_theo[i]/(1-X2_theo[i])
#             KGAv_theo[i] = (Q_G_mole*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_theo[i])+Y1-Y2_theo[i])
            
#         return KGAv_theo


# if Var == 'Gas Flow Rate':
#     Q_g_flux_array = np.zeros(len(x2_exp))
#     Q_g_flux_array[:] = x2_exp[:]
    
#     Q_L = Q_l_flux*A_D/3600
    
#     hm_array, U_array, d_array = hm_calc_QGarray(Q_g_flux_array, Q_l_flux, A_D, K_factor, d_orifice, P_CO2_in, P_tot, T, nozzle_type)

#     def KGAv_Theo_for_Qg(Q_g_flux_array, hm_array, U_array, d_array):
#         Q_G_array = np.zeros(len(Q_g_flux_array))
#         Q_G_mole_array = np.zeros(len(Q_g_flux_array))
#         X2_theo = np.zeros(len(Q_g_flux_array))
#         Y2_theo = np.zeros(len(Q_g_flux_array))
#         KGAv_theo = np.zeros(len(Q_g_flux_array))
        
#         for i in range(len(Q_g_flux_array)):
#             Q_G_array[i] = Q_g_flux_array[i]*A_D/3600
#             Q_G_mole_array[i] = 3600*(Q_G_array[i]*1000/(22.4*1000))*(((273+T)/298)/(P_tot/101.325))/(A_D)
            
#             d = d_array[i]
#             A_1d = 4*3.14*d*d/4
#             Vol_1d = (4/3)*3.14*(d**3)/8
            
#             X2_theo[i] = X1*math.exp(-hm_array[i]*Q_L*A_1d*H/(Q_G_array[i]*Vol_1d*U_array[i]))
#             #X2_calc[i] = X1*math.exp(-c*H*((Q_L_array[i]*d_orifice)**.5)/Q_G)
#             Y2_theo[i] = X2_theo[i]/(1-X2_theo[i])
#             KGAv_theo[i] = (Q_G_mole_array[i]*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_theo[i])+Y1-Y2_theo[i])
            
#         return KGAv_theo


# def KGAv_Calculator_for_Ql(Q_l_flux_array, c, K):
#     Q_L_array = np.zeros(len(Q_l_flux_array))
#     X2_calc = np.zeros(len(Q_l_flux_array))
#     Y2_calc = np.zeros(len(Q_l_flux_array))
#     KGAv_calc = np.zeros(len(Q_l_flux_array))
    
#     for i in range(len(Q_l_flux_array)):
#         Q_L_array[i] = Q_l_flux_array[i]*A_D/3600 
#         X2_calc[i] = X1*math.exp(-c*hm_array[i]*H*((Q_L_array[i]*d_orifice)**.5)/Q_G)
#         #X2_calc[i] = X1*math.exp(-c*H*((Q_L_array[i]*d_orifice)**.5)/Q_G)
#         Y2_calc[i] = X2_calc[i]/(1-X2_calc[i])
#         #KGAv_calc[i] = (Q_G_mole*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_calc[i])+Y1-Y2_calc[i])
#         KGAv_calc[i] = (Q_G_mole*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_calc[i])+Y1-Y2_calc[i]) + K
   
#     return KGAv_calc

# #def KGAv_Calculator_for_Qg(Q_g_flux_array, c, K):
# def KGAv_Calculator_for_Qg(Q_g_flux_array, c, K, c2):
#     Q_G_array = np.zeros(len(Q_g_flux_array))
#     X2_calc = np.zeros(len(Q_g_flux_array))
#     Y2_calc = np.zeros(len(Q_g_flux_array))
#     KGAv_calc = np.zeros(len(Q_g_flux_array))
    
#     for i in range(len(Q_g_flux_array)):
#         Q_G_array[i] = Q_g_flux_array[i]*A_D/3600 
#         X2_calc[i] = X1*math.exp(-c*hm_array[i]*H*((Q_L*d_orifice)**.5)/Q_G_array[i]**c2)
#         #X2_calc[i] = X1*math.exp(-c*6*hm_array[i]*H*Q_L/(d_array[i]*Q_G_array[i]*U_array[i]))

#         Y2_calc[i] = X2_calc[i]/(1-X2_calc[i])
#         print(Y2_calc[i])
#         Q_G_mole = 3600*(Q_G_array[i]*1000/(22.4*1000))*(((273+T)/298)/(P_tot/101.325))/(A_D)
#         #KGAv_calc[i] = (Q_G_mole*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_calc[i])+Y1-Y2_calc[i])
#         KGAv_calc[i] = (Q_G_mole*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_calc[i])+Y1-Y2_calc[i]) + K
   
#     return KGAv_calc

# if Var == 'Liquid Flow Rate':
#     params, extras = curve_fit(KGAv_Calculator_for_Ql, x_exp, y_exp)
#     c = params[0]
#     K = params[1]
    
#     plt.plot(x_exp, KGAv_Theo_for_Ql(Q_l_flux_array, hm_array, U_array, d_array), 'go-', label = 'Analytical Model')
#     plt.plot(x_exp, KGAv_Calculator_for_Ql(Q_l_flux_array, c, K), 'ro-', label = 'Semi-Analytical-Empirical Fit')
#     plt.plot(x_exp, y_exp, 'bo-', label = 'Experimental data')
#     plt.xlabel('Liquid flow rate (m3/m2-h)')


# if Var == 'Gas Flow Rate':
#     params, extras = curve_fit(KGAv_Calculator_for_Qg, x2_exp, y_exp)
#     c = params[0]
#     K = params[1]
#     c2 = params[2]
    
#     #plt.plot(x2_exp, KGAv_Theo_for_Qg(Q_g_flux_array, hm_array, U_array, d_array), 'g', label = 'Theoretical Model')
#     plt.plot(x2_exp, KGAv_Calculator_for_Qg(Q_g_flux_array, c, K, c2), 'ro-', label = 'Analytical Fit')
#     #plt.plot(x2_exp, KGAv_Calculator_for_Qg(Q_g_flux_array, c, K), 'ro', label = 'Analytical Fit')
#     plt.plot(x2_exp, y_exp, 'bo-', label = 'Experimental data')
#     plt.xlabel('Gas flow rate (m3/m2-h)')

# plt.ylabel('KgAe (kmol/m3-h-kPa)')
# plt.legend(loc='best')
# #plt.title('Kuntz fit 2f and 2c')
# #plt.savefig()#'%s/Concentration_timehistory.pdf'%WritetoFolder, bbox_inches='tight') 



# if nozzle_type == 'P-40':
#     c_P40 = c
#     K_P40 = K
# if nozzle_type == 'P-28':
#     c_P28 = c
#     K_P28 = K


# print(c, K)

######plt.plot(x_exp, KGAv_Calculator_for_Ql(x_exp, c), 'r')
# plt.plot(x_exp, KGAv_Calculator_for_Ql(x_exp, c, K), 'ro', label = 'Analytical Fit')
# #####plt.plot(x_exp, KGAv_Calculator_for_Ql(x_exp, 504.1583049223172, -1.6742542161512164), 'g')
# #####plt.plot(x_exp, KGAv_Calculator_for_Ql(x_exp, 246.1363698486537, -0.17978853997655872), 'g')
# ######plt.plot(x_exp, KGAv_Calculator_for_Ql(x_exp, 1737430449.8144524, 1.39), 'g')

# plt.plot(x_exp, y_exp, 'bo', label = 'Experimental data')

#plt.plot(x2_exp, KGAv_Calculator_for_Qg(x2_exp, c,c2, K), 'ro', label = 'Analytical Fit')
#plt.plot(x2_exp, y_exp, 'bo', label = 'Experimental data')


#



# r_probe=R/8
# t_probe=1
# C_probe=DropletDiffusion_Functions.calc_probe(r_probe,t_probe,R,Rbar,kg,Bi,Henry,DCO2,Cinit,PCO2)

# tstart=1
# tfinal=100
# Ntime=20
# tsample=np.linspace(tstart,tfinal,Ntime)
# C_all_Rover4=np.zeros((Ntime,))
# C_all_Rover2=np.zeros((Ntime,))
# C_all_Rover1=np.zeros((Ntime,))
# for i in range(0,Ntime):
#     C_all_Rover4[i]=DropletDiffusion_Functions.calc_probe(R/4,tsample[i],R,Rbar,kg,Bi,Henry,DCO2,Cinit,PCO2)
#     C_all_Rover2[i]=DropletDiffusion_Functions.calc_probe(R/2,tsample[i],R,Rbar,kg,Bi,Henry,DCO2,Cinit,PCO2)
#     C_all_Rover1[i]=DropletDiffusion_Functions.calc_probe(R/1,tsample[i],R,Rbar,kg,Bi,Henry,DCO2,Cinit,PCO2)


# Ngrid=20
# rstart=R/Ngrid
# rsample=np.linspace(rstart,R,Ngrid)
# C_all_t10=np.zeros((Ngrid,))
# C_all_t20=np.zeros((Ngrid,))
# C_all_t30=np.zeros((Ngrid,))
# for i in range(0,Ngrid):
#     C_all_t10[i]=DropletDiffusion_Functions.calc_probe(rsample[i],10,R,Rbar,kg,Bi,Henry,DCO2,Cinit,PCO2)
#     C_all_t20[i]=DropletDiffusion_Functions.calc_probe(rsample[i],20,R,Rbar,kg,Bi,Henry,DCO2,Cinit,PCO2)
#     C_all_t30[i]=DropletDiffusion_Functions.calc_probe(rsample[i],30,R,Rbar,kg,Bi,Henry,DCO2,Cinit,PCO2)


# plt.figure()
# plt.plot(tsample,C_all_Rover4,label='R/4')
# plt.plot(tsample,C_all_Rover2,label='R/2')
# plt.plot(tsample,C_all_Rover1,label='R')
# plt.xlabel('time (s)')
# plt.ylabel('Concentration (mol/m3)')
# plt.legend(loc='best')
# plt.title('Time history of concentration at various locations of the droplet')
# # plt.savefig('%s/Concentration_timehistory.pdf'%WritetoFolder, bbox_inches='tight') 

# plt.figure()
# plt.plot(rsample,C_all_t10,label='t=10s')
# plt.plot(rsample,C_all_t20,label='t=20s')
# plt.plot(rsample,C_all_t30,label='t=30s')
# plt.xlabel('r (m)')
# plt.ylabel('Concentration (mol/m3)')
# plt.legend(loc='best')
# plt.title('Radial distribution of concentration in the droplet at various times')
# # plt.savefig('%s/Concentration_timehistory.pdf'%WritetoFolder, bbox_inches='tight') 
```


```python
# -*- coding: utf-8 -*-
"""
Created on Fri Feb  4 17:50:33 2022

@author: ab65462
"""

#### Includes design parameters, droplet size, add more..

import os
import math
import numpy as np
import pandas as pd
import scipy.interpolate 

### Universal Constants-
g = 9.8                               # Gravitational acceleration [m/s2]
Rbar = 8.3145                           # Universal gas constant [J/mol/K]

### Pressure and Temperature at Inlet and ambient/operating conditions

P_tot = 101.325                      #Ambient Pressure [kPa]
T = 25 + 273.15                     #Inlet Temperature [K]
P_CO2_in = 15.000                    #Inlet CO2 partial pressure [kPa]
X1= P_CO2_in/P_tot
Y1 = X1/(1-X1)

C_init = 5                              #Initial Concentration in [M] or [kmol/m3]

###Reactor/Channel Design Specs-
D = 0.1          #Inner Diameter of Channel [m]
H = 0.55         #Height of channel [m]
A_D = np.pi*D*D/4   #Cross sectional area of channel


#nozzle_type = 'P-40'
#nozzle_type = 'P-28'
#nozzle_type = 'MPL30N'

folder = r"C:\Users\ab65462\Desktop\CO2 Musk Foundation\Spray Model\Kuntz_Aroonwilas_Data"

#fig_validate = '2b-1.91.csv'
#fig_validate = '2b-4.59.csv'
#fig_validate = '2fc-P28.csv'
fig_validate = '2f-P40.csv'

filename = os.path.join(folder,fig_validate)

df = pd.read_csv(filename)
kG_v_Q = df.to_numpy()
x_exp = np.array(kG_v_Q[:,0])
y_exp = np.array(kG_v_Q[:,1])


### Need to create filename for df so that these if conditions become cleaner

if fig_validate == '2b-1.91.csv':
    nozzle_type = 'P-28'
    Q_l_flux = 1.91                         # Flux of liquid flow rate in m3/m2-h
    Var = 'Gas Flow Rate'

if fig_validate == '2b-4.59.csv':
    nozzle_type = 'P-28'
    Q_l_flux = 4.59                         # Flux of liquid flow rate in m3/m2-h
    Var = 'Gas Flow Rate'

if fig_validate == '2fc-P28.csv':
    nozzle_type = 'P-28'
    Q_g_flux = 382                         # Flux of gas flow rate in m3/m2-h
    Var = 'Liquid Flow Rate'

if fig_validate == '2f-P40.csv':
    nozzle_type = 'P-40'
    Q_g_flux = 382                         # Flux of gas flow rate in m3/m2-h
    Var = 'Liquid Flow Rate'

#nozzle_type = 'MPL30N'

####Nozzle Specs-
if nozzle_type == 'P-40':
    df_Q_vs_delP = pd.read_csv(r'C:\Users\ab65462\Desktop\CO2 Musk Foundation\Spray Model\Manufacturer Data\P-40-gpm_vs_psi.csv')
    df_d_vs_delP = pd.read_csv(r'C:\Users\ab65462\Desktop\CO2 Musk Foundation\Spray Model\Manufacturer Data\P-40-micron_vs_psi.csv')
    QvP = df_Q_vs_delP.to_numpy()
    dvP = df_d_vs_delP.to_numpy()
    d_orifice = 0.042*0.254     #[m]
    K_factor = .0443            #Manufacturer's design parameter for nozzle that relates deltaP and GPM


if nozzle_type == 'P-28':
    df_Q_vs_delP = pd.read_csv(r'C:\Users\ab65462\Desktop\CO2 Musk Foundation\Spray Model\Manufacturer Data\P-28-gpm_vs_psi.csv')
    df_d_vs_delP = pd.read_csv(r'C:\Users\ab65462\Desktop\CO2 Musk Foundation\Spray Model\Manufacturer Data\P-28-micron_vs_psi.csv')
    QvP = df_Q_vs_delP.to_numpy()
    dvP = df_d_vs_delP.to_numpy()
    d_orifice = 0.028*0.254     #[m]
    K_factor = .0206             #Manufacturer's design parameter for nozzle that relates deltaP and GPM
    
if nozzle_type == 'MPL30N':
    df_Q_vs_delP = pd.read_csv(r'C:\Users\ab65462\Desktop\CO2 Musk Foundation\Spray Model\Manufacturer Data\MPL30N-gpm_vs_psi.csv')
    df_d_vs_delP = pd.read_csv(r'C:\Users\ab65462\Desktop\CO2 Musk Foundation\Spray Model\Manufacturer Data\MPL30N-micron_vs_psi.csv')
    QvP = df_Q_vs_delP.to_numpy()
    dvP = df_d_vs_delP.to_numpy()
    d_orifice = 1.1*1e-3     #[m] 
    #K_factor = ?             #Manufacturer's design parameter for nozzle that relates deltaP and GPM
    

def d_from_Q(Q_GPM):                #Returns droplet diameter in um from flowrate in gpm
    f_PvQ = scipy.interpolate.interp1d(QvP[:,1], QvP[:,0])
    f_dvP = scipy.interpolate.interp1d(dvP[:,0], dvP[:,1])
    delP = f_PvQ(Q_GPM)
    #print(delP, 'psi')
    droplet_dia = f_dvP(delP)
    return droplet_dia

#print(d_from_Q(.3))
```


```python
# -*- coding: utf-8 -*-
"""
Created on Wed Apr  6 21:05:44 2022

@author: ab65462
"""

### Includes simple side functions like viscosity, henry constant etc calculators


import numpy as np
import math
import matplotlib.pyplot as plt
import scipy.interpolate 
from Input_Params import *
from scipy.integrate import *


def D_and_H_CO2_l_calc(T):                   # Temp input is in K, returns D and H for CO2 in 5M MEA at 30C
    D_N2O_w = 5.07*1e-6*math.exp(-2371/T)    # m2/s ADDD REFERENCE FOR THESE if it is not Tsai etal 2000!!!!
    D_CO2_w = 2.35*1e-6*math.exp(-2119/T)    # m2/s
    D_N2O_MEA = 1.41*1e-9                    # 30C and 5M MEA, data from Jium-Jie Ko et al 2001, [m2/s]
    D_CO2_MEA = D_N2O_MEA*D_CO2_w/D_N2O_w
    
    H_N2O_w = 8.547*1e6*math.exp(-2284/T)   #[kPa-m3/kmol]
    H_CO2_w = 2.8249*1e6*math.exp(-2044/T)  #[kPa-m3/kmol]
    H_N2O_MEA = 4924                        #30C and 5M,  data from Tsai et al 2000, [kPa-m3/kmol]
    H_CO2_MEA = H_N2O_MEA*H_CO2_w/H_N2O_w
    ##### need to add a function of temp to H_N2O_MEA and D_N2O_MEA, then later add effect of concentration on all this too
    return D_CO2_MEA, H_CO2_MEA

def k2_calc(T):                             #Temp input is in K
    k2 = 1.72*1e11*math.exp(-4915.5/T)        #Ramezaini etal 2021, 2nd order reaction rate used here, [m3/mols]
    return k2

def calc_viscosity(T):          #T in Kelvin, and dynamic viscosity in [kg/m-s]
    #The following relation is based on Sutherland's formula for 3 coefficients (values are taken for air) ...ref- https://www.afs.enea.it/project/neptunius/docs/fluent/html/ug/node294.htm
    T0 = 273.11     #K
    mu_0 = 1.716*1e-5   #kg/m-s
    Suth = 110.56        #Sutherland Constant (K)
    mu = mu_0*((T/T0)**1.5)*((T0+Suth)/(T+Suth))
    return mu

#plt.figure()

def U_calc(d, Ujet, nu_gas, U_G, H, rho_gas, rho_drop):
    tau = (d**2)*(rho_drop/rho_gas)/(18*nu_gas)
    #print(tau)
    #print((U_G+Ujet)*d/nu_gas, 'Re')
    U0 = Ujet + U_G           # Adjusts the jet velocity to relative velocity wrt air (counterflow)
    n_points = 1001
    tf = 5
    t = np.linspace(0,tf,n_points)
    U = np.zeros(n_points)      # Absolute Velocity of fluid in m/s
    s = np.zeros(n_points)
    reached_column_end = 0

    for i in range(n_points):
        U[i] = U0*math.exp(-t[i]/tau) + g*tau*(1-math.exp(-t[i]/tau)) - U_G
        s[i] = U0*tau*(1-math.exp(-t[i]/tau)) + g*tau*t[i] - g*(tau**2)*(1-math.exp(-t[i]/tau)) - U_G*t[i]
        #print(t[i], s[i], U[i])
        

        if s[i] > H and reached_column_end == 0:
            print('Here!!')
            reached_column_end = 1
            n_points_needed = i+1
            residence_time = (t[i] +t[i-1])/2
            plt.plot(t[:], s[:], 'b-', label = 'U vs s')
            #print(s[i])
            #residence_time = (t[i-1])/2
            print(residence_time, 't_res')
            break
    
    n_points_needed = i+1
    print(n_points_needed)
    plt.plot(t[:], s[:], 'g-', label = 'U vs s')

    U_fn_t = np.zeros(n_points_needed)
    s_fn_t = np.zeros(n_points_needed)
    t_fn_t = np.zeros(n_points_needed)

    for j in range(n_points_needed):
        U_fn_t[j] = U[j]
        s_fn_t[j] = s[j]
        t_fn_t[j] = t[j]
    
    z = np.linspace(0,H,n_points_needed)
    U_z = scipy.interpolate.interp1d(s_fn_t, U_fn_t)
    U_fn_z = U_z(z)
    
    t_z = scipy.interpolate.interp1d(s_fn_t, t_fn_t)
    t_fn_z = t_z(z)
    #s_z = scipy.interpolate.interp1d(s_fn_t, U_fn_t)
    
    #plt.plot(z[:], U_fn_z[:], 'go', label = 'U vs s')
    return U_fn_z, s_fn_t, t_fn_z, residence_time, z

def calc_gas_prop(P_CO2_in,T):
    X1 = P_CO2_in/P_tot
    mu_air=calc_viscosity(T)        # named as mu_air because it is currently evaluating for air only and not for flue gas
    mu_gas = mu_air                 # Assumed viscosity of flue gas is same as that of air
    M_gas = 0.028*(1-X1) + 0.044*(X1)
    rho_gas = P_tot*1000/(Rbar/M_gas)/(T)
    nu_gas = mu_gas/rho_gas
    return nu_gas, rho_gas

def U_calc2(Q_l_flux, Q_g_flux, T):
    Q_L = Q_l_flux*A_D/3600        #Liquid flow rate [m3/s] when Q_l_flux is in m3/m2-h
    Q_L_GPM = Q_L*15850.3               #Liquid Flow Rate in Gallons/minute
    d = d_from_Q(Q_L_GPM)*1e-6      #droplet SMD in [m]
    Q_G =Q_g_flux*A_D/3600   # Gas Flow Rate [m3/s] (For aroonwilas paper, the flow rate is assumed total gas flow rate for now as info is not clerarly mentioned in paper)
    U_G = Q_G/A_D            # Gas Velocity [m/s]
    U_jet = Q_L*4/(3.14*d_orifice*d_orifice) # Velocity of the droplet.. jet velocity is used here [m/s]
    
    nu_gas, rho_gas = calc_gas_prop(P_CO2_in,T)
    rho_drop = 1000          # Density of the droplet [kg/m3]. Both MEA and water have a density around 1000 
    U_fn_z, s_fn_z, t_fn_z, t_res, z = U_calc(d, U_jet, nu_gas, U_G, H, rho_gas, rho_drop)
    U = H/t_res         #Average velocity
    return U


def hm_calc(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading, U):    
    
    Q_L = Q_l_flux*A_D/3600        #Liquid flow rate [m3/s] when Q_l_flux is in m3/m2-h
    Q_L_GPM = Q_L*15850.3               #Liquid Flow Rate in Gallons/minute
    
    d = d_from_Q(Q_L_GPM)*1e-6      #droplet SMD in [m]
    R = d/2
    Vol_1d = (4/3)*3.14*(d**3)/8
    
    X1 = P_CO2_in/P_tot          # Mole fraction of CO2 at inlet    
    ####some gas quantities are for air, rest for flue gas, I need to incorporate the air part to flue gas as well later...
    mu_air=calc_viscosity(T)        # named as mu_air because it is currently evaluating for air only and not for flue gas
    mu_gas = mu_air                 # Assumed viscosity of flue gas is same as that of air
    M_gas = 0.028*(1-X1) + 0.044*(X1)
    rho_gas = P_tot*1000/(Rbar/M_gas)/(T)
    nu_gas = mu_gas/rho_gas
    nu_gas, rho_gas = calc_gas_prop(P_CO2_in,T)
    
    ######remove####DCO2_w=1.92e-9      # Diffusivity of CO2 in the droplet [m2/s]. This value is for MEA wt30%
    D_CO2_air = 1.6e-5   # Diffusivity of CO2 in the air [m2/s]
    ######remove#####Henry=4219          # Henry's coefficient in [m3Pa/mol]. This value is for water (AbuArabi 2001)
    rho_drop = 1000       # Density of the droplet [kg/m3]. Both MEA and water have a density around 1000
    mu_drop = 2.48*1e-3        # Dynamic viscosity of the droplet at about 30C and 30wt% or 5M MEA [Pas] Ref: Arachchigeetal 2013
    nu_drop = mu_drop/rho_drop
    
    D_CO2_MEA, H_CO2_MEA = D_and_H_CO2_l_calc(T)
    #print(D_CO2_MEA, H_CO2_MEA)
    
    
        #Cinit=0  # Initial CO2 concentration in the droplet [mol/m3]   .. To be used later when C(z) will be added 
    k2 = k2_calc(T) #m3/mol-s
    C_drop = C_init      #kmol/m3 or M 
    #print(k2)
    print(CO2_loading)
    1/0
    # if CO2_loading > 0.5:
        
    # if CO2_loading < 0.5:
    k_liq_w_kin = ((k2*D_CO2_MEA*C_init*(1-2*CO2_loading)*1000)**0.5/H_CO2_MEA)           #[mol/(Pa-m2-s)] or [mol-s/(kg-m)]
    #print(k_liq_w_kin)
    ###Flow Parameters-    
    
    # Q_G =Q_g_flux*A_D/3600   #Gas Flow Rate [m3/s] (For aroonwilas paper, the flow rate is assumed total gas flow rate for now as info is not clerarly mentioned in paper)
    # U_G = Q_G/A_D            #Gas Velocity [m/s]
    U_jet = Q_L*4/(3.14*d_orifice*d_orifice) # Velocity of the droplet.. jet velocity is used here [m/s]
    # # U_fn_z, s_fn_z, t_fn_z, t_res, z = U_calc(d, U_jet, nu_gas, U_G, H, rho_gas, rho_drop) ### Gives velocity of droplet at each location of height, and droplet residence time
    # U = H/t_res         #Average velocity
    # #U = U_jet
    
    #print('Gas velocity is', U_G, 'and droplet velocity is', U)
    # tau = d*d/(18*nu_drop)      #Time constant for velocity variation
    # U_rel_terminal = g*tau
    # U_terminal = U_rel_terminal-U_G
    # #print( 'Gas velocity is', U_G, 'and droplet velocity is', U, 'and terminal velocity is', U_terminal)
    # U_array[j] = U_terminal
    # print(U_terminal, U, U_G)


    
    Re_drop = U*d/nu_gas            ###Should it be relative velocity or absolute fluid velocity
    Sc_drop = nu_gas/D_CO2_air        # Check these!! Also check if the diffusion represented here is co2 or amine    
    Sh_drop = 2 + 0.552*(Re_drop**(1/2))*(Sc_drop**(1/3))   #Froessling Equation
    #print(Sh_drop, 'is sherwood number')
    #print('nu_gas=',nu_gas)
    Re_jet = U_jet*d/nu_gas
    Stop_dist = (rho_drop*d/rho_gas)*(Re_jet**(1/3)-math.sqrt(6)*math.atan(Re_jet**(1/3)/math.sqrt(6)))
    #print('S=', Stop_dist)
    
    kg=Sh_drop*D_CO2_air/(Rbar*T*d)               # [mol/(Pa-m2-s)] or [mol-s/(kg-m)]
    #print(kg)
    #kG = kg
    kG = 1/(1/kg + 1/k_liq_w_kin)
    hm = kG*Rbar*T                                  #[m/s]
    # print('kL = ', k_liq_w_kin, 'kg = ', kg, 'and kG = ', kG)
    
    #print('Kgav in mol/Pa-m3-hr assuming area for all droplets = ', kG*3600*H*Q_L/(U*Vol_1d*A_D*H) )
    #print(hm, 'is hm')
    
    Bi=kg*(d)*Rbar*T/D_CO2_air
    #print("Biot is %2.2f" %(Bi))
    return hm, d
```


```

```