From de9c9d4bf91af6f03276e3b0981e81045fb02e6f Mon Sep 17 00:00:00 2001 From: Jannis Teunissen Date: Wed, 22 Nov 2023 10:35:53 +0100 Subject: [PATCH] Improve tools absorption_function.py - Can now fit analytic functions or numerical absorption function - Clean up code --- tools/absorption_function.py | 153 +++++++++++++++++++++++------------ 1 file changed, 101 insertions(+), 52 deletions(-) diff --git a/tools/absorption_function.py b/tools/absorption_function.py index 72c773d3..fcad794b 100755 --- a/tools/absorption_function.py +++ b/tools/absorption_function.py @@ -5,7 +5,7 @@ # Behnaz Bagheri # from sympy import * -from scipy.integrate import quad, simps +from scipy.integrate import quad, simpson from scipy.optimize import curve_fit from scipy.interpolate import interp1d import numpy as np @@ -26,12 +26,17 @@ def get_args(): help='Partial pressures of gases (bar)') p.add_argument('-fit_range', nargs=2, type=float, default=[1e-4, 3e-3], help='Distance range for fit of coefficients') - p.add_argument('-n_modes', type=int, default=2, + p.add_argument('-n_modes', type=int, default=3, help='Number of Helmholtz modes') + p.add_argument('-H2O_model', type=str, choices=['Naidis', 'Aints'], + default='Naidis', help='Type of H2O absorption model') p.add_argument('-guess_amplitudes', type=float, help='Initial guess for amplitude of Helmholtz modes') p.add_argument('-guess_lambdas', type=float, help='Initial guess for lambdas of Helmholtz modes') + p.add_argument('-fit_what', type=str, default='numerical', + choices=['numerical', 'Naidis', 'Aints'], + help='What type of data/function to fit') p.add_argument('-fit_type', type=str, default='least_squares', choices=['least_squares', 'relative', 'log'], help='What type of errors to use in fit') @@ -72,10 +77,37 @@ def get_args(): # convert to 1/(m * bar) Aints Plasma Processes and Polymers 2008 5, 672-680 K_H2O_MAX = 0.57e2 / ONE_TORR + +def zheleznyak_f(r): + """Analytic absorption function at distance r (m) for dry air, unit: 1/m + From Zheleznyak et al. 1982 + """ + return (np.exp(-MU_MIN * p['O2'] * r) - + np.exp(-MU_MAX * p['O2'] * r)) / (r * np.log(MU_MAX / MU_MIN)) + + +def naidis_moist_zheleznyak_f(r): + """Analytic absorption function at distance r (m) for moist air, unit: 1/m + This approximation is based on Naidis 2006 DOI:10.1088/0963-0252/15/2/010 + """ + nom = (np.exp(-(MU_MIN * p['O2'] + K_H2O * p['H2O']) * r) - + np.exp(-(MU_MAX * p['O2'] + K_H2O * p['H2O']) * r)) + denom = (r * np.log(MU_MAX / MU_MIN)) + return nom/denom + + +def aints_moist_zheleznyak_f(r): + """Analytic absorption function at distance r (m) for most air, unit: 1/m + This approximation is based on Aints 2008 DOI:10.1002/ppap.200800031 + """ + nom = (np.exp(-(MU_MIN * p['O2'] + K_H2O_MIN * p['H2O']) * r) - + np.exp(-(MU_MAX * p['O2'] + K_H2O_MAX * p['H2O']) * r)) + denom = (r * np.log((MU_MAX * p['O2'] + K_H2O_MAX * p['H2O']) / + (MU_MIN * p['O2'] + K_H2O_MIN * p['H2O']))) + return nom/denom + + class Absorption_O2: - """Zheleznyak approximation for O2 absorption coefficient Reference: - "Photoionization of nitrogen and oxygen mixtures by radiation from a gas - discharge" by Zheleznyak and Mnatsakanian, 1982""" gas = "O2" ionizing = True @@ -85,6 +117,8 @@ def set_pressure(self, p): def mu(self, x): """Absorption function at wavelength x (nm). Unit: 1 / (m * bar)""" + + # From Zheleznyak et al. 1982 t = (1/x - 1/LAMBDA_MAX)/(1/LAMBDA_MIN - 1/LAMBDA_MAX) return self.p * MU_MIN * (MU_MAX/MU_MIN)**t @@ -108,17 +142,36 @@ class Absorption_H2O: gas = "H2O" ionizing = False + def __init__(self, model='Naidis'): + self.model = model + print(self.model) + def set_pressure(self, p): """Set partial pressure in bar""" self.p = p def mu(self, x): """Absorption function at wavelength x (nm). Unit: 1 / (m * bar)""" - return self.p * 0.3e2 / ONE_TORR + + if self.model == 'Naidis': + # This approximation is based on Naidis 2006 + # DOI:10.1088/0963-0252/15/2/010 + mu_x = self.p * 0.26e2 / ONE_TORR + elif self.model == 'Aints': + # This approximation is based on Aints 2008 + # DOI:10.1002/ppap.200800031 + t = (1/x - 1/LAMBDA_MAX)/(1/LAMBDA_MIN - 1/LAMBDA_MAX) + mu_x = self.p * K_H2O_MIN * (K_H2O_MAX/K_H2O_MIN)**t + else: + raise ValueError('Unknown H2O absorption model') + + return mu_x # List of all gases -all_gases = [Absorption_O2(), Absorption_CO2(), Absorption_H2O()] +all_gases = [Absorption_O2(), + Absorption_CO2(), + Absorption_H2O(args.H2O_model)] all_names = [x.gas for x in all_gases] # Construct lists of absorbing and ionizing gases present @@ -154,18 +207,29 @@ def integrand(x, r): r = np.linspace(args.fit_range[0], args.fit_range[1], args.n_points) -f = np.zeros(args.n_points) +f_numerical = np.zeros(args.n_points) # Construct absorption function numerically integrated over wavelength for i in range(args.n_points): - f[i] = quad(integrand, LAMBDA_MIN, LAMBDA_MAX, - args=(r[i]))[0] / (LAMBDA_MAX - LAMBDA_MIN) + tmp = quad(integrand, LAMBDA_MIN, LAMBDA_MAX, args=(r[i])) + f_numerical[i] = tmp[0] / (LAMBDA_MAX - LAMBDA_MIN) -lambda_guess = -np.log(f[-2]/f[-1])/(r[-2] - r[-1]) + +if args.fit_what == 'numerical': + f_to_fit = f_numerical +elif args.fit_what == 'Naidis': + f_to_fit = naidis_moist_zheleznyak_f(r) +elif args.fit_what == 'Aints': + f_to_fit = aints_moist_zheleznyak_f(r) +else: + raise ValueError('Invalid argument for fit_what') + + +lambda_guess = -np.log(f_to_fit[-2]/f_to_fit[-1])/(r[-2] - r[-1]) amplitude_guess = lambda_guess**2 / args.n_modes # Construct interpolating function -f_func = interp1d(r, f) +f_func = interp1d(r, f_to_fit) # Guesses for the fit fit_guess = np.ones((2 * args.n_modes)) @@ -179,6 +243,7 @@ def integrand(x, r): else: fit_guess[1::2] = lambda_guess + # Function for fitting with a variable number of exponential functions def fit_func(x, *fitargs): amplitudes = fitargs[0::2] @@ -188,9 +253,11 @@ def fit_func(x, *fitargs): val += x * c1 * np.exp(-c2 * x) return val + def fit_func_relative(x, *fitargs): return fit_func(x, *fitargs)/f_func(x) + def fit_func_log(x, *fitargs): return np.log(fit_func(x, *fitargs)) @@ -200,9 +267,10 @@ def fit_func_log(x, *fitargs): popt, pcov = curve_fit(fit_func_relative, r, np.ones(args.n_points), p0=fit_guess) if args.fit_type == 'log': - popt, pcov = curve_fit(fit_func_log, r, np.log(f), p0=fit_guess) + popt, pcov = curve_fit(fit_func_log, r, + np.log(f_to_fit), p0=fit_guess) else: - popt, pcov = curve_fit(fit_func, r, f, p0=fit_guess) + popt, pcov = curve_fit(fit_func, r, f_to_fit, p0=fit_guess) except RuntimeError as e: print('No convergence, adjust guess_amplitudes and/or guess_lambdas') print(e) @@ -221,30 +289,19 @@ def fit_func_log(x, *fitargs): np.mean(np.sqrt(np.diag(pcov)) / np.abs(popt)))) print('Fit range (m): {:.5e} -- {:.5e}'.format( args.fit_range[0], args.fit_range[1])) -print('Integral of absorption function over fit range: {:.5e}'.format( - simps(f, r))) - - -def zheleznyak_f(r): - """Analytic absorption function at distance r (m) for dry air, unit: 1/m""" - return (np.exp(-MU_MIN * p['O2'] * r) - - np.exp(-MU_MAX * p['O2'] * r)) / \ - (r * np.log(MU_MAX / MU_MIN)) +f_to_display = [] +if args.show_Zheleznyak: + f_to_display.append(['Zheleznyak air', zheleznyak_f]) +if args.show_Naidis_moist: + f_to_display.append(['Naidis moist air', naidis_moist_zheleznyak_f]) +if args.show_Aints_moist: + f_to_display.append(['Aints moist air', aints_moist_zheleznyak_f]) -def naidis_moist_zheleznyak_f(r): - """Analytic absorption function at distance r (m) for moist air, unit: 1/m""" - return (np.exp(-(MU_MIN * p['O2'] + K_H2O * p['H2O']) * r) - - np.exp(-(MU_MAX * p['O2'] + K_H2O * p['H2O']) * r)) / \ - (r * np.log(MU_MAX / MU_MIN)) - -def aints_moist_zheleznyak_f(r): - """Analytic absorption function at distance r (m) for most air, unit: 1/m""" - return (np.exp(-(MU_MIN * p['O2'] + K_H2O_MIN * p['H2O']) * r) - - np.exp(-(MU_MAX * p['O2'] + K_H2O_MAX * p['H2O']) * r)) / \ - (r * np.log( (MU_MAX * p['O2'] + K_H2O_MAX * p['H2O']) / - (MU_MIN * p['O2'] + K_H2O_MIN * p['H2O']))) - +print('Integrals of absorption functions over fit range:') +print(f'{"Numerical":<20} {simpson(f_numerical, r):12.5f}') +for name, f in f_to_display: + print(f'{name:<20} {simpson(f(r), r):12.5f}') fig, ax = plt.subplots(1, 2) @@ -252,16 +309,12 @@ def aints_moist_zheleznyak_f(r): plt.xlabel('r (m)') plt.ylabel('absorption function (1/m)') plt.title('Logarithmic scale') -plt.semilogy(r, f, '.', label='numerical') +plt.semilogy(r, f_numerical, '.', label='numerical') +plt.semilogy(r, fit_func(r, *popt), label='fit ({}-term)'.format(args.n_modes)) -if args.show_Zheleznyak: - plt.semilogy(r, zheleznyak_f(r), '--', label='Zheleznyak air') -if args.show_Naidis_moist: - plt.semilogy(r, naidis_moist_zheleznyak_f(r), '--', label='Naidis moist air') -if args.show_Aints_moist: - plt.semilogy(r, aints_moist_zheleznyak_f(r), '--', label='Aints moist air') +for name, f in f_to_display: + plt.semilogy(r, f(r), '--', label=name) -plt.semilogy(r, fit_func(r, *popt), label='fit ({}-term)'.format(args.n_modes)) if args.show_curve: plt.semilogy(r, fit_func(r, *args.show_curve), label='custom') plt.legend() @@ -270,16 +323,12 @@ def aints_moist_zheleznyak_f(r): plt.xlabel('r (m)') plt.ylabel('absorption function (1/m)') plt.title('Linear scale') -plt.plot(r, f, '.', label='numerical') +plt.plot(r, f_numerical, '.', label='numerical') +plt.plot(r, fit_func(r, *popt), label='fit ({}-term)'.format(args.n_modes)) -if args.show_Zheleznyak: - plt.plot(r, zheleznyak_f(r), '--', label='Zheleznyak air') -if args.show_Naidis_moist: - plt.plot(r, naidis_moist_zheleznyak_f(r), '--', label='Naidis moist air') -if args.show_Aints_moist: - plt.plot(r, aints_moist_zheleznyak_f(r), '--', label='Aints moist air') +for name, f in f_to_display: + plt.plot(r, f(r), '--', label=name) -plt.plot(r, fit_func(r, *popt), label='fit ({}-term)'.format(args.n_modes)) if args.show_curve: plt.plot(r, fit_func(r, *args.show_curve), label='custom') plt.legend()