Skip to content

Commit

Permalink
Improve tools absorption_function.py
Browse files Browse the repository at this point in the history 
- Can now fit analytic functions or numerical absorption function
- Clean up code
  • Loading branch information
@jannisteunissen
jannisteunissen committed Nov 22, 2023
1 parent b94432f commit de9c9d4
Showing 1 changed file with 101 additions and 52 deletions.
153 changes: 101 additions & 52 deletions tools/absorption_function.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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')
Expand Down Expand Up @@ -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

Expand All @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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))
Expand All @@ -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]
Expand All @@ -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))

Expand All @@ -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)
Expand All @@ -221,47 +289,32 @@ 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)

plt.subplot(121)
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()
Expand All @@ -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()
Expand Down

0 comments on commit de9c9d4

Please sign in to comment.