diff --git a/SPiCE.py b/SPiCE.py index 6e4beeb..0195231 100755 --- a/SPiCE.py +++ b/SPiCE.py @@ -21,7 +21,7 @@ import processes -__version__ = "0.0.1-alpha" +__version__ = "0.0.2-alpha" def get_class(base_module, class_name): @@ -41,8 +41,11 @@ def read_config_file(self, name): self.context = settings.validate(user_settings) def __init__(self, config_file): + #The code crash with 'FloatingPointError' for the cases: "x/0", overflow, underflow and nan ("0/0" or "inf/inf") + np.seterr(divide='raise',over='raise',under='raise',invalid='raise') if config_file is None: - config_file = 'model.yml.example' + config_file = 'gas_example.yml' + self.read_config_file(config_file) self.phases = {} @@ -72,8 +75,11 @@ def run(self): timesteps_Gyr.append(phase.get_timestep_Gyr()) timestep_Gyr = np.nanmax([np.nanmin(timesteps_Gyr), self.integrator['minimum_timestep_Gyr']]) + #total_mass_sum= 0 for phase in self.phases.values(): + # total_mass_sum += phase.current_mass_Msun() phase.update_mass(timestep_Gyr) + #print(total_mass_sum) self.time_Gyr.append(current_time_Gyr + timestep_Gyr) def update_derivatives(self, term): @@ -81,8 +87,6 @@ def update_derivatives(self, term): raise(-1) - - if __name__ == "__main__": parser = argparse.ArgumentParser( prog="SPiCE", diff --git a/gas_example.yml b/gas_example.yml new file mode 100644 index 0000000..736c839 --- /dev/null +++ b/gas_example.yml @@ -0,0 +1,67 @@ +imf: kroupa +binary_star_rates: 0.40 +dtd_sn: rlp + +integrator: + relative_accuracy: 1.0e-6 + initial_time_Gyr: 0. + final_time_Gyr: 13.7 + minimum_timestep_Gyr: 1.0e-4 #Should not be higher with the current integrator we have + +phases: + gas_HII: + type: ionized_hydrogen + params: + initial_mass_Msun: 100. + initial_temperature_K: 1.0e4 + initial_pressure_cgs: 1.38e-12 #erg/cm^3 + gas_HI: + type: neutral_hydrogen + params: + initial_mass_Msun: 0.0 + initial_temperature_K: 100. + initial_pressure_cgs: 1.38e-12 #erg/cm^3 + gas_H2: + type: molecular_hydrogen + params: + initial_mass_Msun: 0.0 + initial_temperature_K: 10. + initial_pressure_cgs: 1.38e-12 #erg/cm^3 + stars: + type: star + params: + initial_mass_Msun: 0.0 + +processes: + star_formation: + type: timescale_constant + params: + input_phase: gas_H2 + output_phase: stars + timescale_Gyr: 3.478 + recombination: + type: H_recombination + params: + input_phase: gas_HII + output_phase: gas_HI + molecule_formation: + type: H2_formation + params: + input_phase: gas_HI + output_phase: gas_H2 + agent_phase: gas_HII + metallicity: 0.02 #Placeholder! one should put something like 'agent_phase2: dust' or something for better models + photoionization: + type: basic_photoionization + params: + input_phase: gas_HI + output_phase: gas_HII + agent_phase: stars + efficiency: 955.29 + photodissociation: + type: basic_photoionization + params: + input_phase: gas_H2 + output_phase: gas_HI + agent_phase: stars + efficiency: 380.93 \ No newline at end of file diff --git a/phases/GasGenerics.py b/phases/GasGenerics.py new file mode 100644 index 0000000..32bd441 --- /dev/null +++ b/phases/GasGenerics.py @@ -0,0 +1,358 @@ +''' +Ideal gas classes + +Mario Romero May 2019 +''' + +import numpy as np +import astropy.units as u +import astropy.constants as cte +from . import basic + +''' +PARENT CLASSES +classes defined here should not be defined as objects in the main code. +''' +class Gas(basic.Phase): + + #--------------------- + #DEFAULT SETTINGS + #--------------------- + def __init__(self, model, params): + #What everyone does + self.model = model + self.params = {**self.default_settings(), **params} + self.mass_history_Msun = [float(self.params['initial_mass_Msun'])] + self.temperature_history_K = [float(self.params['initial_temperature_K'])] + self.pressure_history_cgs = [float(self.params['initial_pressure_cgs'])] + + self._constants() + self._set() + self._state() + + def reset_timestep(self): + self.dm_dt_Msun_Gyr = 0. + self.dP_dt_cgs_Gyr = 0. + + def current_temperature_K(self): + return self.temperature_history_K[-1] + def current_pressure_cgs(self): + return self.pressure_history_cgs[-1] + + def _constants(self): + #Useful conversions and constants (this will save us LOTS of parentesis, trust me!) + self._solMass_to_g = ((1*u.solMass).to(u.g)).value + self._kB = ((cte.k_B).to(u.erg/u.K)).value #Boltzmann constant in erg/K + self._hP = ((cte.h).to(u.erg*u.s)).value #Planck constant in erg*s + self._hbar = self._hP/(2.*np.pi) + self._c = ((cte.c).to(u.cm/u.s)).value #Speed of light + self._G = ((cte.G).to(u.cm*u.cm*u.cm/(u.g*u.s*u.s))).value #Gravitational constant + + def _set(self): + #All internal variables will be in cgs + + self._temperature = self.temperature_history_K[-1] + self._pressure = self.pressure_history_cgs[-1] + self._gas_mass = self.mass_history_Msun[-1]*self._solMass_to_g + + def update_mass(self, timestep_Gyr): + #Euler integrator is becoming a bottleneck, it should be replaced with a more efficient one. + #Actually, this should be called 'update_all', but is called this way to keep compatibility + #------------------------------ + + #First, I deal with Elmegreen89 Formula to get the pressure (or the box area if first iteration) + new_pressure = self.current_pressure_cgs() + try: + self._box_area_cm2 #It will go to the except block for the first iteration + #self._compute_pressure_derivative() + new_pressure = self._compute_pressure_directly() + except AttributeError: + self._set_box_area() + #Update mass + self.mass_history_Msun.append(self.current_mass_Msun() + self.dm_dt_Msun_Gyr*timestep_Gyr) + + #Update pressure + #self.pressure_history_cgs.append(self.current_pressure_cgs() + self.dP_dt_cgs_Gyr*timestep_Gyr) + self.pressure_history_cgs.append(new_pressure) + #print(self.current_pressure_cgs()/self._kB,self.current_mass_Msun(),type(self)) + + #If you don't want to store the pressure history because, for example, it does not change over time + # you only have to comment the relevant line of these two. + self.temperature_history_K.append(self.current_temperature_K()) + + self._set() + self._state() + + def default_settings(self): + return{ + 'initial_mass_Msun':0.0, + 'initial_temperature_K': 0.0, + 'initial_pressure_cgs':0.0 + } + + #--------------------- + #COMMON ATRIBUTES/METHODS TO ALL DERIVED CLASSES + #--------------------- + def _compute_variables_indepedent_of_state(self): + + try: + self._number_density = self._number_particles / self._volume + self._thermal_energy_density = self._thermal_energy / self._volume + self._gamma = 1.+(self._pressure/self._thermal_energy_density) + self._chemical_potential = self._kB*self._temperature*np.log(self._fugacity) + except (FloatingPointError,ZeroDivisionError): + self._number_density = 0.0 + self._thermal_energy_density = 0.0 + self._gamma = np.Infinity + self._chemical_potential = -np.Infinity + self._mass_density = self._number_density * self._particle_mass_g + #print(self._thermal_energy_density , self._thermal_energy , self._volume) + + def _set_box_area(self): + #Here we use the Elmegreen89 (approximate) formula for pressure: + #P_gas = (pi/2) * G * (S_gas + S_Total) + #Where S is the superficial density = Mass/Area. We don't have the area, so we compute it here as a 'constant' + #Here we get the Area + + #Get total mass, and total GAS mass + total_mass = 0.0 + total_gas_mass = 0.0 + total_gas_pressure = 0.0 + for phase in self.model.phases.values(): + total_mass += phase.mass_history_Msun[0] #We need only the first value, and it may be updated, so using 'current_mass' is a bad idea + if( issubclass(type(phase),Gas) ): + total_gas_mass += phase.mass_history_Msun[0] + total_gas_pressure += phase.pressure_history_cgs[0] + #Convert to cgs + total_mass *= self._solMass_to_g + total_gas_mass *= self._solMass_to_g + #Get the result + self._box_area_cm2 = np.sqrt( np.pi*self._G*total_gas_mass*total_mass / (2.*total_gas_pressure) ) + + def _compute_pressure_directly(self): + #Here we use the Elmegreen89 (approximate) formula for pressure: + #P_i = (pi/2) * G * (S_i + S_Total) + #Where S is the superficial density = Mass/Area. We don't have the area, so we compute it here as a 'constant' + #Here we compute P_i using the formula. + + #Here we need the total mass again. + #Total mass of a previous phase may be updated before in the main routine. We have to cover this! + total_mass = 0.0 + thisPhase_historylength = len(self.mass_history_Msun) + for phase in self.model.phases.values(): + currentPhase_historylength = len(phase.mass_history_Msun) + if( thisPhase_historylength == currentPhase_historylength ): + #currentPhase has not been updated yet, use current value + total_mass += phase.mass_history_Msun[-1] + elif( thisPhase_historylength == currentPhase_historylength-1 ): + #currentPhase has been updated in this timestep before this call + total_mass += phase.mass_history_Msun[-2] + else: + print("This should not happen!") + raise(-1) + #Make the conversion Msun->g + total_mass = total_mass * self._solMass_to_g + current_mass = self.current_mass_Msun() * self._solMass_to_g + #Get the new pressure + return current_mass*total_mass*( np.pi*self._G / (2.*self._box_area_cm2*self._box_area_cm2) ) + + + ''' + def _compute_pressure_derivative(self): + #Here we use the Elmegreen89 (approximate) formula for pressure: + #P_i = (pi/2) * G * (S_i + S_Total) + #Where S is the superficial density = Mass/Area. We don't have the area, so we compute it here as a 'constant' + #Here we compute the derivative of P_i + + #print( (np.sqrt(self._box_area_cm2/np.pi)*u.cm).to(u.pc) ) #Need to check if we have the correct numbers! + + #Here we need the sum of all derivatives of dm_dt, and again the total mass. + #Total mass of a previous phase may be updated before in the main routine. We have to cover this! + total_dm_dt_Msun_Gyr = 0.0 + total_mass = 0.0 + thisPhase_historylength = len(self.mass_history_Msun) + for phase in self.model.phases.values(): + currentPhase_historylength = len(phase.mass_history_Msun) + if( thisPhase_historylength == currentPhase_historylength ): + #currentPhase has not been updated yet, use current value + total_mass += phase.mass_history_Msun[-1] + elif( thisPhase_historylength == currentPhase_historylength-1 ): + #currentPhase has been updated in this timestep before this call + total_mass += phase.mass_history_Msun[-2] + else: + print("This should not happen!") + raise(-1) + #Sum derivatives + total_dm_dt_Msun_Gyr += phase.dm_dt_Msun_Gyr + #Make the conversion Msun->g + total_dm_dt_g_Gyr = total_dm_dt_Msun_Gyr * self._solMass_to_g + curr_dm_dt_g_Gyr = self.dm_dt_Msun_Gyr * self._solMass_to_g + total_mass = total_mass * self._solMass_to_g + current_mass = self.current_mass_Msun() * self._solMass_to_g + #Get the pressure derivative + self.dP_dt_cgs_Gyr = curr_dm_dt_g_Gyr*total_mass + current_mass*total_dm_dt_g_Gyr + self.dP_dt_cgs_Gyr *= np.pi*self._G / (2.*self._box_area_cm2*self._box_area_cm2) + ''' + #--------------------- + #OUTPUTS + #--------------------- + #def mass(...): comes from basic phase. So I do not redefine here + def pressure(self,units=u.erg/(u.cm*u.cm*u.cm)): + return ((self._pressure)*units).value + def temperature(self,units=u.K): + return ((self._temperature)*units).value + def number_particles(self,units=u.m/u.m): + return ((self._number_particles)*units).value + def thermal_energy(self,units=u.erg): + return ((self._thermal_energy)*units).value + def fugacity(self,units=u.m/u.m): + return ((self._fugacity)*units).value + def volume(self,units=(u.cm*u.cm*u.cm)): + return ((self._volume)*units).value + def chemical_potential(self,units=u.erg): + return (self._chemical_potential).to(units) + def number_density(self,units=1./(u.cm*u.cm*u.cm)): + return ((self._number_density)*units).value + def mass_density(self,units=u.g/(u.cm*u.cm*u.cm)): + return ((self._mass_density)*units).value + def thermal_energy_density(self,units=u.erg/(u.cm*u.cm*u.cm)): + return ((self._thermal_energy_density(self))*(units)).value + def adiabatic_constant(self,units=u.m/u.m): + return ((self._gamma)*units).value + + #--------------------- + #POLYMORPHIC METHODS + #--------------------- + #I.e.: These are the ONLY functions you have to redefine in each subclass + def _state(self,P=-1,T=-1,M=-1,N=-1):#inputs: pressure, temperature, number of particles, thermal energy, mass + #print("Warning: Using a parent class. Should not happen!") + raise NameError("Called by parent class. Should not happen!") + + #--------------------- + #DEBUGGING + #--------------------- + def debug(self): + print([self._gas_mass,self._particle_mass_g]) + +''' +CHILDS OF GAS +These classes must have 'state()' well defined and not polymorphic +''' +class Monoatomic(Gas): + + #--------------------- + #STATE METHOD + #--------------------- + def _state(self):#inputs: pressure, temperature, number of particles, thermal energy, mass + #Here we compute the other thermodynamic variables + + self._number_particles = (self._gas_mass)/(self._particle_mass_g) + self._thermal_energy = 1.5*self._number_particles * self._kB * self._temperature + try: + self._volume = (2.*self._thermal_energy) / (3.*self._pressure) + except (FloatingPointError,ZeroDivisionError): + self._volume = 0.0 + + #Fugacity = exp(chemical_potential / kT) . It is another way to describe the chemical potential + self._fugacity = self._pressure * (2.*np.pi*self._particle_mass_g / (self._hP*self._hP) )**(-1.5) * (self._kB*self._temperature)**(-2.5) + #Compute the other variables + self._compute_variables_indepedent_of_state() + +class Diatomic(Gas): + + #--------------------- + #AUXILIARY METHODS + #--------------------- + def _molecule_mass(self): + return self._atom_mass_g[0]+self._atom_mass_g[1] #ESTOY AQUÍ + + def _moment_of_inertia(self): + #Get reduced mass + m_r = self._atom_mass_g[0]*self._atom_mass_g[1] / (self._atom_mass_g[0]+self._atom_mass_g[1]) + return m_r*self._atom_distance_cm*self._atom_distance_cm + + def _angular_frequency(self): + return self._wavenumber_cm_minus1 * self._c + + #--------------------- + #PARTITION FUNCTION METHODS + #--------------------- + def _Z_rot(self): + #We compute Z for only rotations + + #Compute a parameter to check if quantum effects are important + theta = 0.5 * self._hbar*self._hbar / (self._kB * self._temperature * self._moment_of_inertia()) + #Now do stuff + Z = None + dZ_db = None #Partial derivative of Z with respect to beta = 1./kT (units: erg) + if theta <= 1: #All rotational modes are enabled (=classical) + Z = 1./theta + dZ_db = - self._kB * self._temperature / theta + else: #All rotational modes are disabled + Z = 1.0 + dZ_db = 0.0 + + return [Z , dZ_db] + ''' + Note that this is a simplification. + The exact method has Z = sum (2l+1)*exp(-theta*l(l+1)) ; for l=0 to infinity. + Be aware that it's computationally expensive. + ''' + + def _Z_vib(self): + #We compute Z for only molecule vibrations + + #no_units = u.cm/u.cm + #Compute a parameter to check if quantum effects are important + xi = self._hbar * self._angular_frequency() / (self._kB * self._temperature) + #Now do stuff + Z = None + dZ_db = None #Partial derivative of Z with respect to beta = 1./kT (units: erg) + if xi > 1: #All vibrational modes are disabled + Z = 1.0 + dZ_db = 0.0 + else: + Z = 1./xi + dZ_db = - self._kB * self._temperature / xi + + return [Z , dZ_db] + ''' + Again, this is a simplification. + The exact method has Z = sum exp(-n*xi) = 1./(1-exp(-xi)) . + Also note that the zero of energies is set as all vibrations at the fundamental state of the harmonic oscillator. + ''' + + #--------------------- + #STATE METHOD + #--------------------- + def _state(self):#inputs: pressure, temperature, number of particles, mass + + #uV = (u.cm*u.cm*u.cm) #volume units + #no_units = u.m/u.m #dimensionless units + + #Because you give each atom mass individually, you need to construct the particle mass first. + self._particle_mass_g = self._molecule_mass() + #Get the total number of particles + self._number_particles = (self._gas_mass)/(self._particle_mass_g) + #Let's define an useful variable + NkT = self._number_particles * (self._kB) * self._temperature + + #Partition functions (they are arrays, [0] is the Z itself, and [1] the derivative with respect to beta) + Zrot = self._Z_rot() + Zvib = self._Z_vib() + + #Now we compute the remaining variables + try: + self._volume = NkT / self._pressure + except (FloatingPointError,ZeroDivisionError): + self._volume = 0.0 + #Fugacity = exp(chemical_potential / kT) . It is another way to describe the chemical potential + self._fugacity = self._pressure * (2.*np.pi*self._particle_mass_g / (self._hP*self._hP) )**(-1.5) * (self._kB*self._temperature)**(-2.5) / (Zrot[0]*Zvib[0]) + #Thermal energy + #const_factor = 1.5 - (Zrot[1]/( Zrot[0] * self._kB * self._temperature )) - (Zvib[1]/( Zvib[0] * self._kB * self._temperature )) + const_factor = 1.5 #Translational part + const_factor -= (Zrot[1]/( Zrot[0] * self._kB * self._temperature )) #Rotational part + const_factor -= (Zvib[1]/( Zvib[0] * self._kB * self._temperature )) #Vibrational part + self._thermal_energy = NkT * const_factor + + self._compute_variables_indepedent_of_state() diff --git a/phases/GasSpecies.py b/phases/GasSpecies.py new file mode 100644 index 0000000..c88a74f --- /dev/null +++ b/phases/GasSpecies.py @@ -0,0 +1,41 @@ +''' +Ideal gas classes + +Mario Romero May 2019 +''' + +#import numpy as np +import astropy.units as u +import astropy.constants as cte +from . import GasGenerics as gas + +''' +SOME CONSTANTS +''' + +atomic_mass_g = ((cte.u).to(u.g)).value +proton_mass_g = ((cte.m_p).to(u.g)).value + +''' +SPECIFIC GASES +This script contains all realistic gas species (e.g. Hydrogen, Helium, etc). They are subclasses of classes from 'GasGenerics.py' +''' + +#-------------- +#MONOATOMIC SPECIES +#Required data = particle mass +#-------------- +class Neutral_Hydrogen(gas.Monoatomic): + _particle_mass_g = 1.008*atomic_mass_g + +class Ionized_Hydrogen(gas.Monoatomic): + _particle_mass_g = proton_mass_g #Proton mass + +#-------------- +#DIATOMIC SPECIES +#Needed data: particle masses (array), distance between particles and their fundamental vibration mode. +#-------------- +class Molecular_Hydrogen(gas.Diatomic): + _atom_mass_g = [1.008*atomic_mass_g , 1.008*atomic_mass_g ] #Two hydrogen atoms + _atom_distance_cm = 74.14*1e-10 #Taken from wikipedia + _wavenumber_cm_minus1 = 4342.0 #wavenumber (k) of fundamental vibration. Taken from here: https://www.chem.purdue.edu/gchelp/vibs/h2.html // diff --git a/phases/StarGenerics.py b/phases/StarGenerics.py new file mode 100644 index 0000000..cf25106 --- /dev/null +++ b/phases/StarGenerics.py @@ -0,0 +1,26 @@ +''' +Star classes + +Mario Romero May 2019 +''' + +import numpy as np +import astropy.units as u +import astropy.constants as cte +from . import basic + +''' +PLACEHOLDER CLASS +Defined here to make photo-ionization and photo-dissociation processes work +''' +class Star(basic.Phase): + + #--------------------- + #DEFAULT SETTINGS + #--------------------- + def __init__(self, model, params): + #What everyone does + self.model = model + self.params = {**self.default_settings(), **params} + self.mass_history_Msun = [float(self.params['initial_mass_Msun'])] + self.SFR_history_Msun_per_Gyr = [] #It updates when a star formation process is called \ No newline at end of file diff --git a/phases/__init__.py b/phases/__init__.py index da2cb55..ad2f126 100644 --- a/phases/__init__.py +++ b/phases/__init__.py @@ -1,8 +1,16 @@ from . import basic +# from . import GasGenerics +from . import StarGenerics +from . import GasSpecies from . import example + registry = { - "basic": basic.Phase, - "multiphase": basic.MultiphaseMedium, - "dust": example.dust.Dust, - } + "basic": basic.Phase, + "multiphase": basic.MultiphaseMedium, + "star": StarGenerics.Star, #PLACEHOLDER! + "ionized_hydrogen": GasSpecies.Ionized_Hydrogen, + "neutral_hydrogen": GasSpecies.Neutral_Hydrogen, + "molecular_hydrogen": GasSpecies.Molecular_Hydrogen, + "dust": example.dust.Dust, +} diff --git a/phases/basic.py b/phases/basic.py index e35ac3c..225f96c 100644 --- a/phases/basic.py +++ b/phases/basic.py @@ -28,7 +28,7 @@ def current_mass_Msun(self): return self.mass_history_Msun[-1] def mass(self, t): - return np.interp(t, self.model.time, self.mass_history[-1]) + return np.interp(t, self.model.time, self.mass_history_Msun[-1]) def reset_timestep(self): self.dm_dt_Msun_Gyr = 0. @@ -37,10 +37,15 @@ def update_derivatives(self, term): self.dm_dt_Msun_Gyr += term def get_timestep_Gyr(self): - return np.abs((self.model.integrator['relative_accuracy'] + try: + return np.abs((self.model.integrator['relative_accuracy'] * self.current_mass_Msun() / self.dm_dt_Msun_Gyr)) + except (FloatingPointError,ZeroDivisionError): + #How should the previous expression be solved if gives 0/0? + return np.Infinity def update_mass(self, timestep_Gyr): + #Euler integrator is becoming a bottleneck, it should be replaced with a more efficient one. self.mass_history_Msun.append(self.current_mass_Msun() + self.dm_dt_Msun_Gyr*timestep_Gyr) diff --git a/processes/Ionization.py b/processes/Ionization.py new file mode 100644 index 0000000..fa76177 --- /dev/null +++ b/processes/Ionization.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Processes related with ionization proccesses + +Mario Romero July 2019 +""" + +import numpy as np +import astropy.units as u +from . import basic + +''' +Hydrogen Photoionization +HI -> HII process +''' +class Basic_Photoionization(basic.Process): + + #--------------------- + #INIT + #--------------------- + def __init__(self, model, params): + self.input = model.phases[params['input_phase']] + self.output = model.phases[params['output_phase']] + self.agent = model.phases[params['agent_phase']] #Stars + self._eff = float(params['efficiency']) #=Mass of photoionized gas per Mass of stars + + def compute_derivatives(self): + #flux = self.input.current_mass_Msun()/self.tau_Gyr + try: + flux = self._eff * self.agent.SFR_history_Msun_per_Gyr[-1] + except IndexError: + flux = 0.0 + #print(self.agent.SFR_history_Msun_per_Gyr[0]) + self.input.update_derivatives(-flux) + self.output.update_derivatives(flux) + + try: + self.tau_Gyr = self.input.current_mass_Msun() / flux + except (FloatingPointError,ZeroDivisionError): + self.tau_Gyr = np.Infinity + #print(self.model.input) \ No newline at end of file diff --git a/processes/MoleculeFormation.py b/processes/MoleculeFormation.py new file mode 100644 index 0000000..5db9272 --- /dev/null +++ b/processes/MoleculeFormation.py @@ -0,0 +1,63 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Processes related with formation of molecules + +Mario Romero July 2019 +""" + +import numpy as np +import astropy.units as u +from . import basic + +''' +Hydrogen molecule creation +H -> H2 process +''' +class MolecularHydrogen_Formation(basic.Process): + + #--------------------- + #INIT + #--------------------- + def __init__(self, model, params): + self.input = model.phases[params['input_phase']] #Neutral Hydrogen + self.output = model.phases[params['output_phase']] #Molecular Hydrogen + self.agent = model.phases[params['agent_phase']] #Ionized Hydrogen + self._Z = float(params['metallicity']) + + self._constants() + self.tau_Gyr = self._formation_timescale() + + def compute_derivatives(self): + flux = self.input.current_mass_Msun()/self.tau_Gyr + self.input.update_derivatives(-flux) + self.output.update_derivatives(flux) + self.tau_Gyr = self._formation_timescale() + + #--------------------- + #GETTING PHYSICAL PARAMETERS + #--------------------- + def _dust_density(self): + #total_H_density_cm3 = 0.5*self.agent.number_density() + self.input.number_density() + 2.*self.output.number_density() + total_V_cm3 = self.agent.volume() + self.input.volume() + self.output.volume() + total_H_density_cm3 = (0.5*self.agent.volume()/total_V_cm3)*self.agent.number_density() + (self.input.volume()/total_V_cm3)*self.input.number_density()+ (2.*self.output.volume()/total_V_cm3)*self.output.number_density() + #print(total_H_density_cm3) + return self._Z * total_H_density_cm3 + def _formation_timescale(self): + #Taking Ascasibar+(in prep) formula again, ignoring the effective metallicity thing + mean_ov_cm3_per_s = 6.e-17 * ((self.input.temperature()/100.)**(0.5)) #Neutral H or H2 temperature???? + #mean_ov_cm3_Gyr = 0.189*(self.input.temperature()/100.)**(0.5) + + tau_s = None #Not the correct tau unit-wise + try: + tau_s = 0.5/(mean_ov_cm3_per_s * self._dust_density()) + except ZeroDivisionError: + tau_s = np.Infinity + + return tau_s*self._s_to_Gyr #Correct units. + + #--------------------- + #DEFINING USEFUL VARIABLES + #--------------------- + def _constants(self): + self._s_to_Gyr = ((1.*u.s).to(u.Gyr)).value \ No newline at end of file diff --git a/processes/Recombination.py b/processes/Recombination.py new file mode 100644 index 0000000..fce0347 --- /dev/null +++ b/processes/Recombination.py @@ -0,0 +1,54 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Processes related with recombination + +Mario Romero July 2019 +""" + +import numpy as np +import astropy.units as u +from . import basic + +''' +Hydrogen Recombination +HI -> H process +''' +class Hydrogen_Recombination(basic.Process): + + #--------------------- + #INIT + #--------------------- + def __init__(self, model, params): + self.input = model.phases[params['input_phase']] #Ionized Hydrogen + self.output = model.phases[params['output_phase']] #Neutral Hydrogen + + self._constants() + self.tau_Gyr = self._recombination_timescale() + + def compute_derivatives(self): + flux = self.input.current_mass_Msun()/self.tau_Gyr + #print(self.input.current_mass_Msun(),self.tau_Gyr,flux,'\n') + self.input.update_derivatives(-flux) + self.output.update_derivatives(flux) + self.tau_Gyr = self._recombination_timescale() + #print(self.tau_Gyr) + #--------------------- + #GETTING THE TAU + #--------------------- + def _recombination_timescale(self): + #Taking Ascasibar+(in prep) formula + mean_ov_cm3_s = 4.1e-10 * ((self.input.temperature())**(-0.8)) #Case B recombination cross-section (Verner&Ferland 96) + + tau_s = None #Not the correct tau unit-wise + try: + tau_s = 2./(mean_ov_cm3_s * self.input.number_density()) + except ZeroDivisionError: + tau_s = np.Infinity + + return tau_s*self._s_to_Gyr #Correct units. + #--------------------- + #DEFINING USEFUL VARIABLES + #--------------------- + def _constants(self): + self._s_to_Gyr = ((1.*u.s).to(u.Gyr)).value \ No newline at end of file diff --git a/processes/__init__.py b/processes/__init__.py index ab526f2..a34cf86 100644 --- a/processes/__init__.py +++ b/processes/__init__.py @@ -1,5 +1,11 @@ from . import basic +from . import Recombination +from . import MoleculeFormation +from . import Ionization registry = { "timescale_constant": basic.Constant_timescale, + "H_recombination": Recombination.Hydrogen_Recombination, + "H2_formation": MoleculeFormation.MolecularHydrogen_Formation, + "basic_photoionization": Ionization.Basic_Photoionization, } diff --git a/processes/basic.py b/processes/basic.py index cd87e46..c62119b 100644 --- a/processes/basic.py +++ b/processes/basic.py @@ -26,3 +26,5 @@ def compute_derivatives(self): flux = self.input.current_mass_Msun()/self.tau_Gyr self.input.update_derivatives(-flux) self.output.update_derivatives(flux) + #PLACEHOLDER! This will fail if two processes give an SFR as main result + self.output.SFR_history_Msun_per_Gyr.append(flux) diff --git a/processes/star_formation.py b/processes/star_formation.py index 2ba27fd..b9708af 100644 --- a/processes/star_formation.py +++ b/processes/star_formation.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Created on Wed Mar 27 10:30:28 2019 @@ -21,3 +19,4 @@ def compute_derivatives(self): SFR = self.gas.mass()/self.tau self.gas.update_derivatives(-SFR) self.stars.update_derivatives(SFR) + diff --git a/settings.py b/settings.py index 146f44d..753c194 100644 --- a/settings.py +++ b/settings.py @@ -22,6 +22,7 @@ valid_values = { 'imf': ['salpeter', 'starburst', 'chabrier', 'ferrini', 'kroupa', 'miller_scalo', 'maschberger'], 'dtd_sn': ['rlp', 'mdvp'], # rlp = Ruiz-Lapuente, mdvp = Mannucci, Della Valle, Panagia (2006) + #'phases': ['gas','ionized_hydrogen', 'neutral_hydrogen', 'molecular_hydrogen'], } def validate(params):