From c4c327aedcd63e15e43da0b6ad7eccda2db83bc4 Mon Sep 17 00:00:00 2001
From: Jannis Teunissen <jannis@teunissen.net>
Date: Wed, 22 Nov 2023 11:51:50 +0100
Subject: [PATCH] absorption_function.py: show quenching pressure for Aints
 model

---
 tools/absorption_function.py | 18 +++++++++++++++---
 1 file changed, 15 insertions(+), 3 deletions(-)

diff --git a/tools/absorption_function.py b/tools/absorption_function.py
index fcad794b..a2269be2 100755
--- a/tools/absorption_function.py
+++ b/tools/absorption_function.py
@@ -35,11 +35,13 @@ def get_args():
     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'],
+                   choices=['numerical', 'Zheleznyak', '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')
+    p.add_argument('-ptot_for_quenching', type=float,
+                   help='Total gas pressure (bar) to show quenching info')
     p.add_argument('-show_Zheleznyak', action='store_true',
                    help='Show Zheleznyak curve for air')
     p.add_argument('-show_Naidis_moist', action='store_true',
@@ -69,6 +71,7 @@ def get_args():
 # Zheleznyak coefficients for air
 MU_MAX = 2.0e2 / ONE_TORR       # convert to 1/(m * bar)
 MU_MIN = 0.035e2 / ONE_TORR     # convert to 1/(m * bar)
+pq_air_Zheleznyak = 30 * ONE_TORR  # Quenching parameter
 
 # convert to 1/(m * bar) Naidis 2006 PSST 15 253
 K_H2O = 0.26e2 / ONE_TORR
@@ -76,6 +79,7 @@ def get_args():
 K_H2O_MIN = 0.13e2 / ONE_TORR
 # convert to 1/(m * bar) Aints Plasma Processes and Polymers 2008 5, 672-680
 K_H2O_MAX = 0.57e2 / ONE_TORR
+pq_H2O_Aints = 0.5 * ONE_TORR   # Quenching parameter
 
 
 def zheleznyak_f(r):
@@ -177,7 +181,8 @@ def mu(self, x):
 # Construct lists of absorbing and ionizing gases present
 f_absorption = []
 f_ionizing = []
-p = {}                          # Pressures
+# Pressures
+p = {'O2': 0., 'CO2': 0., 'H2O': 0.}
 
 for gas, pressure in zip(args.gases, args.pressures):
     ix = all_names.index(gas)
@@ -217,10 +222,17 @@ def integrand(x, r):
 
 if args.fit_what == 'numerical':
     f_to_fit = f_numerical
-elif args.fit_what == 'Naidis':
+elif args.fit_what == 'Zheleznyak':
     f_to_fit = naidis_moist_zheleznyak_f(r)
 elif args.fit_what == 'Aints':
     f_to_fit = aints_moist_zheleznyak_f(r)
+    # Show quenching information
+    if args.ptot_for_quenching is not None:
+        Q = (1 + (args.ptot_for_quenching - p['H2O'])/pq_air_Zheleznyak +
+             p['H2O']/pq_H2O_Aints)**-1
+        pq_effective = -Q * args.ptot_for_quenching/(Q - 1)
+        print('Effective quenching pressure for Aints model: '
+              f'{pq_effective:.5e} bar')
 else:
     raise ValueError('Invalid argument for fit_what')