diff --git a/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/BoundBoundTransitions.py b/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/BoundBoundTransitions.py new file mode 100755 index 0000000000..f1766733da --- /dev/null +++ b/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/BoundBoundTransitions.py @@ -0,0 +1,211 @@ +#!/usr/bin/env python +""" +atomicPhysics(FLYonPIC) reference rate calculation +This file is part of the PIConGPU. +Copyright 2023-2024 PIConGPU contributors +Authors: Brian Marre, Axel Huebl +License: GPLv3+ +""" + +import numpy as np +import scipy.constants as const +import scipy.special as scipy + +""" @file reference implementation of the rate calculation for bound-bound transitions """ + + +class BoundBoundTransitions: + @staticmethod + def _gaunt(U, cxin1, cxin2, cxin3, cxin4, cxin5): + """calculate gaunt factor + + @param U float (energy interacting electron)/(delta Energy Transition), unitless + @param cxin1 float gaunt approximation coefficient 1, unitless + @param cxin2 float gaunt approximation coefficient 2, unitless + @param cxin3 float gaunt approximation coefficient 3, unitless + @param cxin4 float gaunt approximation coefficient 4, unitless + @param cxin5 float gaunt approximation coefficient 5, unitless + + @return unitless + """ + + if U < 1.0: + return 0.0 + else: + return cxin1 * np.log(U) + cxin2 + cxin3 / (U + cxin5) + cxin4 / (U + cxin5) ** 2 + + @staticmethod + def _multiplicity(levelVector): + """get degeneracy of atomicState with given levelVector""" + result = 1.0 + for i, n_i in enumerate(levelVector): + result *= scipy.comb(2 * (i + 1) ** 2, n_i) + return result + + @staticmethod + def collisionalBoundBoundCrossSection( + energyElectron, + energyDiffLowerUpper, + collisionalOscillatorStrength, + cxin1, + cxin2, + cxin3, + cxin4, + cxin5, + lowerStateLevelVector, + upperStateLevelVector, + excitation, + ): + """cross section for de-/excitation transition due to interaction with free electron + + @param excitation bool true =^= excitation, false =^= deexcitation + @param energyElectron float energy of interacting free electron, [eV] + + @param energyDiffLowerUpper upperStateEnergy - lowerStateEnergy, [eV] + @param collisionalOscillatorStrength of the transition, unitless + @param cxin1 float gaunt approximation coefficient 1, unitless + @param cxin2 float gaunt approximation coefficient 2, unitless + @param cxin3 float gaunt approximation coefficient 3, unitless + @param cxin4 float gaunt approximation coefficient 4, unitless + @param cxin5 float gaunt approximation coefficient 5, unitless + + @return unit 1e6b, 10^-22 m^2 + """ + if energyDiffLowerUpper > energyElectron: + return 0.0 + + U = energyElectron / energyDiffLowerUpper # unitless + E_Rydberg = const.physical_constants["Rydberg constant times hc in eV"][0] # eV + + a0 = const.physical_constants["Bohr radius"][0] # m + c0 = 8 * (np.pi * a0) ** 2 / np.sqrt(3.0) # m^2 + + crossSection_butGaunt = ( + c0 + * (E_Rydberg / energyDiffLowerUpper) ** 2 + * collisionalOscillatorStrength + * (energyDiffLowerUpper / energyElectron) + / 1e-22 + ) + # 1e6b + #! @attention differs from original publication + + if excitation: + return crossSection_butGaunt * BoundBoundTransitions._gaunt(U, cxin1, cxin2, cxin3, cxin4, cxin5) # 1e6b + else: + statisticalRatio = BoundBoundTransitions._multiplicity( + lowerStateLevelVector + ) / BoundBoundTransitions._multiplicity(upperStateLevelVector) # unitless + return ( + statisticalRatio + * crossSection_butGaunt + * BoundBoundTransitions._gaunt(U + 1.0, cxin1, cxin2, cxin3, cxin4, cxin5) + ) # 1e6b + + @staticmethod + def rateCollisionalBoundBoundTransition( + energyElectron, + energyElectronBinWidth, + densityElectrons, + energyDiffLowerUpper, + collisionalOscillatorStrength, + cxin1, + cxin2, + cxin3, + cxin4, + cxin5, + lowerStateLevelVector, + upperStateLevelVector, + excitation=True, + ): + """rate of collisional de-/excitation + + @param excitation bool true =^= excitation, false =^= deexcitation + + @param energyElectron float central energy of electron bin, [eV] + @param energyElectronBinWidth width of energy bin, [eV] + @param densityElectrons number density of physical electrons in bin, 1/(m^3 * eV) + + @param energyDiffLowerUpper upperStateEnergy - lowerStateEnergy, [eV] + @param collisionalOscillatorStrength of the transition, unitless + @param cxin1 float gaunt approximation coefficient 1 + @param cxin2 float gaunt approximation coefficient 2 + @param cxin3 float gaunt approximation coefficient 3 + @param cxin4 float gaunt approximation coefficient 4 + @param cxin5 float gaunt approximation coefficient 5 + + @return unit 1/s + """ + + sigma = BoundBoundTransitions.collisionalBoundBoundCrossSection( + energyElectron, + energyDiffLowerUpper, + collisionalOscillatorStrength, + cxin1, + cxin2, + cxin3, + cxin4, + cxin5, + lowerStateLevelVector, + upperStateLevelVector, + excitation, + ) # 1e6b + + if sigma < 0.0: + sigma = 0.0 + + electronRestMassEnergy = const.physical_constants["electron mass energy equivalent in MeV"][0] * 1e6 # eV + + # derivation for relativist kinetic energy to velocity + # E = gamma * m*c^2 E_kin + m*c^2 = gamma * m*c^2 => E_kin = (gamma-1) * m*c^2 + # => gamma = E_kin / (m*c^2) + 1 + # 1/sqrt(1-beta^2) = E_kin / (m*c^2) + 1 + # 1/(1-beta^2) = (E_kin/(m*c^2) + 1)^2 + # 1-beta^2 = 1/(E_kin/(m*c^2) + 1)^2 + # beta^2 = 1 - 1/(E_kin/(m*c^2) + 1)^2 + # => beta = sqrt(1 - 1/(1 + E_kin/(m*c^2))^2) + # v = c * beta = c * sqrt(1 - 1/(1 + E_kin/(m*c^2))^2) + + # dE * sigma(E) * rho_e * v + # eV * 1e6b * m^2/(1e6b) * 1/(m^3 * eV) * m/s * unitless + # = eV/(eV) * m^2 * 1/m^3 * m/s = 1/s + return ( + energyElectronBinWidth + * sigma + * 1e-22 + * densityElectrons + * const.physical_constants["speed of light in vacuum"][0] + * np.sqrt(1.0 - 1.0 / (1.0 + energyElectron / electronRestMassEnergy) ** 2) + ) # 1/s + + @staticmethod + def rateSpontaneousDeexcitation( + absorptionOscillatorStrength, frequencyPhoton, lowerStateLevelVector, upperStateLevelVector + ): + """rate of spontaneous deexcitation under photon emission + + @param absorptionOscillatorStrength unitless + @param frequencyPhoton [1/s] + @param lowerStateLevelVector occupation number level vector for lower state of transition + i.e. how many electron occupy states in each principal quantum number shell + @param upperStateLevelVector same as lowerStateLevelVector but upper state of transition + + @return unit 1/s + """ + # (As)^2 * N/A^2 * / (kg * m/s ) = A^2 * s^2 kg*m/s^2 * 1/A^2 /(kg *m/s ) + # = A^2/A^2 * s^2/s^2 * (kg*m)/(kg*m) * 1/(1/s) = s + scalingConstant = ( + 2 + * np.pi + * const.physical_constants["elementary charge"][0] ** 2 + * const.physical_constants["vacuum mag. permeability"][0] + / (const.value("electron mass") * const.physical_constants["speed of light in vacuum"][0]) + ) # s + + ratio = BoundBoundTransitions._multiplicity(lowerStateLevelVector) / BoundBoundTransitions._multiplicity( + upperStateLevelVector + ) + # untiless + + # s * 1/s^2 * untiless * unitless = 1/s + return scalingConstant * frequencyPhoton**2 * ratio * absorptionOscillatorStrength # 1/s diff --git a/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/BoundFreeCollisionalTransitions.py b/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/BoundFreeCollisionalTransitions.py new file mode 100755 index 0000000000..d3917c7c80 --- /dev/null +++ b/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/BoundFreeCollisionalTransitions.py @@ -0,0 +1,139 @@ +#!/usr/bin/env python +""" +atomicPhysics(FLYonPIC) reference rate calculation +This file is part of the PIConGPU. +Copyright 2023-2024 PIConGPU contributors +Authors: Brian Marre, Axel Huebl +License: GPLv3+ +""" + +import numpy as np +import scipy.constants as const +import scipy.special as scipy + +""" @file reference implementation of the rate calculation for bound-free collisional transitions """ + + +class BoundFreeCollisionalTransitions: + @staticmethod + def _betaFactor(screenedCharge): + return 0.25 * (np.sqrt((100.0 * screenedCharge + 91.0) / (4.0 * screenedCharge + 3.0)) - 5) # unitless + + @staticmethod + def _wFactor(U, screenedCharge): + return np.power(np.log(U), BoundFreeCollisionalTransitions._betaFactor(screenedCharge) / U) # unitless + + @staticmethod + def _multiplicity(lowerStateLevelVector, upperStateLevelVector): + result = 1.0 + for i in range(len(lowerStateLevelVector)): + result *= scipy.comb(lowerStateLevelVector[i], lowerStateLevelVector[i] - upperStateLevelVector[i]) + return result + + @staticmethod + def collisionalIonizationCrossSection( + energyElectron, + ionizationEnergy, + excitationEnergyDifference, + screenedCharge, + lowerStateLevelVector, + upperStateLevelVector, + ): + """get cross section for ionization by interaction with a free electron + + @param energyElectron float energy of interacting free electron, [eV] + @param ionizationEnergy float sum of ionization energies between initial and final charge state, [eV] + @param excitationEnergyDifference float excitation energy upper state - excitation energy lower state + all relative to each respective's charge state ground state + @param screenedCharge screenedCharge of the ionizing shell + @param lowerStateLevelVector occupation number level vector for lower state of transition + i.e. how many electron occupy states in each principal quantum number shell + @param upperStateLevelVector same as lowerStateLevelVector but upper state of transition + + @return unit: 1e6b + """ + + energyDifference = ionizationEnergy + excitationEnergyDifference + U = energyElectron / energyDifference + combinatorialFactor = BoundFreeCollisionalTransitions._multiplicity( + lowerStateLevelVector, upperStateLevelVector + ) + + # m^2 * (eV/(eV))^2 * 1/(eV/eV) * unitless * unitless / (m^2/1e6b) = 1e6b + if U > 1: + return ( + np.pi + * const.value("Bohr radius") ** 2 + * 2.3 + * combinatorialFactor + * (const.value("Rydberg constant times hc in eV") / energyDifference) ** 2 + * 1.0 + / U + * np.log(U) + * BoundFreeCollisionalTransitions._wFactor(U, screenedCharge) + ) / 1e-22 # 1e6b, 1e-22 m^2 + else: + return 0.0 + + @staticmethod + def rateCollisionalIonization( + energyElectron, + energyElectronBinWidth, + densityElectrons, + ionizationEnergy, + excitationEnergyDifference, + screenedCharge, + lowerStateLevelVector, + upperStateLevelVector, + ): + """rate of ionization due to interaction with free electron + + @param energyElectron float central energy of electron bin, [eV] + @param energyElectronBinWidth float width of energy bin, [eV] + @param densityElectrons float number density of physical electrons in bin, 1/(m^3 * eV) + + @param ionizationEnergy sum of ionization energies between initial and final charge state, [eV] + @param excitationEnergyDifference float excitation energy upper state - excitation energy lower state + all relative to each respective's charge state ground state + @param screenedCharge screenedCharge of the ionizing shell + @param lowerStateLevelVector occupation number level vector for lower state of transition + i.e. how many electron occupy states in each principal quantum number shell + @param upperStateLevelVector same as lowerStateLevelVector but upper state of transition + + @return unit: 1/s + """ + sigma = BoundFreeCollisionalTransitions.collisionalIonizationCrossSection( + energyElectron, + ionizationEnergy, + excitationEnergyDifference, + screenedCharge, + lowerStateLevelVector, + upperStateLevelVector, + ) # 1e6b + + if sigma < 0.0: + sigma = 0.0 + + electronRestMassEnergy = const.value("electron mass energy equivalent in MeV") * 1e6 # eV + + # derivation for relativist kinetic energy to velocity + # E = gamma * m*c^2 E_kin + m*c^2 = gamma * m*c^2 => E_kin = (gamma-1) * m*c^2 + # => gamma = E_kin / (m*c^2) + 1 + # 1/sqrt(1-beta^2) = E_kin / (m*c^2) + 1 + # 1/(1-beta^2) = (E_kin/(m*c^2) + 1)^2 + # 1-beta^2 = 1/(E_kin/(m*c^2) + 1)^2 + # beta^2 = 1 - 1/(E_kin/(m*c^2) + 1)^2 + # => beta = sqrt(1 - 1/(1 + E_kin/(m*c^2))^2) + # v = c * beta = c * sqrt(1 - 1/(1 + E_kin/(m*c^2))^2) + + # dE * sigma(E) * rho_e * v + # eV * 1e6b * m^2/(1e6b) * 1/(m^3 * eV) * m/s * unitless + # = eV/(eV) * m^2 * 1/m^3 * m/s = 1/s + return ( + energyElectronBinWidth + * sigma + * 1e-22 + * densityElectrons + * const.value("speed of light in vacuum") + * np.sqrt(1.0 - 1.0 / (1.0 + energyElectron / electronRestMassEnergy) ** 2) + ) # 1/s diff --git a/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/BoundFreeFieldTransitions.py b/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/BoundFreeFieldTransitions.py new file mode 100755 index 0000000000..e84a7e1e78 --- /dev/null +++ b/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/BoundFreeFieldTransitions.py @@ -0,0 +1,127 @@ +#!/usr/bin/env python +""" +atomicPhysics(FLYonPIC) reference rate calculation +This file is part of the PIConGPU. +Copyright 2023-2024 PIConGPU contributors +Authors: Brian Marre, Marco Garten +License: GPLv3+ +""" + +import numpy as np +import scipy.constants as sc +import mpmath as mp + +""" @file reference implementation of the rate calculation for bound-free field based transitions """ + +# set decimal precision of mpmath to be used in calculation +mp.mp.dps = 200 + +# dictionary of atomic units (AU) - values in SI units +atomic_unit = { + "electric field": sc.m_e**2 * sc.e**5 / (sc.hbar**4 * (4 * sc.pi * sc.epsilon_0) ** 3), + "intensity": sc.m_e**4 / (8 * sc.pi * sc.alpha * sc.hbar**9) * sc.e**12 / (4 * sc.pi * sc.epsilon_0) ** 6, + "energy": sc.m_e * sc.e**4 / (sc.hbar**2 * (4 * sc.pi * sc.epsilon_0) ** 2), + "time": sc.hbar**3 * (4 * sc.pi * sc.epsilon_0) ** 2 / sc.m_e / sc.e**4, +} + + +class BoundFreeFieldTransitions: + @staticmethod + def F_crit_BSI(Z, E_Ip): + """Classical barrier suppression field strength. + + :param Z: charge state of the resulting ion + :param E_Ip: ionization potential [unit: AU] + + :returns: critical field strength [unit: AU] + """ + return E_Ip**2 / (4.0 * Z) + + @staticmethod + def n_eff_mpmath(Z, E_Ip): + """Effective principal quantum number as defined in the ADK rate model. + + @detail version using arbitrary precision library mpmath for calculation + + :param Z: charge state of the resulting ion + :param E_Ip: ionization potential [unit: AU] + + :returns: effective principal quantum number + """ + return Z / mp.sqrt(mp.mpf(2) * E_Ip) + + @staticmethod + def n_eff_numpy(Z, E_Ip): + """Effective principal quantum number as defined in the ADK rate model. + + @detail version using fixed numpy functions for calculation, same precision as PIConGPU implementation + + :param Z: charge state of the resulting ion + :param E_Ip: ionization potential [unit: AU] + + :returns: effective principal quantum number + """ + return np.float32(Z / np.sqrt(2.0 * E_Ip)) + + @staticmethod + def ADKRate_numpy(Z, E_Ip, F): + """Ammosov-Delone-Krainov ionization rate. + + A rate model, simplified by Stirling's approximation and setting the + magnetic quantum number m=0 like in publication [DeloneKrainov1998]. + + @detail version using fixed numpy functions for calculation, same precision as PIConGPU implementation + + :param Z: charge state of the resulting ion + :param E_Ip: ionization potential [unit: AU] + :param F: field strength [unit: AU] + + :returns: ionization rate [unit: 1/AU(time)] + """ + + nEff = BoundFreeFieldTransitions.n_eff_numpy(Z, E_Ip) + + dBase = np.float64(4.0 * np.exp(1.0) * Z**3.0 / (F * nEff**4.0)) + D = dBase ** np.float64(nEff) + + print("\t\t nEff: {:.4e}".format(nEff)) + print("\t\t dBase: {:.4e}".format(dBase)) + print("\t\t D: {:.4e}".format(D)) + + gamowFactor = np.exp(-(2.0 * Z**3.0) / (3.0 * nEff**3.0 * F)) + print("\t\t gamowFactor: {:.4e}".format(gamowFactor)) + + rate = ( + F * np.float32(D**2.0 * gamowFactor) * np.sqrt((3.0 * nEff**3.0 * F) / (np.pi * Z**3.0)) / (8.0 * np.pi * Z) + ) + return rate + + @staticmethod + def ADKRate_mpmath(Z, E_Ip, F): + """Ammosov-Delone-Krainov ionization rate. + + A rate model, simplified by Stirling's approximation and setting the + magnetic quantum number m=0 like in publication [DeloneKrainov1998]. + + @detail version using arbitrary precision library mpmath for calculation + + :param Z: charge state of the resulting ion + :param E_Ip: ionization potential [unit: AU] + :param F: field strength [unit: AU] + + :returns: ionization rate [unit: 1/AU(time)] + """ + + nEff = BoundFreeFieldTransitions.n_eff_mpmath(Z, E_Ip) + + dBase = (mp.mpf(4) * mp.e * Z ** mp.mpf(3)) / (F * nEff ** mp.mpf(4)) + D = dBase**nEff + + rate = ( + (F * D ** mp.mpf(2)) + / (mp.mpf(8) * mp.pi * Z) + * mp.exp(-(mp.mpf(2) * Z ** mp.mpf(3)) / (mp.mpf(3) * nEff ** mp.mpf(3) * F)) + * mp.sqrt((mp.mpf(3) * nEff ** mp.mpf(3) * F) / (mp.pi * Z ** mp.mpf(3))) + ) + + return rate diff --git a/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/README.md b/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/README.md new file mode 100644 index 0000000000..d963c63fe8 --- /dev/null +++ b/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/README.md @@ -0,0 +1,11 @@ +FLYonPIC Rate Calculation Reference Implementation +================================================== + +Rate- and cross-section calculation in AtomicPhysics(FLYonPIC) are tested by comparison to reference values in the FLYonPIC RateCalculation unit test. + +The reference values were generated by this python implementation, to regenerate the reference values execute the following: + +```bash + python calculatorMain.py + +``` diff --git a/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/calculatorMain.py b/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/calculatorMain.py new file mode 100755 index 0000000000..ca920a29da --- /dev/null +++ b/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/calculatorMain.py @@ -0,0 +1,217 @@ +#!/usr/bin/env python +""" +atomicPhysics(FLYonPIC) reference rate calculation +This file is part of the PIConGPU. +Copyright 2023-2024 PIConGPU contributors +Authors: Brian Marre +License: GPLv3+ +""" + +import BoundBoundTransitions as boundbound +import BoundFreeCollisionalTransitions as boundfreecollisional +import BoundFreeFieldTransitions as boundfreefield + +import scipy.constants as const +import numpy as np +import mpmath as mp + +""" @file call to get print out of reference rates/crossections """ + +if __name__ == "__main__": + # electron histogram info + # eV + energyElectron = 1000.0 + # eV + energyElectronBinWidth = 10.0 + # 1/(eV * m^3) + densityElectrons = 1e28 + + # bound-free collisional transition data + # eV + ionizationEnergy = 100.0 + # eV + excitationEnergyDifference = 5.0 + # e + screenedCharge = 5.0 + lowerStateLevelVectorBoundFree = (1, 1, 0, 0, 0, 0, 1, 0, 1, 0) + upperStateLevelVectorBoundFree = (1, 1, 0, 0, 0, 0, 1, 0, 0, 0) + + # bound-free field transition data + # screened Charge lower state - 1 + screenedCharge = 4 + # in eV + Hartree = 27.211386245981 + # in Hartree + ionizationEnergy = 5.0 / Hartree + # in atomic_unit["electric field"] ~ 5.1422e11 V/m + fieldStrength = 0.0126 + + # bound-bound transition data + # eV + energyDiffLowerUpper = 5.0 + cxin1 = 1.0 + cxin2 = 2.0 + cxin3 = 3.0 + cxin4 = 4.0 + cxin5 = 5.0 + collisionalOscillatorStrength = 1.0 + absorptionOscillatorStrength = 1.0e-1 + lowerStateLevelVectorBoundBound = (1, 0, 2, 0, 0, 0, 1, 0, 0, 0) + upperStateLevelVectorBoundBound = (1, 0, 1, 0, 0, 0, 1, 0, 1, 0) + # 1/s + frequencyPhoton = energyDiffLowerUpper / const.physical_constants["Planck constant in eV/Hz"][0] + + print("cross sections:") + print("- bound-free") + print( + "\t collisional ionization cross section: \t {0:.12e} 1e6*barn".format( + boundfreecollisional.BoundFreeCollisionalTransitions.collisionalIonizationCrossSection( + energyElectron, + ionizationEnergy, + excitationEnergyDifference, + screenedCharge, + lowerStateLevelVectorBoundFree, + upperStateLevelVectorBoundFree, + ) + ) + ) + + print("- bound-bound") + print( + "\t collisional excitation cross section: \t {0:.12e} 1e6*barn".format( + boundbound.BoundBoundTransitions.collisionalBoundBoundCrossSection( + energyElectron, + energyDiffLowerUpper, + collisionalOscillatorStrength, + cxin1, + cxin2, + cxin3, + cxin4, + cxin5, + lowerStateLevelVectorBoundBound, + upperStateLevelVectorBoundBound, + excitation=True, + ) + ) + ) + print( + "\t collisional deexcitation cross section: {0:.12e} 1e6*barn".format( + boundbound.BoundBoundTransitions.collisionalBoundBoundCrossSection( + energyElectron, + energyDiffLowerUpper, + collisionalOscillatorStrength, + cxin1, + cxin2, + cxin3, + cxin4, + cxin5, + lowerStateLevelVectorBoundBound, + upperStateLevelVectorBoundBound, + excitation=False, + ) + ) + ) + + print("rates:") + print("- bound-free collsional") + print( + "\t collisional ionization rate: \t\t {0:.12e} 1/s".format( + boundfreecollisional.BoundFreeCollisionalTransitions.rateCollisionalIonization( + energyElectron, + energyElectronBinWidth, + densityElectrons, + ionizationEnergy, + excitationEnergyDifference, + screenedCharge, + lowerStateLevelVectorBoundFree, + upperStateLevelVectorBoundFree, + ) + ) + ) + + print("- bound-free field") + # in unit electric field atomic units + nEff = boundfreefield.BoundFreeFieldTransitions.n_eff_numpy(screenedCharge, ionizationEnergy) + fieldStrengthMaxADKRate = 4.0 * screenedCharge**3 / (3.0 * nEff**3 * (4.0 * nEff - 3.0)) + print( + "\t ADK fieldStrength: {0:.4e}, F_maxADK: {1:.4e}, F_crit_BSI: {2:.4e}".format( + fieldStrength, + fieldStrengthMaxADKRate, + boundfreefield.BoundFreeFieldTransitions.F_crit_BSI(screenedCharge, ionizationEnergy), + ) + ) + + print( + "\t ADK rate(numpy) : {0:.9e} * 1/(3.3e-17s)".format( + np.float32( + boundfreefield.BoundFreeFieldTransitions.ADKRate_numpy( + np.float32(screenedCharge), np.float32(ionizationEnergy), np.float32(fieldStrength) + ) + * 3.3e-17 + / boundfreefield.atomic_unit["time"] + ) + ) + ) + print( + "\t ADK rate(mpmath): " + + mp.nstr( + boundfreefield.BoundFreeFieldTransitions.ADKRate_mpmath( + mp.mpf(screenedCharge), mp.mpf(ionizationEnergy), mp.mpf(fieldStrength) + ) + * mp.mpf(3.3e-17) + / mp.mpf(boundfreefield.atomic_unit["time"]), + 10, + max_fixed=1, + ) + + " * 1/(3.3e-17s)" + ) + + print("- bound-bound") + print( + "\t collisional excitation rate: \t\t {0:.12e} 1/s".format( + boundbound.BoundBoundTransitions.rateCollisionalBoundBoundTransition( + energyElectron, + energyElectronBinWidth, + densityElectrons, + energyDiffLowerUpper, + collisionalOscillatorStrength, + cxin1, + cxin2, + cxin3, + cxin4, + cxin5, + lowerStateLevelVectorBoundBound, + upperStateLevelVectorBoundBound, + excitation=True, + ) + ) + ) + print( + "\t collisional deexcitation rate: \t {0:.12e} 1/s".format( + boundbound.BoundBoundTransitions.rateCollisionalBoundBoundTransition( + energyElectron, + energyElectronBinWidth, + densityElectrons, + energyDiffLowerUpper, + collisionalOscillatorStrength, + cxin1, + cxin2, + cxin3, + cxin4, + cxin5, + lowerStateLevelVectorBoundBound, + upperStateLevelVectorBoundBound, + excitation=False, + ) + ) + ) + print( + "\t spontaneous radiative deexcitation rate: {0:.12e} 1/s".format( + boundbound.BoundBoundTransitions.rateSpontaneousDeexcitation( + absorptionOscillatorStrength, + frequencyPhoton, + lowerStateLevelVectorBoundBound, + upperStateLevelVectorBoundBound, + ) + ) + ) diff --git a/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/requirements.txt b/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/requirements.txt new file mode 100644 index 0000000000..103f0aa66f --- /dev/null +++ b/lib/python/picongpu/extra/utils/FLYonPICRateCalculationReference/requirements.txt @@ -0,0 +1,3 @@ +scipy >= 1.12.0 +numpy >= 1.26.2 +mpmath >= 1.3.0 diff --git a/lib/python/picongpu/extra/utils/__init__.py b/lib/python/picongpu/extra/utils/__init__.py index d75b85d0aa..cc9df056e9 100644 --- a/lib/python/picongpu/extra/utils/__init__.py +++ b/lib/python/picongpu/extra/utils/__init__.py @@ -1,6 +1,7 @@ from .find_time import FindTime from .memory_calculator import MemoryCalculator from .field_ionization import FieldIonization +from . import FLYonPICRateCalcualtionReference -__all__ = ["FindTime", "MemoryCalculator", "FieldIonization"] +__all__ = ["FindTime", "MemoryCalculator", "FieldIonization", "FLYonPICRateCalcualtionReference"] diff --git a/lib/python/picongpu/extra/utils/requirements.txt b/lib/python/picongpu/extra/utils/requirements.txt index 24ce15ab7e..5fd814c0af 100644 --- a/lib/python/picongpu/extra/utils/requirements.txt +++ b/lib/python/picongpu/extra/utils/requirements.txt @@ -1 +1,2 @@ numpy +-r FLYonPICRateCalculationReference/requirements.txt