diff --git a/tardis/plasma/properties/atomic.py b/tardis/plasma/properties/atomic.py index 49caac342f6..fcbef7bcbcb 100644 --- a/tardis/plasma/properties/atomic.py +++ b/tardis/plasma/properties/atomic.py @@ -10,7 +10,7 @@ logger = logging.getLogger(__name__) __all__ = ['Levels', 'Lines', 'LinesLowerLevelIndex', 'LinesUpperLevelIndex', - 'AtomicMass', 'IonizationData', 'ZetaData', 'NLTEData'] + 'AtomicMass', 'IonizationData', 'ZetaData', 'NLTEData', 'Chi0'] class BaseAtomicDataProperty(ProcessingPlasmaProperty): __metaclass__ = ABCMeta @@ -232,3 +232,9 @@ def calculate(self, atomic_data): return (getattr(self, self.outputs[0]),) else: return atomic_data.nlte_data + +class Chi0(ProcessingPlasmaProperty): + outputs = ('chi_0',) + + def calculate(self, atomic_data): + return atomic_data.ionization_data.ionization_energy.ix[20].ix[2] diff --git a/tardis/plasma/properties/ion_population.py b/tardis/plasma/properties/ion_population.py index dee616b9f43..74e10d11f0e 100644 --- a/tardis/plasma/properties/ion_population.py +++ b/tardis/plasma/properties/ion_population.py @@ -94,21 +94,17 @@ class RadiationFieldCorrection(ProcessingPlasmaProperty): outputs = ('delta',) latex_name = ('\\delta',) - def __init__(self, plasma_parent, departure_coefficient=None, - chi_0_species=(20,2)): + def __init__(self, plasma_parent, departure_coefficient=None): super(RadiationFieldCorrection, self).__init__(plasma_parent) self.departure_coefficient = departure_coefficient - self.chi_0_species = chi_0_species def calculate(self, w, ionization_data, beta_rad, t_electrons, t_rad, - beta_electron, delta_input): + beta_electron, delta_input, chi_0): if delta_input is None: if self.departure_coefficient is None: departure_coefficient = 1. / w else: departure_coefficient = self.departure_coefficient - chi_0_species=self.chi_0_species - chi_0 = ionization_data.ionization_energy.ix[chi_0_species] radiation_field_correction = -np.ones((len(ionization_data), len( beta_rad))) less_than_chi_0 = ( diff --git a/tardis/plasma/properties/property_collections.py b/tardis/plasma/properties/property_collections.py index f4a73fbdc6f..cb7a8b769a9 100644 --- a/tardis/plasma/properties/property_collections.py +++ b/tardis/plasma/properties/property_collections.py @@ -10,7 +10,7 @@ RadiationFieldCorrection, RadiationFieldCorrectionInput, LevelBoltzmannFactorNoNLTE, LevelBoltzmannFactorNLTE, NLTEData, NLTESpecies, PreviousBetaSobolevs, LTEJBlues, - PreviousElectronDensities) + PreviousElectronDensities, Chi0) class PlasmaPropertyCollection(list): pass @@ -29,7 +29,7 @@ class PlasmaPropertyCollection(list): macro_atom_properties = PlasmaPropertyCollection([BetaSobolev, TransitionProbabilities]) nebular_ionization_properties = PlasmaPropertyCollection([PhiSahaNebular, - ZetaData, BetaElectron, RadiationFieldCorrection]) + ZetaData, BetaElectron, RadiationFieldCorrection, Chi0]) dilute_lte_excitation_properties = PlasmaPropertyCollection([ LevelBoltzmannFactorDiluteLTE]) non_nlte_properties = PlasmaPropertyCollection([LevelBoltzmannFactorNoNLTE]) diff --git a/tardis/plasma/standard_plasmas.py b/tardis/plasma/standard_plasmas.py index 15f2b567ae3..b0b789679d3 100644 --- a/tardis/plasma/standard_plasmas.py +++ b/tardis/plasma/standard_plasmas.py @@ -8,6 +8,7 @@ macro_atom_properties, dilute_lte_excitation_properties, nebular_ionization_properties, non_nlte_properties, nlte_properties) +from tardis.io.util import parse_abundance_dict_to_dataframe logger = logging.getLogger(__name__) @@ -47,7 +48,7 @@ def initial_w(self, number_densities): def update_radiationfield(self, t_rad, ws, j_blues, nlte_config, t_electrons=None, n_e_convergence_threshold=0.05, initialize_nlte=False): - if nlte_config.species: + if nlte_config is not None and nlte_config.species: self.store_previous_properties() self.update(t_rad=t_rad, w=ws, j_blues=j_blues) @@ -74,7 +75,7 @@ def __init__(self, number_densities, atomic_data, time_explosion, raise NotImplementedError('Sorry ' + ionization_mode + ' not implemented yet.') - if nlte_config.species: + if nlte_config is not None and nlte_config.species: plasma_modules += nlte_properties else: plasma_modules += non_nlte_properties @@ -90,16 +91,23 @@ def __init__(self, number_densities, atomic_data, time_explosion, abundance, density = self.from_number_densities(number_densities, atomic_data) - initial_beta_sobolevs = np.ones((len(atomic_data.lines), - len(number_densities.columns))) + try: + initial_beta_sobolevs = np.ones((len(atomic_data.lines), + len(number_densities.columns))) + except: + initial_beta_sobolevs = np.ones((len(atomic_data._lines), + len(number_densities.columns))) initial_electron_densities = number_densities.sum(axis=0) - self.nlte_config = nlte_config + if nlte_config is not None and nlte_config.species: + nlte_species = nlte_config.species + else: + nlte_species = None super(LegacyPlasmaArray, self).__init__(plasma_properties=plasma_modules, t_rad=t_rad, abundance=abundance, density=density, atomic_data=atomic_data, time_explosion=time_explosion, j_blues=None, w=w, link_t_rad_t_electron=link_t_rad_t_electron, - delta_input=delta_treatment, nlte_species=nlte_config.species, + delta_input=delta_treatment, nlte_species=nlte_species, previous_electron_densities=initial_electron_densities, previous_beta_sobolevs=initial_beta_sobolevs) diff --git a/tardis/plasma/tests/conftest.py b/tardis/plasma/tests/conftest.py index 0df4284a63f..4aa8f703a71 100644 --- a/tardis/plasma/tests/conftest.py +++ b/tardis/plasma/tests/conftest.py @@ -14,7 +14,8 @@ from tardis.plasma.properties.partition_function import ( LevelBoltzmannFactorLTE, PartitionFunction, LevelBoltzmannFactorDiluteLTE) from tardis.plasma.properties.atomic import (Levels, Lines, AtomicMass, - IonizationData, LinesUpperLevelIndex, LinesLowerLevelIndex, ZetaData) + IonizationData, LinesUpperLevelIndex, LinesLowerLevelIndex, ZetaData, + Chi0) from tardis.plasma.properties.level_population import LevelNumberDensity from tardis.plasma.properties.radiative_properties import (TauSobolev, StimulatedEmissionFactor, BetaSobolev, TransitionProbabilities) @@ -165,6 +166,11 @@ def lines_lower_level_index(lines, levels): lower_level_index_module = LinesLowerLevelIndex(None) return lower_level_index_module.calculate(levels, lines) +@pytest.fixture +def chi_0(atomic_data): + chi_0_module = Chi0(None) + return chi_0_module.calculate(atomic_data) + # PARTITION FUNCTION PROPERTIES @pytest.fixture @@ -223,12 +229,12 @@ def electron_densities(phi_saha_lte, partition_function, number_density): @pytest.fixture def delta(w, ionization_data, beta_rad, t_electrons, t_rad, beta_electron, - levels): + levels, chi_0): delta_input = None delta_module = RadiationFieldCorrection(None) delta_module.chi_0_species = (2, 2) return delta_module.calculate(w, ionization_data, beta_rad, t_electrons, - t_rad, beta_electron, delta_input) + t_rad, beta_electron, delta_input, chi_0) # LEVEL POPULATION PROPERTIES diff --git a/tardis/plasma/tests/test_property_ion_population.py b/tardis/plasma/tests/test_property_ion_population.py index 675ca4d3f67..6c3bfc411b5 100644 --- a/tardis/plasma/tests/test_property_ion_population.py +++ b/tardis/plasma/tests/test_property_ion_population.py @@ -20,7 +20,8 @@ def test_electron_densities(electron_densities): assert np.allclose(electron_densities, 1.181197e+09) def test_radiation_field_correction(delta): + print delta.ix[2].ix[2] assert np.allclose(delta.ix[2].ix[2], 0.000807200897) def test_phi_saha_nebular(phi_saha_nebular): - assert np.allclose(phi_saha_nebular.ix[2].ix[1], 1449588387.9727359) + assert np.allclose(phi_saha_nebular.ix[2].ix[1], 54784028.964375) diff --git a/tardis/tests/test_plasma_simple.py b/tardis/tests/test_plasma_simple.py index 5197d65166d..a06d50f7ceb 100644 --- a/tardis/tests/test_plasma_simple.py +++ b/tardis/tests/test_plasma_simple.py @@ -1,25 +1,37 @@ from astropy import constants as const, units as u import os +import pandas as pd import tardis from tardis import atomic import pytest +from tardis.plasma.standard_plasmas import LegacyPlasmaArray +from tardis.io.util import parse_abundance_dict_to_dataframe #from numpy.testing import assert_allclose data_path = os.path.join(tardis.__path__[0], 'tests', 'data') helium_test_db = os.path.join(data_path, 'chianti_he_db.h5') -''' class TestNebularPlasma(object): def setup(self): atom_data = atomic.AtomData.from_hdf5(helium_test_db) - self.plasma = plasma_array.BasePlasmaArray.from_abundance( - {'He':1.0}, 1e-15*u.Unit('g/cm3'), atom_data, 10 * u.day, - ionization_mode='nebular', excitation_mode='dilute-lte') + density = 1e-15 * u.Unit('g/cm3') + abundance = parse_abundance_dict_to_dataframe({'He':1.0}) + abundance = pd.DataFrame({0:abundance}) + number_densities = abundance * density.to('g/cm^3').value + number_densities = number_densities.div( + atom_data.atom_data.mass.ix[number_densities.index], axis=0) + self.plasma = LegacyPlasmaArray(number_densities, + atomic_data=atom_data, time_explosion=10 * u.day, + ionization_mode='nebular') + +# self.plasma = plasma_array.BasePlasmaArray.from_abundance( +# {'He':1.0}, 1e-15*u.Unit('g/cm3'), atom_data, 10 * u.day, +# ionization_mode='nebular', excitation_mode='dilute-lte') def test_high_temperature(self): with pytest.raises(ValueError) as excinfo: - self.plasma.update_radiationfield([100000.], [1.]) + self.plasma.update_radiationfield(t_rad=[100000.], + ws=[0.5], j_blues=None, nlte_config=None) assert str(excinfo.value).startswith('t_rads outside of zeta ' - 'factor interpolation') -''' \ No newline at end of file + 'factor interpolation') \ No newline at end of file