From 49acc8810a5bf7c29cbcdc251e94c29bda217739 Mon Sep 17 00:00:00 2001 From: Clement Moissard Date: Wed, 13 Jul 2022 21:23:49 +0100 Subject: [PATCH 1/4] Added fit for ionisation from Fermi model --- hack/Thomas_Fermi_ionisation.py | 41 +++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 hack/Thomas_Fermi_ionisation.py diff --git a/hack/Thomas_Fermi_ionisation.py b/hack/Thomas_Fermi_ionisation.py new file mode 100644 index 0000000..d38fb28 --- /dev/null +++ b/hack/Thomas_Fermi_ionisation.py @@ -0,0 +1,41 @@ +import numpy as np +import astropy.units as u +from plasmapy.particles import * + + +def Zav_TF(particle: Particle, rho: u.g/u.cm , T: u.eV=0) -> u.dimensionless_unscaled: + """ + Finite Temperature Thomas Fermi Charge State using + R.M. More (1985), "Pressure Ionization, Resonances, and the + Continuity of Bound and Free States", Adv. in Atomic + Mol. Phys., Vol. 21, p. 332 (Table IV). + + This is a light-weight fit of the Thomas-Fermi model, not the model itself. + """ + + Z = atomic_number(particle) + A = mass_number(particle) + + alpha = 14.3139 + beta = 0.6624 + a1 = 0.003323 + a2 = 0.9718 + a3 = 9.26148e-5 + a4 = 3.10165 + b0 = -1.7630 + b1 = 1.43175 + b2 = 0.31546 + c1 = -0.366667 + c2 = 0.983333 + + rho1 = rho / A * Z + T1 = T / Z ** (4. / 3.) + Tf = T1 / (1 + T1) + Ac = a1 * T1 ** a2 + a3 * T1 ** a4 + B = -np.exp(b0 + b1 * Tf + b2 * Tf ** 7) + C = c1 * Tf + c2 + Q1 = Ac * rho1 ** B + Q = (rho1 ** C + Q1 ** C) ** (1 / C) + x = alpha * Q ** beta + + return Z * x / (1 + x + np.sqrt(1 + 2. * x)) \ No newline at end of file From deb325a84a3789efded3f9f2b2798dacf32402bd Mon Sep 17 00:00:00 2001 From: Clement Moissard Date: Wed, 13 Jul 2022 21:38:02 +0100 Subject: [PATCH 2/4] Added fit for ionisation from Fermi model --- hack/{ => Ionisation}/Thomas_Fermi_ionisation.py | 0 hack/Ionisation/__init__.py | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename hack/{ => Ionisation}/Thomas_Fermi_ionisation.py (100%) create mode 100644 hack/Ionisation/__init__.py diff --git a/hack/Thomas_Fermi_ionisation.py b/hack/Ionisation/Thomas_Fermi_ionisation.py similarity index 100% rename from hack/Thomas_Fermi_ionisation.py rename to hack/Ionisation/Thomas_Fermi_ionisation.py diff --git a/hack/Ionisation/__init__.py b/hack/Ionisation/__init__.py new file mode 100644 index 0000000..e69de29 From da0321a6c93213a4e8c42051ab24c3a0d53bb66c Mon Sep 17 00:00:00 2001 From: Clement Moissard Date: Thu, 14 Jul 2022 20:10:46 +0100 Subject: [PATCH 3/4] Corrected fit for ionisation from Fermi model - now this works --- hack/Ionisation/Thomas_Fermi_ionisation.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/hack/Ionisation/Thomas_Fermi_ionisation.py b/hack/Ionisation/Thomas_Fermi_ionisation.py index d38fb28..7e58ff0 100644 --- a/hack/Ionisation/Thomas_Fermi_ionisation.py +++ b/hack/Ionisation/Thomas_Fermi_ionisation.py @@ -1,9 +1,11 @@ import numpy as np import astropy.units as u -from plasmapy.particles import * +from plasmapy.particles import Particle, particle_input +from plasmapy.utils.decorators import validate_quantities - -def Zav_TF(particle: Particle, rho: u.g/u.cm , T: u.eV=0) -> u.dimensionless_unscaled: +@validate_quantities(density={"can_be_negative": False}, equivalencies=u.temperature_energy()) +@particle_input +def Zav_TF(particle: Particle, rho: u.g/u.cm**-3 , T: u.eV=0) -> u.dimensionless_unscaled: """ Finite Temperature Thomas Fermi Charge State using R.M. More (1985), "Pressure Ionization, Resonances, and the @@ -13,8 +15,8 @@ def Zav_TF(particle: Particle, rho: u.g/u.cm , T: u.eV=0) -> u.dimensionless_uns This is a light-weight fit of the Thomas-Fermi model, not the model itself. """ - Z = atomic_number(particle) - A = mass_number(particle) + Z = particle.atomic_number + A = particle.mass_number alpha = 14.3139 beta = 0.6624 @@ -28,8 +30,8 @@ def Zav_TF(particle: Particle, rho: u.g/u.cm , T: u.eV=0) -> u.dimensionless_uns c1 = -0.366667 c2 = 0.983333 - rho1 = rho / A * Z - T1 = T / Z ** (4. / 3.) + rho1 = rho.value / A * Z + T1 = T.value / Z ** (4. / 3.) Tf = T1 / (1 + T1) Ac = a1 * T1 ** a2 + a3 * T1 ** a4 B = -np.exp(b0 + b1 * Tf + b2 * Tf ** 7) From 6eb94cb65259368e110d482bf8d92dcb8af70d44 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9ment=20Moissard?= Date: Fri, 15 Jul 2022 14:27:40 +0100 Subject: [PATCH 4/4] Update hack/Ionisation/Thomas_Fermi_ionisation.py Co-authored-by: Erik Everson --- hack/Ionisation/Thomas_Fermi_ionisation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hack/Ionisation/Thomas_Fermi_ionisation.py b/hack/Ionisation/Thomas_Fermi_ionisation.py index 7e58ff0..9752bfe 100644 --- a/hack/Ionisation/Thomas_Fermi_ionisation.py +++ b/hack/Ionisation/Thomas_Fermi_ionisation.py @@ -40,4 +40,4 @@ def Zav_TF(particle: Particle, rho: u.g/u.cm**-3 , T: u.eV=0) -> u.dimensionless Q = (rho1 ** C + Q1 ** C) ** (1 / C) x = alpha * Q ** beta - return Z * x / (1 + x + np.sqrt(1 + 2. * x)) \ No newline at end of file + return Z * x / (1 + x + np.sqrt(1 + 2. * x))