From 8d3c47affd68f01501bae0de9b5ddb3a62aa82fe Mon Sep 17 00:00:00 2001 From: "Michael.Klauser" Date: Thu, 7 Mar 2013 17:25:13 +0100 Subject: [PATCH 1/4] add bf calc pra yaml merge --- tardis/atomic.py | 34 ++++++++++++++++++++++++++++++++-- tardis/config_reader.py | 1 - tardis/plasma.py | 23 +++++++++++++++++++++++ 3 files changed, 55 insertions(+), 3 deletions(-) diff --git a/tardis/atomic.py b/tardis/atomic.py index 6b983b2830b..59900e34fe7 100644 --- a/tardis/atomic.py +++ b/tardis/atomic.py @@ -214,6 +214,17 @@ def read_collision_data(fname): return collision_data, collision_temperatures +def read_ion_cx_data(fname): + try: + h5_file = h5py.File(fname) + ion_cx_th_data = h5_file['ionization_cx_threshold'] + ion_cx_sp_data = h5_file['ionization_cx_support'] + return ion_cx_th_data, ion_cx_sp_data + except IOError, err: + print(err.errno) + print(err) + logger.critical('Cannot import. Error opening the file to read ionization_cx') + def read_macro_atom_data(fname): if fname is None: @@ -317,9 +328,15 @@ def from_hdf5(cls, fname=None): else: synpp_refs = None + + if 'ion_cx_data' in h5_datasets and 'ion_cx_data' in h5_datasets: + ion_cx_data = read_ion_cx_data(fname) + else: + ion_cx_data = None + atom_data = cls(atom_data=atom_data, ionization_data=ionization_data, levels_data=levels_data, lines_data=lines_data, macro_atom_data=macro_atom_data, zeta_data=zeta_data, - collision_data=(collision_data, collision_data_temperatures), synpp_refs=synpp_refs) + collision_data=(collision_data, collision_data_temperatures), synpp_refs=synpp_refs, ion_cx_data = ion_cx_data) with h5py.File(fname) as h5_file: atom_data.uuid1 = h5_file.attrs['uuid1'] @@ -328,7 +345,7 @@ def from_hdf5(cls, fname=None): return atom_data def __init__(self, atom_data, ionization_data, levels_data, lines_data, macro_atom_data=None, zeta_data=None, - collision_data=None, synpp_refs=None): + collision_data=None, synpp_refs=None,ion_cx_data = None): if macro_atom_data is not None: @@ -339,6 +356,17 @@ def __init__(self, atom_data, ionization_data, levels_data, lines_data, macro_at else: self.has_macro_atom = False + if ion_cx_data is not None: + self.has_ion_cx_data = True + #TODO:Farm a panda here + self.ion_cx_th_data = DataFrame(np.array(ion_cx_data[0])) + self.ion_cx_th_data.set_index(['atomic_number', 'ion_number', 'level_id'], inplace=True) + + self.ion_cx_sp_data = DataFrame(np.array(ion_cx_data[1])) + self.ion_cx_sp_data.set_index(['atomic_number', 'ion_number', 'level_id']) + else: + self.has_ion_cx_data = False + if zeta_data is not None: self.zeta_data = zeta_data self.has_zeta_data = True @@ -351,6 +379,8 @@ def __init__(self, atom_data, ionization_data, levels_data, lines_data, macro_at self.collision_data.set_index(['atomic_number', 'ion_number', 'level_number_lower', 'level_number_upper'], inplace=True) + + self.has_collision_data = True else: self.has_collision_data = False diff --git a/tardis/config_reader.py b/tardis/config_reader.py index ab6bc6f9b28..507fced2bd7 100644 --- a/tardis/config_reader.py +++ b/tardis/config_reader.py @@ -214,7 +214,6 @@ def parse_general_section(self, config_dict): self.velocities, self.time_explosion) elif density_set == 'exponential': - #TODO:Add here the function call which generates the exponential density profile. The easy way from tonight don't work as expected!! if not (('exponential_n_factor' in config_dict) and ('exponential_rho0' in config_dict)): raise ValueError( 'If density=exponential is set the exponential_n_factor(float) and exponential_rho_0 have to be specified.') diff --git a/tardis/plasma.py b/tardis/plasma.py index 2b75b503a1b..e08bf290eef 100644 --- a/tardis/plasma.py +++ b/tardis/plasma.py @@ -437,6 +437,29 @@ def calculate_tau_sobolev(self): n_lower = self.level_populations.values[self.atom_data.lines_lower2level_idx] self.tau_sobolevs = sobolev_coefficient * f_lu * wavelength * self.time_explosion * n_lower + def calculate_bound_free(self): + """ + :return: + """ + nu_bins = range(1000, 10000, 1000) #TODO: get the binning from the input file. + try: + bf = np.zeros(len(self.atom_data.levels), len(self.atom_data.selected_atomic_numbers), len(nu_bins)) + except AttributeError: + logger.critical("Err creating the bf array.") + + phis = self.calculate_saha() + nnlevel = self.level_populations + for nu in nu_bins: + for i,(level_id,level) in enumerate(self.atom_data.levels.iterrows()): + atomic_number = level.name[0] + ion_number = level.name[1] + level_number = level.name[2] + sigma_bf_th = self.atom_data.ion_cx_th.ix[atomic_number,ion_number,level_number] + phi = phis.ix[atomic_number,ion_number] + + + + def update_macro_atom(self): """ Updating the Macro Atom computations From 2ff5228f9b864105a222b6b15c027705c08ea87c Mon Sep 17 00:00:00 2001 From: "Michael.Klauser" Date: Thu, 7 Mar 2013 20:12:27 +0100 Subject: [PATCH 2/4] Merge remote-tracking branch 'remotes/origin/master' into Bound-Free and add the exponential. Conflicts: tardis/config_reader.py --- tardis/config_reader.py | 107 ++++++++++++++++++++++++++++------------ 1 file changed, 76 insertions(+), 31 deletions(-) diff --git a/tardis/config_reader.py b/tardis/config_reader.py index 6a5c6824901..aaa65e553a0 100644 --- a/tardis/config_reader.py +++ b/tardis/config_reader.py @@ -39,6 +39,56 @@ def calculate_density_after_time(densities, time_0, time_explosion): return densities * (time_explosion / time_0) ** -3 +def calculate_exponential_densities(velocities, velocity_0, rho_0, exponent): + """ + + This function computes a descret exponential density profile. + :math:`\\rho = \\rho_0 \\times \\left( \\frac{v_0}{v} \\right)^n` + + Parameters + ---------- + + velocities : Array like list + velocities in km/s + + velocity_0 : ~float + Velocity at the inner boundary + + + rho_0 : ~float + density at velocity_0 + + exponent : ~float + exponent used in the powerlaw + + Returns + ------- + + Array like density structure + + """ + densities = rho_0 * (velocity_0 / velocities) ** exponent + return densities[1:] + +def calculate_velocities( v_inner=None, v_outer=None, v_sampling='linear'): + """ + calculates the velocity shells. At the moment only linear v sampling is supported. + + :param velocities: + :param v_inner: + :param v_outer: + :param v_sampling: + + + """ + + if v_sampling == 'linear': + velocities = np.linspace( + v_inner, v_outer, self.no_of_shells + 1) + else: + raise ValueError('Currently only v_sampling = linear is possible') + return velocities + def parse_density_file_section(density_file_dict, time_explosion): density_file_parser = {} @@ -159,9 +209,35 @@ def parse_branch85(density_dict, no_of_shells, v_inner, v_outer, time_explosion) densities = calculate_density_after_time(densities, time_0, time_explosion) + + return densities density_parser['branch85_w7'] = parse_branch85 + def parse_exponential(density_dict,no_of_shells_v_inner, v_outer, time_explosion): + time_0 = density_dict.pop('time_0',19.9999584) + if isinstance(time_0, basestring): + time_0 = parse2quantity(time_0).to('s').value + else: + logger.debug('time_0 not supplied for density branch85 - using sensible default %g', time_0) + try: + rho_0 = density_dict.pop('rho') + except KeyError: + rho_0 = 1e-2 + logger.warning('rho_o was not given in the config! Using %g', rho_0) + try: + exponent = density_dict.pop('exponent') + except: + exponent = 2 + logger.warning('exponent was not given in the config file! Using %f', exponent) + + + velocities = calculate_velocities(v_inner=v_inner, v_outer=v_outer) + densities = calculate_exponential_densities(velocities, v_inner, rho_0, exponent) + + return densities + density_parser['exponential'] = parse_exponential + try: parser = density_parser[density_dict['type']] except KeyError: @@ -236,36 +312,6 @@ def calculate_w7_branch85_densities(velocities, time_explosion, time_0=19.999958 return densities[1:] -def calculate_exponential_densities(velocities, velocity_0, rho_0, exponent): - """ - - This function computes a descret exponential density profile. - :math:`\\rho = \\rho_0 \\times \\left( \\frac{v_0}{v} \\right)^n` - - Parameters - ---------- - - velocities : Array like list - velocities in km/s - - velocity_0 : ~float - Velocity at the inner boundary - - - rho_0 : ~float - density at velocity_0 - - exponent : ~float - exponent used in the powerlaw - - Returns - ------- - - Array like density structure - - """ - densities = rho_0 * (velocity_0 / velocities) ** exponent - return densities[1:] def read_w7_densities(fname=None): @@ -384,7 +430,6 @@ def from_yaml(cls, fname, args=None): if structure_dict != {}: logger.warn('Accepted file for structure (density/velocity) structure ignoring all other arguments: \n%s\n', pprint.pformat(structure_dict, indent=4)) - else: #requiring all keys: no_of_shells, velocity, density if not all([item in structure_dict.keys() for item in ('no_of_shells', 'velocity', 'density')]): From 14489debedc83140187f055e8ee72d65d0a9102c Mon Sep 17 00:00:00 2001 From: "Michael.Klauser" Date: Thu, 7 Mar 2013 20:44:58 +0100 Subject: [PATCH 3/4] the exponential density profile is fixed --- tardis/config_reader.py | 23 ++--------------------- 1 file changed, 2 insertions(+), 21 deletions(-) diff --git a/tardis/config_reader.py b/tardis/config_reader.py index aaa65e553a0..4f3c7aa8d87 100644 --- a/tardis/config_reader.py +++ b/tardis/config_reader.py @@ -70,24 +70,6 @@ def calculate_exponential_densities(velocities, velocity_0, rho_0, exponent): densities = rho_0 * (velocity_0 / velocities) ** exponent return densities[1:] -def calculate_velocities( v_inner=None, v_outer=None, v_sampling='linear'): - """ - calculates the velocity shells. At the moment only linear v sampling is supported. - - :param velocities: - :param v_inner: - :param v_outer: - :param v_sampling: - - - """ - - if v_sampling == 'linear': - velocities = np.linspace( - v_inner, v_outer, self.no_of_shells + 1) - else: - raise ValueError('Currently only v_sampling = linear is possible') - return velocities def parse_density_file_section(density_file_dict, time_explosion): @@ -227,12 +209,11 @@ def parse_exponential(density_dict,no_of_shells_v_inner, v_outer, time_explosion logger.warning('rho_o was not given in the config! Using %g', rho_0) try: exponent = density_dict.pop('exponent') - except: + except KeyError: exponent = 2 logger.warning('exponent was not given in the config file! Using %f', exponent) - - velocities = calculate_velocities(v_inner=v_inner, v_outer=v_outer) + velocities = 0.5 * (v_inner + v_outer) densities = calculate_exponential_densities(velocities, v_inner, rho_0, exponent) return densities From febb56af2ce9475200afbae470c745c650b56eb7 Mon Sep 17 00:00:00 2001 From: "Michael.Klauser" Date: Thu, 7 Mar 2013 20:49:28 +0100 Subject: [PATCH 4/4] the exponential density profile is fixed --- tardis/config_reader.py | 22 +--------------------- 1 file changed, 1 insertion(+), 21 deletions(-) diff --git a/tardis/config_reader.py b/tardis/config_reader.py index 4f3c7aa8d87..e1ec882de98 100644 --- a/tardis/config_reader.py +++ b/tardis/config_reader.py @@ -214,7 +214,7 @@ def parse_exponential(density_dict,no_of_shells_v_inner, v_outer, time_explosion logger.warning('exponent was not given in the config file! Using %f', exponent) velocities = 0.5 * (v_inner + v_outer) - densities = calculate_exponential_densities(velocities, v_inner, rho_0, exponent) + densities = calculate_exponential_densities(velocities, v_inner[0], rho_0, exponent) return densities density_parser['exponential'] = parse_exponential @@ -574,26 +574,6 @@ def from_yaml(cls, fname, args=None): return cls(config_dict) - def set_velocities(self, velocities=None, v_inner=None, v_outer=None, v_sampling='linear'): - """ - Setting the velocities - - :param velocities: - :param v_inner: - :param v_outer: - :param v_sampling: - - - """ - if self.no_of_shells is None: - raise ValueError('Can not set abundances before number of shells have been set') - - if v_sampling == 'linear': - self.velocities = np.linspace( - v_inner, v_outer, self.no_of_shells + 1) - else: - raise ValueError('Currently only v_sampling = linear is possible') - def set_abundances(self, abundances): """ Setting the abundances