diff --git a/tardis/io/atom_data/base.py b/tardis/io/atom_data/base.py index cf6c10a3394..f991e5ac577 100644 --- a/tardis/io/atom_data/base.py +++ b/tardis/io/atom_data/base.py @@ -44,24 +44,24 @@ class AtomData(object): levels: pandas.DataFrame A DataFrame containing the *levels data* with: - index: no index; + index: numerical index; columns: atomic_number, ion_number, level_number, energy[eV], g[1], metastable. lines: pandas.DataFrame A DataFrame containing the *lines data* with: - index: no index; + index: numerical index; columns: line_id, atomic_number, ion_number, level_number_lower, level_number_upper, wavelength[angstrom], nu[Hz], f_lu[1], f_ul[1], B_ul[?], B_ul[?], A_ul[1/s]. macro_atom_data: A DataFrame containing the *macro atom data* with: - index: no_index; + index: numerical index; columns: atomic_number, ion_number, source_level_number, destination_level_number, transition_line_id, transition_type, transition_probability; macro_atom_references: A DataFrame containing the *macro atom references* with: - index: no_index; + index: numerical index; columns: atomic_number, ion_number, source_level_number, count_down, count_up, count_total. Refer to the docs: http://tardis.readthedocs.io/en/latest/physics/plasma/macroatom.html @@ -74,10 +74,20 @@ class AtomData(object): collision_data_temperatures: np.array An array with the collision temperatures. - zeta_data: ? + zeta_data: + A DataFrame containing the *zeta data* for the + nebular ionization calculation + (i.e., the fraction of recombinations that go directly to the + ground state) with: + index: atomic_number, ion_charge; + columns: temperatures[K] + synpp_refs: ? - ion_cx_tx_data: ? - ion_cx_sp_data: ? + + photoionization_data: pandas.DataFrame + A DataFrame containing the *photoionization data* with: + index: numerical index; + columns: atomic_number, ion_number, level_number, nu[Hz], x_sect[cm^2] Attributes ------------- @@ -91,10 +101,9 @@ class AtomData(object): collision_data_temperatures: numpy.array zeta_data: pandas.DataFrame synpp_refs: pandas.DataFrame - ion_cx_tx_data: pandas.DataFrame - ion_cx_sp_data: pandas.DataFrame symbol2atomic_number: OrderedDict atomic_number2symbol OrderedDict + photoionization_data: pandas.DataFrame Methods -------- @@ -118,8 +127,8 @@ class AtomData(object): "collision_data", "collision_data_temperatures", "synpp_refs", - "ion_cx_th_data", - "ion_cx_sp_data"] + "photoionization_data" + ] # List of tuples of the related dataframes. # Either all or none of the related dataframes must be given @@ -179,12 +188,11 @@ def from_hdf(cls, fname=None): return atom_data - def __init__( - self, atom_data, ionization_data, levels=None, lines=None, - macro_atom_data=None, macro_atom_references=None, - zeta_data=None, collision_data=None, - collision_data_temperatures=None, synpp_refs=None, - ion_cx_th_data=None, ion_cx_sp_data=None): + def __init__(self, atom_data, ionization_data, levels=None, lines=None, + macro_atom_data=None, macro_atom_references=None, + zeta_data=None, collision_data=None, + collision_data_temperatures=None, synpp_refs=None, + photoionization_data=None): self.prepared = False @@ -229,8 +237,8 @@ def __init__( self.collision_data_temperatures = collision_data_temperatures self.synpp_refs = synpp_refs - self.ion_cx_th_data = ion_cx_th_data - self.ion_cx_sp_data = ion_cx_sp_data + + self.photoionization_data = photoionization_data self._check_related() @@ -317,9 +325,6 @@ def prepare_atom_data( self.levels_index.loc[tmp_lines_upper2level_idx]. astype(np.int64).values) - self.atom_ion_index = None - self.levels_index2atom_ion_index = None - if ( self.macro_atom_data_all is not None and not line_interaction_type == 'scatter'): diff --git a/tardis/io/schemas/plasma.yml b/tardis/io/schemas/plasma.yml index 8d368bf372d..56a45371cd6 100644 --- a/tardis/io/schemas/plasma.yml +++ b/tardis/io/schemas/plasma.yml @@ -70,6 +70,17 @@ properties: type: boolean default: false description: sets all beta_sobolevs to 1 + continuum_interaction: + type: object + default: {} + additionalProperties: false + properties: + species: + type: array + default: [] + description: Species that are requested to be treated with continuum + interactios (radiative/collisional ionization and recombination) + in the format ['Si II', 'Ca I', etc.] helium_treatment: type: string default: none diff --git a/tardis/plasma/__init__.py b/tardis/plasma/__init__.py index f7807e3b90d..67591ed64b5 100644 --- a/tardis/plasma/__init__.py +++ b/tardis/plasma/__init__.py @@ -1,2 +1 @@ from tardis.plasma.base import BasePlasma -from tardis.plasma.standard_plasmas import LTEPlasma diff --git a/tardis/plasma/properties/__init__.py b/tardis/plasma/properties/__init__.py index 53899f71981..ae1f36d0be6 100644 --- a/tardis/plasma/properties/__init__.py +++ b/tardis/plasma/properties/__init__.py @@ -7,3 +7,4 @@ from tardis.plasma.properties.radiative_properties import * from tardis.plasma.properties.nlte import * from tardis.plasma.properties.j_blues import * +from tardis.plasma.properties.continuum_processes import * diff --git a/tardis/plasma/properties/atomic.py b/tardis/plasma/properties/atomic.py index 15bb291f3d3..503daba3856 100644 --- a/tardis/plasma/properties/atomic.py +++ b/tardis/plasma/properties/atomic.py @@ -11,7 +11,8 @@ logger = logging.getLogger(__name__) __all__ = ['Levels', 'Lines', 'LinesLowerLevelIndex', 'LinesUpperLevelIndex', - 'AtomicMass', 'IonizationData', 'ZetaData', 'NLTEData'] + 'AtomicMass', 'IonizationData', 'ZetaData', 'NLTEData', + 'PhotoIonizationData'] class Levels(BaseAtomicDataProperty): """ @@ -65,6 +66,48 @@ def _set_index(self, lines): # lines.set_index('line_id', inplace=True) return lines, lines['nu'], lines['f_lu'], lines['wavelength_cm'] + +class PhotoIonizationData(ProcessingPlasmaProperty): + """ + Attributes + ---------- + photo_ion_cross_sections: Pandas DataFrame (nu, x_sect, + index=['atomic_number', + 'ion_number', + 'level_number']), + dtype float) + Table of photoionization cross sections as a + function of frequency. + photo_ion_block_references: One-dimensional Numpy Array, dtype int + Indices where the photoionization data for + a given level starts. Needed for calculation + of recombination rates. + photo_ion_index: Pandas MultiIndex, dtype int + Atomic, ion and level numbers for which + photoionization data exists. + """ + outputs = ('photo_ion_cross_sections', 'photo_ion_block_references', + 'photo_ion_index') + latex_name = ('\\xi_{\\textrm{i}}(\\nu)', '', '') + + def calculate(self, atomic_data, continuum_interaction_species): + photoionization_data = atomic_data.photoionization_data.set_index( + ['atomic_number', 'ion_number', 'level_number'] + ) + selected_species_idx = pd.IndexSlice[ + continuum_interaction_species.get_level_values('atomic_number'), + continuum_interaction_species.get_level_values('ion_number'), + slice(None) + ] + photoionization_data = photoionization_data.loc[selected_species_idx] + phot_nus = photoionization_data['nu'] + block_references = np.hstack( + [[0], phot_nus.groupby(level=[0,1,2]).count().values.cumsum()] + ) + photo_ion_index = photoionization_data.index.unique() + return photoionization_data, block_references, photo_ion_index + + class LinesLowerLevelIndex(HiddenPlasmaProperty): """ Attributes: @@ -92,18 +135,6 @@ def calculate(self, levels, lines): lines_index = lines.index.droplevel('level_number_lower') return np.array(levels_index.loc[lines_index]) -class IonCXData(BaseAtomicDataProperty): - outputs = ('ion_cx_data',) - - def _filter_atomic_property(self, ion_cx_data, selected_atoms): - return ion_cx_data[ion_cx_data.atomic_number.isin([selected_atoms] - if np.isscalar( - selected_atoms) - else selected_atoms)] - - def _set_index(self, ion_cx_data): - return ion_cx_data.set_index(['atomic_number', 'ion_number', - 'level_number']) class AtomicMass(ProcessingPlasmaProperty): """ diff --git a/tardis/plasma/properties/continuum_processes.py b/tardis/plasma/properties/continuum_processes.py new file mode 100644 index 00000000000..5de3f3825e0 --- /dev/null +++ b/tardis/plasma/properties/continuum_processes.py @@ -0,0 +1,92 @@ +import logging + +import numpy as np +import pandas as pd + +from numba import prange, njit +from astropy import constants as const + +from tardis.plasma.properties.base import ProcessingPlasmaProperty + +__all__ = ['SpontRecombRateCoeff'] + +logger = logging.getLogger(__name__) + +njit_dict = {'fastmath': False, 'parallel': False} + + +@njit(**njit_dict) +def integrate_array_by_blocks(f, x, block_references): + """ + Integrates a function f defined at locations x over blocks + given in block_references. + + Parameters + ---------- + f : Two-dimensional Numpy Array, dtype float + x : One-dimensional Numpy Array, dtype float + block_references : One-dimensional Numpy Array, dtype int + + Returns + ------- + integrated : Two-dimensional Numpy Array, dtype float + + """ + integrated = np.zeros((len(block_references) - 1, f.shape[1])) + for i in prange(f.shape[1]): # columns + for j in prange(len(integrated)): # rows + start = block_references[j] + stop = block_references[j + 1] + integrated[j, i] = np.trapz(f[start:stop, i], x[start:stop]) + return integrated + + +def get_ion_multi_index(multi_index_full, next_higher=True): + """ + Integrates a function f defined at locations x over blocks + given in block_references. + + Parameters + ---------- + multi_index_full : Pandas MultiIndex (atomic_number, ion_number, + level_number) + next_higher : bool + If true use ion number of next higher ion, else use ion_number from + multi_index_full. + + Returns + ------- + multi_index : Pandas MultiIndex (atomic_number, ion_number) + + """ + atomic_number = multi_index_full.get_level_values(0) + ion_number = multi_index_full.get_level_values(1) + if next_higher is True: + ion_number += 1 + return pd.MultiIndex.from_arrays([atomic_number, ion_number]) + + +class SpontRecombRateCoeff(ProcessingPlasmaProperty): + """ + Attributes + ---------- + alpha_sp : Pandas DataFrame, dtype float + The rate coefficient for spontaneous recombination. + """ + outputs = ('alpha_sp',) + latex_name = ('\\alpha^{\\textrm{sp}}',) + + def calculate(self, photo_ion_cross_sections, t_electrons, + photo_ion_block_references, photo_ion_index, phi_ik): + x_sect = photo_ion_cross_sections['x_sect'].values + nu = photo_ion_cross_sections['nu'].values + + alpha_sp = (8 * np.pi * x_sect * nu ** 2 / (const.c.cgs.value) ** 2) + alpha_sp = alpha_sp[:, np.newaxis] + boltzmann_factor = np.exp(-nu[np.newaxis].T / t_electrons * + (const.h.cgs.value / const.k_B.cgs.value)) + alpha_sp = alpha_sp * boltzmann_factor + alpha_sp = integrate_array_by_blocks(alpha_sp, nu, + photo_ion_block_references) + alpha_sp = pd.DataFrame(alpha_sp, index=photo_ion_index) + return alpha_sp * phi_ik diff --git a/tardis/plasma/properties/general.py b/tardis/plasma/properties/general.py index 388526d7a1a..cd8717a8d71 100644 --- a/tardis/plasma/properties/general.py +++ b/tardis/plasma/properties/general.py @@ -1,6 +1,7 @@ import logging import numpy as np +import pandas as pd from astropy import units as u from tardis import constants as const @@ -10,7 +11,7 @@ __all__ = ['BetaRadiation', 'GElectron', 'NumberDensity', 'SelectedAtoms', 'ElectronTemperature', 'BetaElectron', 'LuminosityInner', - 'TimeSimulation'] + 'TimeSimulation', 'ThermalGElectron'] class BetaRadiation(ProcessingPlasmaProperty): """ @@ -44,6 +45,22 @@ def calculate(self, beta_rad): return ((2 * np.pi * const.m_e.cgs.value / beta_rad) / (const.h.cgs.value ** 2)) ** 1.5 + +class ThermalGElectron(GElectron): + """ + Attributes + ---------- + thermal_g_electron : Numpy Array, dtype float + """ + outputs = ('thermal_g_electron',) + latex_name = ('g_{\\textrm{electron_thermal}}',) + latex_formula = ('\\Big(\\dfrac{2\\pi m_{e}/\ + \\beta_{\\textrm{electron}}}{h^2}\\Big)^{3/2}',) + + def calculate(self, beta_electron): + return super(ThermalGElectron, self).calculate(beta_electron) + + class NumberDensity(ProcessingPlasmaProperty): """ Attributes @@ -63,7 +80,7 @@ class SelectedAtoms(ProcessingPlasmaProperty): """ Attributes ---------- - selected_atoms : Numpy Array, dtype int + selected_atoms : Pandas Int64Index, dtype int Atomic numbers of elements required for particular simulation """ outputs = ('selected_atoms',) diff --git a/tardis/plasma/properties/ion_population.py b/tardis/plasma/properties/ion_population.py index d4ac016f764..a98a957997f 100644 --- a/tardis/plasma/properties/ion_population.py +++ b/tardis/plasma/properties/ion_population.py @@ -1,6 +1,7 @@ import logging import warnings +import sys import numpy as np import pandas as pd import numexpr as ne @@ -8,13 +9,15 @@ from scipy import interpolate from tardis.plasma.properties.base import ProcessingPlasmaProperty +from tardis.plasma.properties.continuum_processes import get_ion_multi_index from tardis.plasma.exceptions import PlasmaIonizationError logger = logging.getLogger(__name__) __all__ = ['PhiSahaNebular', 'PhiSahaLTE', 'RadiationFieldCorrection', - 'IonNumberDensity', 'IonNumberDensityHeNLTE'] + 'IonNumberDensity', 'IonNumberDensityHeNLTE', 'SahaFactor', + 'ThermalPhiSahaLTE'] def calculate_block_ids_from_dataframe(dataframe): @@ -27,7 +30,8 @@ class PhiSahaLTE(ProcessingPlasmaProperty): """ Attributes: phi : Pandas DataFrame, dtype float - Used for LTE ionization. Indexed by atomic number, ion number. Columns are zones. + Used for LTE ionization (at the radiation temperature). + Indexed by atomic number, ion number. Columns are zones. """ outputs = ('phi',) latex_name = ('\\Phi',) @@ -70,6 +74,28 @@ def _calculate_block_ids(partition_function): partition_function.index.get_level_values(0).unique() +class ThermalPhiSahaLTE(PhiSahaLTE): + """ + Attributes: + phi : Pandas DataFrame, dtype float + Used for LTE ionization (at the electron temperature). + Indexed by atomic number, ion number. Columns are zones. + """ + outputs = ('thermal_phi_lte',) + latex_name = ('\\Phi^{*}(T_\\mathrm{e})',) + latex_formula = ('\\dfrac{2Z_{i,j+1}}{Z_{i,j}}\\Big(\ + \\dfrac{2\\pi m_{e}/\\beta_{\\textrm{electron}}}{h^2}\ + \\Big)^{3/2}e^{\\dfrac{-\\chi_{i,j}}{kT_{\ + \\textrm{electron}}}}',) + + @staticmethod + def calculate(thermal_g_electron, beta_electron, + thermal_lte_partition_function, ionization_data): + return super(ThermalPhiSahaLTE, ThermalPhiSahaLTE).calculate( + thermal_g_electron, beta_electron, thermal_lte_partition_function, + ionization_data + ) + class PhiSahaNebular(ProcessingPlasmaProperty): """ @@ -350,3 +376,40 @@ def calculate(self, phi, partition_function, number_density, ion_number_density.loc[2, 2].update(helium_population_updated.loc[ 2, 0]) return ion_number_density, n_electron, helium_population_updated + + +class SahaFactor(ProcessingPlasmaProperty): + """ + Calculates the 'Saha factor' Phi_ik = n_i* / (n_k* n_e), i.e., + the ratio of the LTE level population n_i*, and the product of + the LTE ion density n_k* and the actual electron density n_e. + + Attributes: + phi_ik: Pandas DataFrame, dtype float + Indexed by atom number, ion number, level number. + Columns are zones. + """ + outputs = ('phi_ik',) + latex_name = ('\\Phi_{i,\\kappa}',) + + def calculate(self, thermal_phi_lte, thermal_lte_level_boltzmann_factor, + thermal_lte_partition_function): + boltzmann_factor = self._prepare_boltzmann_factor( + thermal_lte_level_boltzmann_factor + ) + phi_saha_index = get_ion_multi_index(boltzmann_factor.index) + partition_function_index = get_ion_multi_index(boltzmann_factor.index, + next_higher=False) + phi_saha = thermal_phi_lte.loc[phi_saha_index].values + # Replace zero values in phi_saha to avoid zero division in Saha factor + phi_saha[phi_saha == 0.0] = sys.float_info.min + partition_function = thermal_lte_partition_function.loc[ + partition_function_index].values + return boltzmann_factor / (phi_saha * partition_function) + + @staticmethod + def _prepare_boltzmann_factor(boltzmann_factor): + atomic_number = boltzmann_factor.index.get_level_values(0) + ion_number = boltzmann_factor.index.get_level_values(1) + selected_ions_mask = (atomic_number != ion_number) + return boltzmann_factor[selected_ions_mask] diff --git a/tardis/plasma/properties/partition_function.py b/tardis/plasma/properties/partition_function.py index cb209a123b2..0b31d3c6966 100644 --- a/tardis/plasma/properties/partition_function.py +++ b/tardis/plasma/properties/partition_function.py @@ -11,7 +11,8 @@ __all__ = ['LevelBoltzmannFactorLTE', 'LevelBoltzmannFactorDiluteLTE', 'LevelBoltzmannFactorNoNLTE', 'LevelBoltzmannFactorNLTE', - 'PartitionFunction'] + 'PartitionFunction', 'ThermalLevelBoltzmannFactorLTE', + 'ThermalLTEPartitionFunction'] class LevelBoltzmannFactorLTE(ProcessingPlasmaProperty): @@ -19,8 +20,9 @@ class LevelBoltzmannFactorLTE(ProcessingPlasmaProperty): Attributes ---------- general_level_boltzmann_factor : Pandas DataFrame, dtype float - Level population proportionality values. Indexed - by atomic number, ion number, level number. + Level population proportionality values. + Evaluated at the radiation temperature. + Indexed by atomic number, ion number, level number. Columns corresponding to zones. Does not consider NLTE. """ @@ -41,6 +43,30 @@ def calculate(excitation_energy, g, beta_rad, levels): return level_boltzmann_factor +class ThermalLevelBoltzmannFactorLTE(LevelBoltzmannFactorLTE): + """ + Attributes + ---------- + thermal_lte_level_boltzmann_factor : Pandas DataFrame, dtype float + Level population proportionality values for LTE. + Evaluated at the temperature of the + electron gas (thermal). Indexed + by atomic number, ion number, level number. + Columns corresponding to zones. + """ + outputs = ('thermal_lte_level_boltzmann_factor',) + latex_name = ('bf_{i,j,k}^{\\textrm{LTE}}(T_e)',) + latex_formula = ('g_{i,j,k}e^{\\dfrac{-\\epsilon_{i,j,k}}{k_{\ + \\textrm{B}}T_{\\textrm{electron}}}}',) + + @staticmethod + def calculate(excitation_energy, g, beta_electron, levels): + return super(ThermalLevelBoltzmannFactorLTE, + ThermalLevelBoltzmannFactorLTE).calculate( + excitation_energy, g, beta_electron, levels + ) + + class LevelBoltzmannFactorDiluteLTE(ProcessingPlasmaProperty): """ Attributes @@ -275,3 +301,20 @@ class PartitionFunction(ProcessingPlasmaProperty): def calculate(self, level_boltzmann_factor): return level_boltzmann_factor.groupby( level=['atomic_number', 'ion_number']).sum() + + +class ThermalLTEPartitionFunction(PartitionFunction): + """ + Attributes + ---------- + thermal_lte_partition_function : Pandas DataFrame, dtype float + Indexed by atomic number, ion number. + Columns are zones. + """ + outputs = ('thermal_lte_partition_function',) + latex_name = ('Z_{i,j}(T_\\mathrm{e}',) + + def calculate(self, thermal_lte_level_boltzmann_factor): + return super(ThermalLTEPartitionFunction, self).calculate( + thermal_lte_level_boltzmann_factor + ) diff --git a/tardis/plasma/properties/plasma_input.py b/tardis/plasma/properties/plasma_input.py index 0c27fc96bf1..a9e0cb90f55 100644 --- a/tardis/plasma/properties/plasma_input.py +++ b/tardis/plasma/properties/plasma_input.py @@ -2,7 +2,9 @@ __all__ = ['TRadiative', 'DilutionFactor', 'AtomicData', 'Abundance', 'Density', 'TimeExplosion', 'JBlueEstimator', 'LinkTRadTElectron', - 'HeliumTreatment', 'RInner', 'TInner', 'Volume'] + 'HeliumTreatment', 'RInner', 'TInner', 'Volume', + 'ContinuumInteractionSpecies'] + class TRadiative(ArrayInput): """ @@ -13,6 +15,7 @@ class TRadiative(ArrayInput): outputs = ('t_rad',) latex_name = ('T_{\\textrm{rad}}',) + class DilutionFactor(ArrayInput): """ Attributes @@ -24,6 +27,7 @@ class DilutionFactor(ArrayInput): outputs = ('w',) latex_name = ('W',) + class AtomicData(Input): """ Attributes @@ -42,6 +46,7 @@ class Abundance(Input): """ outputs = ('abundance',) + class Density(ArrayInput): """ Attributes @@ -52,6 +57,7 @@ class Density(ArrayInput): outputs = ('density',) latex_name = ('\\rho',) + class TimeExplosion(Input): """ Attributes @@ -72,6 +78,7 @@ class JBlueEstimator(ArrayInput): outputs = ('j_blue_estimators',) latex_name = ('J_{\\textrm{blue-estimator}}',) + class LinkTRadTElectron(Input): """ Attributes @@ -83,14 +90,29 @@ class LinkTRadTElectron(Input): outputs = ('link_t_rad_t_electron',) latex_name = ('T_{\\textrm{electron}}/T_{\\textrm{rad}}',) + class HeliumTreatment(Input): outputs = ('helium_treatment',) + class RInner(Input): outputs = ('r_inner',) + class TInner(Input): outputs = ('t_inner',) + class Volume(Input): outputs = ('volume',) + + +class ContinuumInteractionSpecies(Input): + """ + Attributes + ---------- + continuum_interaction_species: Pandas MultiIndex, dtype int + Atomic and ion numbers of elements for which continuum interactions + (radiative/collisional ionization and recombination) are treated + """ + outputs = ('continuum_interaction_species',) diff --git a/tardis/plasma/properties/property_collections.py b/tardis/plasma/properties/property_collections.py index cb6c0915862..3f943270b3f 100644 --- a/tardis/plasma/properties/property_collections.py +++ b/tardis/plasma/properties/property_collections.py @@ -5,7 +5,7 @@ class PlasmaPropertyCollection(list): basic_inputs = PlasmaPropertyCollection([TRadiative, Abundance, Density, TimeExplosion, AtomicData, DilutionFactor, LinkTRadTElectron, - HeliumTreatment]) + HeliumTreatment, ContinuumInteractionSpecies]) basic_properties = PlasmaPropertyCollection([BetaRadiation, Levels, Lines, AtomicMass, PartitionFunction, GElectron, IonizationData, NumberDensity, LinesLowerLevelIndex, @@ -36,3 +36,8 @@ class PlasmaPropertyCollection(list): JBluesNormFactor, LuminosityInner, TimeSimulation]) +continuum_interaction_properties = PlasmaPropertyCollection( + [PhotoIonizationData, SpontRecombRateCoeff, + ThermalLevelBoltzmannFactorLTE, ThermalLTEPartitionFunction, BetaElectron, + ThermalGElectron, ThermalPhiSahaLTE, SahaFactor] +) diff --git a/tardis/plasma/standard_plasmas.py b/tardis/plasma/standard_plasmas.py index 5ac131f7045..08a3cae1070 100644 --- a/tardis/plasma/standard_plasmas.py +++ b/tardis/plasma/standard_plasmas.py @@ -1,9 +1,11 @@ import os import logging +import numpy as np import pandas as pd from tardis.io.atom_data import AtomData +from tardis.io.config_reader import ConfigurationError from tardis.util.base import species_string_to_tuple from tardis.plasma import BasePlasma from tardis.plasma.properties.property_collections import (basic_inputs, @@ -11,7 +13,8 @@ macro_atom_properties, dilute_lte_excitation_properties, nebular_ionization_properties, non_nlte_properties, nlte_properties, helium_nlte_properties, helium_numerical_nlte_properties, - helium_lte_properties, detailed_j_blues_properties, detailed_j_blues_inputs) + helium_lte_properties, detailed_j_blues_properties, + detailed_j_blues_inputs, continuum_interaction_properties) from tardis.plasma.exceptions import PlasmaConfigError from tardis.plasma.properties import ( @@ -26,20 +29,6 @@ logger = logging.getLogger(__name__) -class LTEPlasma(BasePlasma): - - def __init__(self, t_rad, abundance, density, time_explosion, atomic_data, - j_blues, link_t_rad_t_electron=0.9, delta_treatment=None): - plasma_modules = basic_inputs + basic_properties + \ - lte_excitation_properties + lte_ionization_properties + \ - non_nlte_properties - - super(LTEPlasma, self).__init__(plasma_properties=plasma_modules, - t_rad=t_rad, abundance=abundance, atomic_data=atomic_data, - density=density, time_explosion=time_explosion, j_blues=j_blues, - w=None, link_t_rad_t_electron=link_t_rad_t_electron, - delta_input=delta_treatment, nlte_species=None) - def assemble_plasma(config, model, atom_data=None): """ @@ -63,6 +52,15 @@ def assemble_plasma(config, model, atom_data=None): nlte_species = [species_string_to_tuple(s) for s in config.plasma.nlte.species] + # Convert the continuum interaction species list to a proper format. + continuum_interaction_species = [ + species_string_to_tuple(s) for s in + config.plasma.continuum_interaction.species + ] + continuum_interaction_species = pd.MultiIndex.from_tuples( + continuum_interaction_species, names=['atomic_number', 'ion_number'] + ) + if atom_data is None: if 'atom_data' in config: if os.path.isabs(config.atom_data): @@ -88,13 +86,28 @@ def assemble_plasma(config, model, atom_data=None): line_interaction_type=config.plasma.line_interaction_type, nlte_species=nlte_species) + # Check if continuum interaction species are in selected_atoms + continuum_atoms = continuum_interaction_species.get_level_values( + 'atomic_number' + ) + continuum_atoms_in_selected_atoms = np.all( + continuum_atoms.isin(atom_data.selected_atomic_numbers) + ) + if not continuum_atoms_in_selected_atoms: + raise ConfigurationError('Not all continuum interaction species ' + 'belong to atoms that have been specified ' + 'in the configuration.') + kwargs = dict(t_rad=model.t_radiative, abundance=model.abundance, density=model.density, atomic_data=atom_data, time_explosion=model.time_explosion, - w=model.dilution_factor, link_t_rad_t_electron=0.9) + w=model.dilution_factor, link_t_rad_t_electron=0.9, + continuum_interaction_species=continuum_interaction_species) plasma_modules = basic_inputs + basic_properties property_kwargs = {} + if config.plasma.continuum_interaction.species: + plasma_modules += continuum_interaction_properties if config.plasma.radiative_rates_type == 'blackbody': plasma_modules.append(JBluesBlackBody) elif config.plasma.radiative_rates_type == 'dilute-blackbody':