import warnings
from .funcsplots import *
[docs]def dic_usual_gadget2_units():
"""Returns a dictionary of usual units used by Gadget-2.
Returns
-------
dict
Usual units used by Gadget-2.
"""
dicunits = {'length_in_cgs': BaseUnitsDefs.CM_PER_PC * 1e3, 'vel_in_cgs': 1e5,
'mass_in_cgs': BaseUnitsDefs.CGS_SOLAR_MASS * 1e10,
'str_length_in_cgs': '3.085678e21', 'str_vel_in_cgs': '1e5', 'str_mass_in_cgs': '1.989e43'}
# G in this set of units
dicunits['gravity_constant_internal'] = BaseUnitsDefs.CGS_GRAVITY / (
dicunits['length_in_cgs'] * dicunits['vel_in_cgs'] ** 2 / (dicunits['mass_in_cgs']))
return dicunits
[docs]def dic_unity_gravity_const_units():
"""Returns a dictionary of set of units so that the gravtational constant is 1.
Returns
-------
dict
Set of units so that the gravtational constant is 1.
"""
dicunits = {'length_in_cgs': BaseUnitsDefs.CM_PER_PC * 1e3, 'vel_in_cgs': 1e7,
'mass_in_cgs': BaseUnitsDefs.CGS_SOLAR_MASS * 2.32519716e9, 'str_length_in_cgs': '3.085678e21',
'str_vel_in_cgs': '1e7', 'str_mass_in_cgs': '4.62481715e+42'}
# G in this set of units
dicunits['gravity_constant_internal'] = BaseUnitsDefs.CGS_GRAVITY / (
dicunits['length_in_cgs'] * dicunits['vel_in_cgs'] ** 2 / (dicunits['mass_in_cgs']))
return dicunits
[docs]class BaseUnitsDefs:
"""Class with basic definitions (values of solar mass, etc in CGS and conversion class methods of CGS to usual astro
units. Units variables are taken from Gadget-2.
"""
def __init__(self):
pass
CGS_GRAVITY = 6.672e-8
CGS_SOLAR_MASS = 1.989e33
CGS_SOLAR_LUM = 3.826e33
CGS_RAD_CONST = 7.565e-15
CGS_AVOGADRO = 6.0222e23
CGS_BOLTZMANN = 1.3806e-16
CGS_GAS_CONST = 8.31425e7
CGS_C = 2.9979e10
CGS_PLANCK = 6.6262e-27
CGS_PROTONMASS = 1.6726e-24
CGS_ELECTRONMASS = 9.10953e-28
CGS_THOMPSON = 6.65245e-25
CGS_ELECTRONCHARGE = 4.8032e-10
CGS_HUBBLE = 3.2407789e-18 # in h/sec
CGS_ELECTRONVOLT = 1.602177e-12 # in erg
CGS_EVONC2MASS = 1.782692e-33 # in g
SEC_PER_MEGAYEAR = 3.155e13
SEC_PER_YEAR = 3.155e7
CM_PER_MP = 3.085678e24
CM_PER_PC = 3.085678e18
HYDROGEN_MASS_FRACTION = 0.76
GAMMA_MONOATOMIC_GAS = 5. / 3.
[docs] @classmethod
def cgs_length_to_kpc(cls, unit_length_in_cgs):
"""Conversion factor from position unit in cgs to kpc.
Parameters
----------
unit_length_in_cgs : float
Length unit in cm.
Returns
-------
float
Conversion factor.
"""
facconv = unit_length_in_cgs / (1e3 * cls.CM_PER_PC)
# Conversion factor set to one if it is already very close to one.
if (np.fabs(facconv - 1) < 1e-5):
facconv = 1.
return facconv
[docs] @classmethod
def cgs_vel_to_kmpers(cls, unit_vel_in_cgs):
"""Conversion factor from velocity unit in cgs to km/s.
Parameters
----------
unit_vel_in_cgs : float
Velocity unit in cm/s.
Returns
-------
float
Conversion factor.
"""
return unit_vel_in_cgs / 1e5
[docs] @classmethod
def cgs_mass_to_sunm(cls, unit_mass_in_cgs):
"""Conversion factor from mass unit in cgs to solar mass.
Parameters
----------
unit_mass_in_cgs : float
Mass unit in g.
Returns
-------
float
Conversion factor.
"""
return unit_mass_in_cgs / cls.CGS_SOLAR_MASS
[docs] @classmethod
def cgs_surface_density_to_sunm_per_pc2(cls, unit_mass_in_cgs, unit_length_in_cgs):
"""Conversion factor from surface density unit in cgs to solar mass per pc^2.
Parameters
----------
unit_mass_in_cgs : float
Mass unit in g.
unit_length_in_cgs : float
Length unit in cm.
Returns
-------
float
Conversion factor.
"""
return unit_mass_in_cgs / cls.CGS_SOLAR_MASS * (cls.CM_PER_PC / unit_length_in_cgs) ** 2
[docs] @classmethod
def cgs_volume_density_to_sunm_per_pc3(cls, unit_mass_in_cgs, unit_length_in_cgs):
"""Conversion factor from volume density unit in cgs to solar mass per pc^3.
Parameters
----------
unit_mass_in_cgs : float
Mass unit in g.
unit_length_in_cgs : float
Length unit in cm.
Returns
-------
float
Conversion factor.
"""
return unit_mass_in_cgs / cls.CGS_SOLAR_MASS * (cls.CM_PER_PC / unit_length_in_cgs) ** 3
[docs] @classmethod
def cgs_volume_nb_density_to_per_pc3(cls, unit_length_in_cgs):
"""Conversion factor from volume number density unit in cgs to per pc^3.
Parameters
----------
unit_length_in_cgs : float
Length unit in cm.
Returns
-------
float
Conversion factor.
"""
return (cls.CM_PER_PC / unit_length_in_cgs) ** 3
[docs] @classmethod
def cgs_time_to_gyr(cls, unit_time_in_cgs):
"""Conversion factor from time unit in cgs to Gyr.
Parameters
----------
unit_time_in_cgs : float
Time unit in s.
Returns
-------
float
Conversion factor.
"""
return unit_time_in_cgs / (1e3 * cls.SEC_PER_MEGAYEAR)
[docs] @classmethod
def get_conversion_factor_from_string(cls, arrstr, dicunits, conv_type=None):
"""Conversion to usual practical astrophysical units.
Parameters
----------
arrstr : str
Member of :attr:`~galaximview.snap.Snap.listquantities`.
dicunits : dict or `None`
Dictionary of simulation units. `None` for no conversion whatever conv_type is set to.
conv_type : {`None`, 'astro', 'cgs'}
Type of units conversion.
Returns
-------
float
conversion factor
"""
fac = 1
if dicunits is not None:
if ((arrstr == 'xpos') | (arrstr == 'ypos') | (arrstr == 'zpos') | (arrstr == 'xpos_ini') | (
arrstr == 'ypos_ini') | (arrstr == 'zpos_ini') | (arrstr == 'radcyl') | (arrstr == 'radspher') | (
arrstr == 'hsml')):
if conv_type == 'cgs':
fac = dicunits['length_in_cgs']
elif conv_type == 'astro':
fac = cls.cgs_length_to_kpc(dicunits['length_in_cgs'])
elif ((arrstr == 'xvel') | (arrstr == 'yvel') | (arrstr == 'zvel') | (arrstr == '-zvel')
| (arrstr == 'v_R') | (arrstr == 'v_phi') | (arrstr == 'velradspher') | (arrstr == 'velthetaspher')):
if conv_type == 'cgs':
fac = dicunits['vel_in_cgs']
elif conv_type == 'astro':
fac = cls.cgs_vel_to_kmpers(dicunits['vel_in_cgs'])
elif ((arrstr == 'mass') | (arrstr == 'massinr')):
if conv_type == 'cgs':
fac = dicunits['mass_in_cgs']
elif conv_type == 'astro':
fac = cls.cgs_mass_to_sunm(dicunits['mass_in_cgs'])
elif ((arrstr == 'jang_norm') | (arrstr == 'zjang')):
if conv_type == 'cgs':
fac = dicunits['mass_in_cgs'] * dicunits['length_in_cgs'] * dicunits['vel_in_cgs']
elif conv_type == 'astro':
fac = cls.cgs_mass_to_sunm(dicunits['mass_in_cgs']) * cls.cgs_length_to_kpc(
dicunits['length_in_cgs']) * cls.cgs_vel_to_kmpers(dicunits['vel_in_cgs'])
elif ((arrstr == 'spec_jang_norm') | (arrstr == 'spec_zjang')):
if conv_type == 'cgs':
fac = dicunits['length_in_cgs'] * dicunits['vel_in_cgs']
elif conv_type == 'astro':
fac = cls.cgs_length_to_kpc(dicunits['length_in_cgs']) * cls.cgs_vel_to_kmpers(
dicunits['vel_in_cgs'])
elif ((arrstr == 'pot') | (arrstr == 'uint')):
if conv_type == 'cgs':
fac = dicunits['vel_in_cgs'] ** 2
elif conv_type == 'astro':
fac = cls.cgs_vel_to_kmpers(dicunits['vel_in_cgs']) ** 2
elif (arrstr == 'energy'):
if conv_type == 'cgs':
fac = dicunits['mass_in_cgs'] * dicunits['vel_in_cgs'] ** 2
elif conv_type == 'astro':
fac = cls.cgs_mass_to_sunm(dicunits['mass_in_cgs']) * cls.cgs_vel_to_kmpers(
dicunits['vel_in_cgs']) ** 2
elif ((arrstr == 'rho') | (arrstr == 'voldens')):
if conv_type == 'cgs':
fac = dicunits['mass_in_cgs'] / dicunits['length_in_cgs'] ** 3
elif conv_type == 'astro':
fac = cls.cgs_volume_density_to_sunm_per_pc3(dicunits['mass_in_cgs'], dicunits['length_in_cgs'])
elif (arrstr == 'nH'):
if conv_type == 'cgs':
fac = 1. / dicunits['length_in_cgs'] ** 3
elif conv_type == 'astro':
fac = cls.cgs_volume_nb_density_to_per_pc3(dicunits['length_in_cgs'])
elif (arrstr == 'entropy'):
if conv_type == 'cgs':
fac = (dicunits['mass_in_cgs'] / dicunits['length_in_cgs'] ** 3) ** (
1 - BaseUnitsDefs.GAMMA_MONOATOMIC_GAS) * dicunits['vel_in_cgs'] ** 2
elif conv_type == 'astro':
fac = cls.cgs_volume_density_to_sunm_per_pc3(dicunits['mass_in_cgs'],
dicunits['length_in_cgs']) ** (
1 - BaseUnitsDefs.GAMMA_MONOATOMIC_GAS) * cls.cgs_vel_to_kmpers(
dicunits['vel_in_cgs']) ** 2
elif (arrstr == 'pressure'):
if conv_type == 'cgs':
fac = dicunits['mass_in_cgs'] / dicunits['length_in_cgs'] ** 3 * dicunits['vel_in_cgs'] ** 2
elif conv_type == 'astro':
fac = cls.cgs_volume_density_to_sunm_per_pc3(dicunits['mass_in_cgs'],
dicunits['length_in_cgs']) * cls.cgs_vel_to_kmpers(
dicunits['vel_in_cgs']) ** 2
elif (arrstr == 'surfdens'):
if conv_type == 'cgs':
fac = dicunits['mass_in_cgs'] / dicunits['length_in_cgs'] ** 3
elif conv_type == 'astro':
fac = cls.cgs_surface_density_to_sunm_per_pc2(dicunits['mass_in_cgs'], dicunits['length_in_cgs'])
elif (arrstr == 'time'):
if conv_type == 'cgs':
fac = dicunits['timeincgs']
elif conv_type == 'astro':
fac = cls.cgs_time_to_gyr(dicunits['timeincgs'])
else:
warnings.warn("No units conversion because unknown array type for unit conversion.")
return fac
[docs] @staticmethod
def get_dic_of_labels():
"""Labels corresponding to usual quantities and used in plots.
Returns
-------
dict
Dictionary of labels corresponding to Members of :attr:`~galaximview.snap.Snap.listquantities` .
"""
diclabs = {'xpos': r'x', 'ypos': r'y', 'zpos': r'z', 'xpos_ini': r'x', 'ypos_ini': r'y',
'zpos_ini': r'z', 'xvel': r'v_x', 'yvel': r'v_y', 'zvel': r'v_z', '-zvel': r'v_{\mathrm{los}}',
'nH': r'n_H', 'entropy': r'A', 'pressure': r'P', 'uint': r'u', 'hsml': r'h',
'stelluvflux': r'UV flux', 'rho': r'\rho', 'pot': r'\phi ', 'energy': r'\mathrm{Energy}',
'radcyl': r'R', 'phi': r'\phi',
'radspher': r'r', 'thetaspher': r'\theta ', 'v_R': r'v_R ', 'v_phi': r'v_{\phi}',
'velradspher': r'v_r', 'velthetaspher': r'v_{\theta} ', 'jang_norm': r'J', 'zjang': r'J_z',
'spec_jang_norm': r'j', 'spec_zjang': r'j_z', 'voldens': r'\rho', 'surfdens': r'\Sigma',
'massinr': r'M(<r)', 'zgas': r'Z_g', 'temperature': r'T', 'time': r'\mathrm{Time}'}
for key in diclabs.keys():
diclabs[key] += r'\,'
return diclabs
[docs] @staticmethod
def get_dic_of_labels_units(conv_type=None):
"""Units to put in labels for usual quantities and used in plots.
Parameters
----------
conv_type : {None, 'astro', 'cgs'}
Type of units conversion.
Returns
-------
dictionary
Dictionary of units corresponding to Members of :attr:`~galaximview.snap.Snap.listquantities` .
"""
if conv_type == 'astro':
diclabsu = {'UL': r' [\mathrm{kpc}]', 'UV': r' [\mathrm{km/s}]', 'UMass': r' [\mathrm{M}_{\odot}]',
'USurfDens': r' [\mathrm{M}_{\odot}/\mathrm{pc}^2]',
'UVolDens': r' [\mathrm{M}_{\odot}/\mathrm{pc}^3]', 'UNumVolDens': r' [\mathrm{pc}^{-3}]',
'UEntropy': r'[(\mathrm{M}_{\odot}/\mathrm{pc}^3)^{1-\gamma} \mathrm{km}^2/\mathrm{'
r's}^2]',
'UPressure': r' [\mathrm{M}_{\odot}/\mathrm{pc}^{3} \mathrm{km}^2/\mathrm{s}^2]',
'UMassFlux': r' [\mathrm{M}_{\odot}/\mathrm{Gyr}]', 'UPOT': r' [\mathrm{km}^2/\mathrm{s}^2]',
'UEnergy': r' [\mathrm{M}_{\odot} \, \mathrm{km}^{2} / s^{2}]',
'UAM': r' [\mathrm{M}_{\odot} \, \mathrm{kpc} \, \mathrm{km}/s]',
'UAMspec': r' [\mathrm{kpc} \, \mathrm{km}/s]', 'UTime': r' [\mathrm{Gyr}]'}
elif conv_type == 'cgs':
diclabsu = {'UL': r' [\mathrm{cm}]', 'UV': r' [\mathrm{cm/s}]', 'UMass': r' [g]',
'USurfDens': r' [\mathrm{g}/\mathrm{cm}^2]', 'UVolDens': r' [\mathrm{g}/\mathrm{cm}^3]',
'UNumVolDens': r' [\mathrm{cm}^{-3}]',
'UEntropy': r' [(\mathrm{g}/\mathrm{cm}^3)^{1-\gamma} \mathrm{cm}^2/\mathrm{s}^2]',
'UPressure': r' [\mathrm{g}/\mathrm{cm}^3 \mathrm{cm}^2/\mathrm{s}^2]',
'UMassFlux': r' [\mathrm{g}/\mathrm{s}]', 'UPOT': r' [\mathrm{cm}^2/\mathrm{s}^2]',
'UEnergy': r' [\mathrm{erg}]',
'UAM': r' [\mathrm{g} \, \mathrm{cm} \, \mathrm{cm}/s]',
'UAMspec': r' [\mathrm{cm} \, \mathrm{cm}/s]', 'UTime': r' [s]'}
else:
diclabsu = {'UL': ' ', 'UV': ' ', 'UMass': ' ', 'USurfDens': ' ', 'UVolDens': ' ', 'UNumVolDens': ' ',
'UEntropy': ' ', 'UPressure': ' ', 'UMassFlux': ' ', 'UPOT': ' ', 'UEnergy': ' ', 'UAM': ' ',
'UAMspec': ' ', 'UTemp': r' [K]', 'UTime': ' '}
diclabs = {'xpos': diclabsu['UL'], 'ypos': diclabsu['UL'], 'zpos': diclabsu['UL'], 'xpos_ini': diclabsu['UL'],
'ypos_ini': diclabsu['UL'], 'zpos_ini': diclabsu['UL'], 'radcyl': diclabsu['UL'],
'radspher': diclabsu['UL'], 'hsml': diclabsu['UL'], 'xvel': diclabsu['UV'], 'yvel': diclabsu['UV'],
'zvel': diclabsu['UV'], '-zvel': diclabsu['UV'], 'v_R': diclabsu['UV'], 'v_phi': diclabsu['UV'],
'velradspher': diclabsu['UV'], 'velthetaspher': diclabsu['UV'], 'jang_norm': diclabsu['UAM'],
'zjang': diclabsu['UAM'], 'spec_jang_norm': diclabsu['UAMspec'], 'spec_zjang': diclabsu['UAMspec'],
'pot': diclabsu['UPOT'], 'phi': '', 'thetaspher': '', 'uint': diclabsu['UPOT'],
'energy': diclabsu['UEnergy'],
'rho': diclabsu['UVolDens'], 'voldens': diclabsu['UVolDens'], 'nH': diclabsu['UNumVolDens'],
'entropy': diclabsu['UEntropy'], 'pressure': diclabsu['UPressure'],
'surfdens': diclabsu['USurfDens'], 'massinr': diclabsu['UMass'], 'zgas': '',
'temperature': r' \, [\mathrm{K}]', 'time': diclabsu['UTime']}
return diclabs
[docs] @classmethod
def get_label_from_string(cls, arrstr, func_for_label=identity_function, label_of_dispersion=False,
conv_type=None):
"""Finds label corresponding to a function func of a usual quantity, with proper unit.
Parameters
----------
arrstr : str
Member of :attr:`~galaximview.snap.Snap.listquantities`.
func_for_label : ufunc, optional
function of array: identity_function (no transformation), np.log, np.log10 or np.sqrt.
Defaults to identity_function.
label_of_dispersion : bool, default = ` False`
`True` if the label is for a dispersion. (Default value = False)
conv_type : {None, 'astro', 'cgs'}
Type of units conversion.
Returns
-------
str
Label corresponding to member arrstr of cls.listquantities.
"""
lab = np.zeros((1), dtype='object')
labbeg = np.zeros((1), dtype='object')
labmid = np.zeros((1), dtype='object')
labend = np.zeros((1), dtype='object')
labbegsig = np.zeros((1), dtype='object')
labendsig = np.zeros((1), dtype='object')
diclabs = cls.get_dic_of_labels()
if not (arrstr in diclabs.keys()):
lab[0] = ''
warnings.warn('No label can be found for the input.')
else:
lab[0] = diclabs[arrstr]
if (func_for_label == np.log10):
labbeg[0] = r'\mathrm{log_{10}}('
labmid[0] = r'/ '
labend[0] = r')'
elif (func_for_label == np.log):
labbeg[0] = r'\mathrm{log}('
labmid[0] = r'/ '
labend[0] = r')'
elif (func_for_label == np.sqrt):
labbeg[0] = r'\sqrt{'
labmid[0] = r'/ '
labend[0] = r'}'
else:
labbeg[0] = ''
labmid[0] = ''
labend[0] = ''
if (label_of_dispersion):
labbegsig[0] = r'\sigma_{'
labendsig[0] = r'}'
else:
labbegsig[0] = ''
labendsig[0] = ''
lab = lab[0]
labbeg = labbeg[0]
labmid = labmid[0]
labend = labend[0]
labbegsig = labbegsig[0]
labendsig = labendsig[0]
lab = labbeg + labbegsig + lab + labendsig + labmid
dicu = cls.get_dic_of_labels_units(conv_type=conv_type)
if not (arrstr in dicu.keys()):
labu = ''
warnings.warn('No units can be found for the input.')
else:
labu = dicu[arrstr]
lab = '$' + lab + labu + labend + '$'
return lab
[docs]def used_params(dirname):
"""Returns dictionary from parameter-used-values file of a simulation.
Parameters
----------
dirname : str
Returns
-------
dict
Dictionary of all parameters of parameters file.
"""
usedparamsfile = dirname + "/" + "parameters-usedvalues"
f = open(usedparamsfile)
params = {}
for line in f.readlines():
sepline = line.split(' ', maxsplit=1)
keydic = sepline[0].split()[0]
valdic = sepline[1].split()[0]
params[keydic] = valdic
f.close()
return params
[docs]def energy_stats(dirname):
"""Returns array from loadtxt of energy file of Gadget-2 simulation.
Parameters
----------
dirname : str
Path to directory fo simulation.
Returns
-------
numpy.loadtxt of file
Energyfile stored in array.
"""
energyfile = dirname + "/" + "energy.txt"
fen = np.loadtxt(energyfile)
return fen
[docs]class SimCaracs(BaseUnitsDefs):
"""Gets units and energy diagnostics of Gadget-2 simulation.
Parameters
----------
dirname: str
Path to simulation directory.
Attributes
----------
dicunits : dict
Dictionary with simulations unis, read for parameters-used-values file if Gadget simulation.
"""
def __init__(self, dirname):
BaseUnitsDefs.__init__(self)
self.dirname = dirname
self.usedparams = used_params(dirname)
self.dicunits = {'length_in_cgs': float(self.usedparams['UnitLength_in_cm']),
'vel_in_cgs': float(self.usedparams['UnitVelocity_in_cm_per_s']),
'mass_in_cgs': float(self.usedparams['UnitMass_in_g']),
'timeincgs': float(self.usedparams['UnitLength_in_cm']) / float(
self.usedparams['UnitVelocity_in_cm_per_s'])}
self.dicunits['length_in_kpc'] = self.dicunits['length_in_cgs'] / (1e3 * self.CM_PER_PC)
self.dicunits['vel_in_km_per_s'] = self.dicunits['vel_in_cgs'] * 1e-5
self.dicunits['massinmsun'] = self.dicunits['mass_in_cgs'] / self.CGS_SOLAR_MASS
self.dicunits['timeinGyr'] = self.dicunits['timeincgs'] / (1e3 * self.SEC_PER_MEGAYEAR)
# G in units used in simulation
self.dicunits['gravity_constant_internal'] = BaseUnitsDefs.CGS_GRAVITY / (
self.dicunits['length_in_cgs'] * self.dicunits['vel_in_cgs'] ** 2 / (self.dicunits['mass_in_cgs']))
[docs] def get_dic_of_component_energy(self, ntypes=6, time_conv_type=None, energy_conv_type=None):
"""
Comp type string amongts: 'total', 'gas', 'DM', 'disc stars', 'bulge stars', 'new stars' and 'all stars'.
"""
fen = energy_stats(self.dirname)
if (np.shape(fen)[1] != 4 * (ntypes + 1)):
raise ValueError('Nb of fields in energy file not consistent with ntypes.')
n_times = np.shape(fen)[0]
time = fen[:, 0]
dic = {'Time': time * BaseUnitsDefs.get_conversion_factor_from_string('time', self.dicunits,
conv_type=time_conv_type),
'eint': np.zeros((n_times, ntypes + 1)), 'epot': np.zeros((n_times, ntypes + 1)),
'ekin': np.zeros((n_times, ntypes + 1))}
for i in range(0, ntypes + 1):
dic['eint'][:, i] = fen[:, 3 * i + 1]
dic['epot'][:, i] = fen[:, 3 * i + 2]
dic['ekin'][:, i] = fen[:, 3 * i + 3]
if np.all(dic['epot'][:, 0]) == 0:
warnings.warn('Potential energy is zero, likely because not computed in simulations, energy conservation '
'and Virial ratio evolution are not computable.')
etot = dic['eint'] + dic['epot'] + dic['ekin']
dic['deltaeone'] = (etot - etot[0]) / etot[0] # likely more correct numerically than after conversion
dic['virrap'] = 2 * dic['ekin'] / np.fabs(dic['epot']) # same
dic['eint'] *= BaseUnitsDefs.get_conversion_factor_from_string('energy', self.dicunits,
conv_type=energy_conv_type)
dic['epot'] *= BaseUnitsDefs.get_conversion_factor_from_string('energy', self.dicunits,
conv_type=energy_conv_type)
dic['ekin'] *= BaseUnitsDefs.get_conversion_factor_from_string('energy', self.dicunits,
conv_type=energy_conv_type)
dic['etot'] = etot * BaseUnitsDefs.get_conversion_factor_from_string('energy', self.dicunits,
conv_type=energy_conv_type)
return dic
[docs] def plot_energy_conservation(self, ntypes=6, time_conv_type=None):
"""
Returns time in Gyr and (E(t)-E(t=0))/E(t=0).
"""
dic = self.get_dic_of_component_energy(ntypes=ntypes, time_conv_type=time_conv_type)
xlab = BaseUnitsDefs.get_label_from_string('time', conv_type=time_conv_type)
plot_1d(dic['Time'], dic['deltaeone'][:, 0], xlab=xlab, ylab=r'$\Delta E / E$')
return dic['Time'], dic['deltaeone']
[docs] def plot_kinetic_and_potential_energy(self, ntypes=6, time_conv_type=None, energy_conv_type=None):
"""
Returns time in gyr, kinetic energy and grav pot energy.
"""
dic = self.get_dic_of_component_energy(ntypes=ntypes, time_conv_type=time_conv_type,
energy_conv_type=energy_conv_type)
xlab = BaseUnitsDefs.get_label_from_string('time', conv_type=time_conv_type)
ylab = BaseUnitsDefs.get_label_from_string('energy', conv_type=energy_conv_type)
plot_1d(dic['Time'], dic['ekin'][:, 0], leglab=r'E$_k$')
plot_1d(dic['Time'], dic['epot'][:, 0], xlab=xlab, ylab=ylab, leglab=r'E$_p$')
return dic['Time'], dic['ekin'][:, 0], dic['epot'][:, 0]
[docs] def plot_virial_ratio(self, ntypes=6, time_conv_type=None):
"""
Returns time in Gyr and Virial ratio = 2*T/abs(W) where T is the kinetic energy and W the grav. potential
energy.
"""
dic = self.get_dic_of_component_energy(ntypes=ntypes, time_conv_type=time_conv_type)
xlab = BaseUnitsDefs.get_label_from_string('time', conv_type=time_conv_type)
plot_1d(dic['Time'], dic['virrap'][:, 0], xlab=xlab, ylab=r'$2 T / |W|$')
return dic['Time'], dic['virrap'][:, 0]