From 5774a38887295a861ffa49428006701ccbdd7fee Mon Sep 17 00:00:00 2001 From: Floor Stikkelbroeck Date: Fri, 6 Dec 2024 13:46:35 +0100 Subject: [PATCH 1/4] update cross section algorithm MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Made the cross section calculation similar to that of López-Puertas et al. 2013 --- amespahdbpythonsuite/transitions.py | 66 +++++++++++++---------------- 1 file changed, 29 insertions(+), 37 deletions(-) diff --git a/amespahdbpythonsuite/transitions.py b/amespahdbpythonsuite/transitions.py index 9c40909..01e7d90 100644 --- a/amespahdbpythonsuite/transitions.py +++ b/amespahdbpythonsuite/transitions.py @@ -923,53 +923,45 @@ def absorption_cross_section(f: np.ndarray) -> np.ndarray: Calculates the PAH absorption cross-section per Li & Draine 2007, ApJ, 657:810-837. - :Params f: frequencies in wavenumber - - :Returns: float array + :param f: Frequencies in wavenumber (can be a single float or np.ndarray). + :return: Array of calculated cross-sections. """ + # Ensure `f` is treated as an array + f = np.atleast_1d(f) + # Convert frequency to wavelength wave = 1e4 / f - A_c = [7.97e-17, 1.23e-17, 20e-21, 14e-21, 80e-24, 84e-24, 46e-24, -322e-24] - W_c = [0.195, 0.217, 0.0805, 0.20, 0.0370, 0.0450, 0.0150, 0.135] - C_c = [0.0722, 0.2175, 1.05, 1.23, 1.66, 1.745, 1.885, 1.90] + # Initialize arrays for different models + cross = np.zeros_like(f) + cross1 = np.zeros_like(f) + cross6 = np.zeros_like(f) - A = np.transpose(np.resize(A_c, (np.size(f), np.size(A_c)))) - W = np.transpose(np.resize(W_c, (np.size(f), np.size(W_c)))) - C = np.transpose(np.resize(C_c, (np.size(f), np.size(C_c)))) + # Set up parameters for Lorentz profile + lamj = np.array([3.3, 6.2, 7.7, 8.6, 11.3, 11.9, 12.7, 16.4, 18.3, 21.2, 23.1, 26.0]) + gammaj = np.array([0.012, 0.032, 0.091, 0.047, 0.018, 0.025, 0.024, 0.010, 0.036, 0.038, 0.046, 0.69]) + aj = np.array([197 * 0.4, 19.6 * 3., 60.9 * 2, 34.7 * 2 * 0.4, 427 / 3. * 0.4, + 72.7 / 3. * 0.4, 167 / 3. * 0.4, 5.52, 6.04, 10.8, 2.78, 15.2]) * 1.0e-20 - # Cutoff wavelength from Salama et al. (1996), over wavelength - # (Eq. (4) in Mattioda et al. (2005)) - y = 1.0 / (0.889 + (2.282 / np.sqrt(0.4 * nc))) / wave + # Calculate Lorentz profile (cross1) for each wavelength + s_values = (2.0 / np.pi) * (gammaj * lamj * aj) / ((wave[:, None] / lamj - lamj / wave[:, None])**2 + gammaj**2) + cross1 = 34.58 * 10**(-18 - 3.431 * wave) + np.sum(s_values, axis=1) - wave_r2 = np.resize(wave, (2, (np.size(f)))) - crosssection = ((1.0 / np.pi) * np.arctan((1e3 * (y - 1.0) ** 3) / y) + 0.5) * ( - 3458e-20 * 10.0 ** (-3.431 * wave) - + (2.0 / np.pi) - * np.sum( - W[:2] - * C[:2] - * A[:2] - / (((wave_r2 / C[:2]) - (C[:2] / wave_r2)) ** 2 + W[:2] ** 2), - axis=0, - ) - ) + # Gaussian model parameters + l0 = np.array([1.05, 1.23, 1.66, 1.745, 1.885, 1.90]) + w = np.array([0.0805, 0.2, 0.037, 0.045, 0.015, 0.135]) + a = np.array([2e-20, 1.4e-20, 8e-23, 8.4e-23, 4.6e-23, 3.22e-22]) - if charge != 0: - wave_r6 = np.resize(wave, (6, np.size(f))) - - crosssection = ( - crosssection - + np.exp(-1e-1 / wave**2) * 1.5e-19 * 10 ** (-wave) - + np.sqrt(2.0 / np.pi) - * np.sum( - A[2:] * np.exp(-2.0 * (wave_r6 - C[2:]) ** 2 / W[2:] ** 2) / W[2:], - axis=0, - ) - ) + # Calculate Gaussian model (cross6) + g_values = (a / (w * np.sqrt(np.pi / 2))) * np.exp(-2 * (wave[:, None] - l0)**2 / w**2) + cross6 = 1.5 * 10**(-19 - wave) + np.sum(g_values, axis=1) - g_values[:, -1] + + # Select final cross-section values based on wavelength threshold + cross = np.where(wave <= 0.8, cross1, cross6) + + return cross - return crosssection @staticmethod def planck(f: np.ndarray) -> np.ndarray: From 5273e65439a6b0345b751c4d0ef1b2168e75660e Mon Sep 17 00:00:00 2001 From: Floor Stikkelbroeck Date: Fri, 10 Jan 2025 13:59:44 +0100 Subject: [PATCH 2/4] fix typos --- amespahdbpythonsuite/transitions.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/amespahdbpythonsuite/transitions.py b/amespahdbpythonsuite/transitions.py index 01e7d90..7b906b0 100644 --- a/amespahdbpythonsuite/transitions.py +++ b/amespahdbpythonsuite/transitions.py @@ -1167,7 +1167,7 @@ def mean_energy(**keywords): def mean_energy_squared(**keywords): """ Callback function to calculate the mean energy in erg^2 for a given - blackbody temperatyre. + blackbody temperature. :Returns: float array @@ -1536,7 +1536,7 @@ def approximate_feature_strength(T: float): def number_of_photons(**keywords): """ Calculate the number of photons given a - blacbody radiation field. + blackbody radiation field. """ From dd3836cd141b2bf7be402d6053acb6f6f569677e Mon Sep 17 00:00:00 2001 From: Floor Stikkelbroeck Date: Wed, 29 Jan 2025 16:14:54 +0100 Subject: [PATCH 3/4] Adjust Absorption_Cross_section adjusted the original code from the nasa ames pahdb database and fixed two small errors such as correct references to Draine & Li 2007 and Li & Draine 2001. implemented the nc>40 requirement for the cutoff frequency, and replaced the charge !=0 with the updated 2008 code. --- amespahdbpythonsuite/transitions.py | 78 ++++++++++++++++++----------- 1 file changed, 48 insertions(+), 30 deletions(-) diff --git a/amespahdbpythonsuite/transitions.py b/amespahdbpythonsuite/transitions.py index 7b906b0..8add085 100644 --- a/amespahdbpythonsuite/transitions.py +++ b/amespahdbpythonsuite/transitions.py @@ -920,47 +920,65 @@ def plot(self, **keywords) -> None: @staticmethod def absorption_cross_section(f: np.ndarray) -> np.ndarray: """ - Calculates the PAH absorption cross-section per Li & Draine 2007, - ApJ, 657:810-837. + Calculates the PAH absorption cross-section per + - Li & Draine 2001b, ApJ, 554:778-802 Eq. [11] + - Draine & 2007, ApJ, 657:810-837. Table 1, same as Li & Draine 2001 - :param f: Frequencies in wavenumber (can be a single float or np.ndarray). - :return: Array of calculated cross-sections. + With NIR corrections from + - Mattioda et al. 2005c + - Mattioda et al. 2008. + + :Params f: frequencies in wavenumber + + :Returns: float array """ - # Ensure `f` is treated as an array - f = np.atleast_1d(f) - # Convert frequency to wavelength wave = 1e4 / f - # Initialize arrays for different models - cross = np.zeros_like(f) - cross1 = np.zeros_like(f) - cross6 = np.zeros_like(f) - - # Set up parameters for Lorentz profile - lamj = np.array([3.3, 6.2, 7.7, 8.6, 11.3, 11.9, 12.7, 16.4, 18.3, 21.2, 23.1, 26.0]) - gammaj = np.array([0.012, 0.032, 0.091, 0.047, 0.018, 0.025, 0.024, 0.010, 0.036, 0.038, 0.046, 0.69]) - aj = np.array([197 * 0.4, 19.6 * 3., 60.9 * 2, 34.7 * 2 * 0.4, 427 / 3. * 0.4, - 72.7 / 3. * 0.4, 167 / 3. * 0.4, 5.52, 6.04, 10.8, 2.78, 15.2]) * 1.0e-20 + #entry 1 and 2 from Draine & Li (2007), (but these are not different from Li & Draine 2001), + #entry 3-5 from Drude Parameters Mattioda et al. (2008) table 3, not used + #entry 6 is one Gaussian Parameters Mattioda et al. (2008) table 3 + C_c = [0.0722, 0.2175, 1.05, 1.26, 1.905, 1.185] + W_c = [0.195, 0.217, 0.055, 0.11, 0.09, 0.2985] + A_c = [7.97e-17, 1.23e-17, 2e-20, 7.8e-21, 1.465e-22, 1.4e-20] - # Calculate Lorentz profile (cross1) for each wavelength - s_values = (2.0 / np.pi) * (gammaj * lamj * aj) / ((wave[:, None] / lamj - lamj / wave[:, None])**2 + gammaj**2) - cross1 = 34.58 * 10**(-18 - 3.431 * wave) + np.sum(s_values, axis=1) + C = np.transpose(np.resize(C_c, (np.size(f), np.size(C_c)))) + W = np.transpose(np.resize(W_c, (np.size(f), np.size(W_c)))) + A = np.transpose(np.resize(A_c, (np.size(f), np.size(A_c)))) - # Gaussian model parameters - l0 = np.array([1.05, 1.23, 1.66, 1.745, 1.885, 1.90]) - w = np.array([0.0805, 0.2, 0.037, 0.045, 0.015, 0.135]) - a = np.array([2e-20, 1.4e-20, 8e-23, 8.4e-23, 4.6e-23, 3.22e-22]) + # Cutoff wavelength from Salama et al. (1996), + # (Eq. (4) in Mattioda et al. (2005c)), adding in the correction for nc > 40 + if nc > 40: + y = 1.0 / (0.889 + (2.282 / np.sqrt(0.4 * nc))) / wave + else: + y = 1.0 / (0.889 + (2.282 / np.sqrt(0.3 * nc))) / wave - # Calculate Gaussian model (cross6) - g_values = (a / (w * np.sqrt(np.pi / 2))) * np.exp(-2 * (wave[:, None] - l0)**2 / w**2) - cross6 = 1.5 * 10**(-19 - wave) + np.sum(g_values, axis=1) - g_values[:, -1] + #continuum from 0.07 - 0.82 μm according to Li & Draine (2001), Eq [11] + wave_r2 = np.resize(wave, (2, (np.size(f)))) + crosssection = ((1.0 / np.pi) * np.arctan((1e3 * (y - 1.0) ** 3) / y) + 0.5) * ( #implement cutoff wavelenght + 3458e-20 * 10.0 ** (-3.431 * wave) #continuum from Li & Draine (2001) Eq. 11 + + (2.0 / np.pi) #Drude profile summation + * np.sum( + W[:2] + * C[:2] + * A[:2] + / (((wave_r2 / C[:2]) - (C[:2] / wave_r2)) ** 2 + W[:2] ** 2), + axis=0, + ) + ) + if charge != 0: + #NIR correction from 0.82 μm according to Mattioda et al. (2008), Eq. (3), table 3 + wave_r6 = np.resize(wave, (6, np.size(f))) - # Select final cross-section values based on wavelength threshold - cross = np.where(wave <= 0.8, cross1, cross6) + crosssection = crosssection + (np.exp(-1e-1 / wave**2) * 1.4e-19 * 10 ** -0.2 * 10 ** (-1.34*wave) + + np.sqrt(2.0 / np.pi) + * np.sum( + A[5:] * np.exp(-2.0 * (wave_r6 - C[5:]) ** 2 / W[5:] ** 2) / W[5:], #Gaussian profile summation + axis=0, + )) - return cross + return crosssection @staticmethod From 5d5bd32e210935f56e63c7246bb758fdf5b568e6 Mon Sep 17 00:00:00 2001 From: Floor Stikkelbroeck Date: Wed, 19 Mar 2025 17:33:24 +0100 Subject: [PATCH 4/4] fix formatting --- amespahdbpythonsuite/transitions.py | 41 ++++++++++++++++------------- 1 file changed, 23 insertions(+), 18 deletions(-) diff --git a/amespahdbpythonsuite/transitions.py b/amespahdbpythonsuite/transitions.py index 8add085..739e98b 100644 --- a/amespahdbpythonsuite/transitions.py +++ b/amespahdbpythonsuite/transitions.py @@ -921,12 +921,12 @@ def plot(self, **keywords) -> None: def absorption_cross_section(f: np.ndarray) -> np.ndarray: """ Calculates the PAH absorption cross-section per - - Li & Draine 2001b, ApJ, 554:778-802 Eq. [11] + - Li & Draine 2001b, ApJ, 554:778-802 Eq. [11] - Draine & 2007, ApJ, 657:810-837. Table 1, same as Li & Draine 2001 - With NIR corrections from - - Mattioda et al. 2005c - - Mattioda et al. 2008. + With NIR corrections from + - Mattioda et al. 2005c + - Mattioda et al. 2008. :Params f: frequencies in wavenumber @@ -935,9 +935,9 @@ def absorption_cross_section(f: np.ndarray) -> np.ndarray: wave = 1e4 / f - #entry 1 and 2 from Draine & Li (2007), (but these are not different from Li & Draine 2001), - #entry 3-5 from Drude Parameters Mattioda et al. (2008) table 3, not used - #entry 6 is one Gaussian Parameters Mattioda et al. (2008) table 3 + # entry 1 and 2 from Draine & Li (2007), (but these are not different from Li & Draine 2001), + # entry 3-5 from Drude Parameters Mattioda et al. (2008) table 3, not used + # entry 6 is one Gaussian Parameters Mattioda et al. (2008) table 3 C_c = [0.0722, 0.2175, 1.05, 1.26, 1.905, 1.185] W_c = [0.195, 0.217, 0.055, 0.11, 0.09, 0.2985] A_c = [7.97e-17, 1.23e-17, 2e-20, 7.8e-21, 1.465e-22, 1.4e-20] @@ -946,19 +946,21 @@ def absorption_cross_section(f: np.ndarray) -> np.ndarray: W = np.transpose(np.resize(W_c, (np.size(f), np.size(W_c)))) A = np.transpose(np.resize(A_c, (np.size(f), np.size(A_c)))) - # Cutoff wavelength from Salama et al. (1996), # (Eq. (4) in Mattioda et al. (2005c)), adding in the correction for nc > 40 - if nc > 40: + if nc > 40: y = 1.0 / (0.889 + (2.282 / np.sqrt(0.4 * nc))) / wave else: y = 1.0 / (0.889 + (2.282 / np.sqrt(0.3 * nc))) / wave - #continuum from 0.07 - 0.82 μm according to Li & Draine (2001), Eq [11] + # continuum from 0.07 - 0.82 μm according to Li & Draine (2001), Eq [11] wave_r2 = np.resize(wave, (2, (np.size(f)))) - crosssection = ((1.0 / np.pi) * np.arctan((1e3 * (y - 1.0) ** 3) / y) + 0.5) * ( #implement cutoff wavelenght - 3458e-20 * 10.0 ** (-3.431 * wave) #continuum from Li & Draine (2001) Eq. 11 - + (2.0 / np.pi) #Drude profile summation + crosssection = ( + (1.0 / np.pi) * np.arctan((1e3 * (y - 1.0) ** 3) / y) + 0.5 + ) * ( # implement cutoff wavelenght + 3458e-20 + * 10.0 ** (-3.431 * wave) # continuum from Li & Draine (2001) Eq. 11 + + (2.0 / np.pi) # Drude profile summation * np.sum( W[:2] * C[:2] @@ -968,19 +970,22 @@ def absorption_cross_section(f: np.ndarray) -> np.ndarray: ) ) if charge != 0: - #NIR correction from 0.82 μm according to Mattioda et al. (2008), Eq. (3), table 3 + # NIR correction from 0.82 μm according to Mattioda et al. (2008), Eq. (3), table 3 wave_r6 = np.resize(wave, (6, np.size(f))) - crosssection = crosssection + (np.exp(-1e-1 / wave**2) * 1.4e-19 * 10 ** -0.2 * 10 ** (-1.34*wave) + crosssection = crosssection + ( + np.exp(-1e-1 / wave**2) * 1.4e-19 * 10**-0.2 * 10 ** (-1.34 * wave) + np.sqrt(2.0 / np.pi) * np.sum( - A[5:] * np.exp(-2.0 * (wave_r6 - C[5:]) ** 2 / W[5:] ** 2) / W[5:], #Gaussian profile summation + A[5:] + * np.exp(-2.0 * (wave_r6 - C[5:]) ** 2 / W[5:] ** 2) + / W[5:], # Gaussian profile summation axis=0, - )) + ) + ) return crosssection - @staticmethod def planck(f: np.ndarray) -> np.ndarray: """