"""
This module contains:
- a class Snap that loads a snapshot of a simulation and has methods to perform various computations on data from the
snapshot.
- a class IterOnSim for a simulation from which snapshots can be obtained."""
import glob
from .coords import *
from .readsnapshot import *
from .units import *
from .classeshists import *
from .funcsplots import *
[docs]class Snap(SimCaracs):
"""This class loads a snapshot of a simulation and has methods to perform various computations on data from the
snapshot.
Parameters
----------
snapshot_directory : str
the directory containing the snapshot file
snapshot_file_name : str
the name of the snapshot file
input_type : {'G1/G2', 'G2_str', 'G4_hdf5'}
Snapshot data type, 'G1/G2' for usual Gadget-2 format.
in_sim : bool, default = `True`
`True` if units are specified in snapshot directory in the case of G2 snapshots
gadget2_header, ntypes:
Only for Gadget1/Gadget2 snapshots types used in Gadget4. See :func:`~galaximview.readgadget.read_gadget`
centre_on_com : bool, default = `False`
`True` to centre positions and velocities on the centre of mass
store_pos_vel : bool, default = `True`
`True` to store positions and velocities before any rotation, useful for viewer
Attributes
----------
data: dictionary
Dictionary read from a functon of :exc:`~galaximview.readgadget`, according to 'input_type'.
listcomps: list of str
List of components that can be used in various methods calls.
listindivcomps: list of str
List of components that can be used in various methods calls and that correspond to only one particle type.
listquantities: list of str
List of quantities that can be used in various methods calls
posini: ndarray
(N,3) Positions as initially read, minus the centre of mass if centre_on_com is set to True.
velini: ndarray
(N,3) Velocities as initially read, minus the velocity of the centre of mass if centre_on_com is set to True.
Raises
------
ValueError
If input type is unknown.
Examples
--------
>>> from galaximview.snap import *
>>> snap = Snap('path/to/simulation', 'snapshot_name')
>>> hm2d, edy, edx = snap.surface_density_2d('total','xpos','ypos')
Notes
-----
Methods often have keywords for units conversion. Possibilities are 'cgs', 'astro' and `None`:
- 'cgs' units are usual cgs units.
- 'astro' units are: lengths in kpc, velocities in km/s, time in Gyr, masses in sun mass, mass surface density \
in sun mass/pc^2, mass volume density in sun mass/pc^3. See \
:meth:`~galaximview.units.BaseUnitsDefs.get_conversion_factor_from_string` for computation of conversion \
factors.
- `None` for no conversion.
"""
def __init__(self, snapshot_directory, snapshot_file_name, input_type='G1/G2', gadget2_header=True, ntypes=6,
multiple_snaps=False, in_sim=True, centre_on_com=False, store_pos_vel=True):
self._sorted_by_ids = False
self._in_sim = in_sim
self.fname = snapshot_directory + '/' + snapshot_file_name
self.input_type = input_type
self.data = read_snapshot(self.fname, input_type=input_type, multiple_snaps=multiple_snaps,
gadget2_header=gadget2_header, ntypes=ntypes)
if self._in_sim:
SimCaracs.__init__(self, snapshot_directory)
else:
self.dicunits = None
npart = self.data['npart']
self.listcomps = ['total']
self.listindivcomps = []
if npart[0]:
self.listindivcomps.append('gas')
if npart[1]:
self.listindivcomps.append('DM')
if npart[2]:
self.listindivcomps.append('disc stars')
if npart[3]:
self.listindivcomps.append('bulge stars')
if npart[4]:
self.listindivcomps.append('new stars')
self.listcomps += self.listindivcomps
if (npart[1] | npart[2]):
self.listcomps.append('old stars')
if (npart[1] | npart[2] | npart[3]):
self.listcomps.append('all stars')
if (npart[0] | npart[1] | npart[2] | npart[3]):
self.listcomps.append('baryons')
self.listquantities = ['xpos', 'ypos', 'zpos', 'xvel', 'yvel', 'zvel']
# additional kinematical/dynamical quantities that may be used
self.listquantities += ['radcyl', 'phi', 'radspher', 'thetaspher', 'v_R', 'v_phi', 'velradspher',
'velthetaspher', 'jang_norm', 'zjang', 'spec_jang_norm', 'spec_zjang']
if (npart[0]):
self.listarrs_gas = ['rho', 'uint', 'hsml', 'nH', 'entropy', 'pressure']
self.listquantities += self.listarrs_gas
else:
self.listarrs_gas = None
# additional quantities usually computed and plotted, additional quantities may be put here
self.listquantities += ['voldens', 'surfdens', 'massinr']
self.listquantities += ['zgas']
if centre_on_com:
self.data['pos'] -= self.get_com(np.arange(0, np.sum(self.data['npart'])))
self.data['vel'] -= self.get_vel_of_com(np.arange(0, np.sum(self.data['npart'])))
if store_pos_vel:
self.posini = self.data['pos']
self.velini = self.data['vel']
self.velini = self.data['vel']
self.listquantities += ['xpos_ini', 'ypos_ini', 'zpos_ini'] # used in viewer
else:
self.posini = None
self.velini = None
self.bool_mask = None
self.ctre_pos = np.zeros((3))
self.ctre_vel = np.zeros((3))
[docs] def info(self):
"""Prints the data dictionnary, components of the snapshot and components that can be used in methods."""
print("The data dictionary is the following:", self.data)
print("The snapshot contains the following individual components (listindivcomps):", self.listindivcomps)
print("The following ensembles of components (listcomps) can be retrieved:", self.listcomps)
[docs] def get_isnap(self):
"""Returns the number of the snapshot, taken from numbers after last underscore of the file name.
Returns
-------
int
Number of the snapshot.
"""
end_str = self.fname.split('_')[-1]
if '.' in end_str:
isnap = int(end_str.split('.')[0])
else:
isnap = int(end_str)
return isnap
[docs] def sort_data_by_id(self):
"""Sorts all data dictionnary by ids (particle identifiers and adds a data['g2_type'] equal to the particle
type. """
if not self._sorted_by_ids:
npart_tot = np.sum(self.data['npart'])
npart_gas = self.data['npart'][0]
npart_new_stars = self.data['npart'][4]
part_g2_type = np.zeros((npart_tot))
for compstr in self.listindivcomps:
type_nb = None
if compstr == 'gas':
type_nb = 0
elif compstr == 'DM':
type_nb = 1
elif compstr == 'disc stars':
type_nb = 2
elif compstr == 'bulge stars':
type_nb = 3
elif compstr == 'new stars':
type_nb = 4
part_g2_type[self.get_indexes_of_component(compstr)] = type_nb
postype = self.data['pos'].dtype
veltype = self.data['vel'].dtype
xpos, ypos, zpos = self.data['pos'].T
xvel, yvel, zvel = self.data['vel'].T
xpos_ini, ypos_ini, zpos_ini = self.posini.T
xvel_ini, yvel_ini, zvel_ini = self.velini.T
dtype_tot = [('xpos', postype), ('ypos', postype), ('zpos', postype), ('xvel', veltype), ('yvel', veltype),
('zvel', veltype), ('xpos_ini', postype), ('ypos_ini', postype), ('zpos_ini', postype),
('xvel_ini', veltype), ('yvel_ini', veltype), ('zvel_ini', veltype),
('g2_type', part_g2_type.dtype)]
dtype_gas = []
dtype_new_stars = []
for k, v in self.data.items():
if np.size(v) == npart_tot:
dtype_tot.append((k, v.dtype))
elif ((np.size(v) == npart_gas) & (npart_gas != npart_tot)):
dtype_gas.append((k, v.dtype))
elif ((np.size(v) == npart_new_stars) & (npart_new_stars != npart_tot)):
dtype_new_stars.append((k, v.dtype))
obj_tot = np.zeros(npart_tot, dtype_tot)
for string, arr in zip(['xpos', 'ypos', 'zpos', 'xvel', 'yvel', 'zvel'],
[xpos, ypos, zpos, xvel, yvel, zvel]):
obj_tot[string] = arr
for string, arr in zip(['xpos_ini', 'ypos_ini', 'zpos_ini', 'xvel_ini', 'yvel_ini', 'zvel_ini'],
[xpos_ini, ypos_ini, zpos_ini, xvel_ini, yvel_ini, zvel_ini]):
obj_tot[string] = arr
obj_tot['g2_type'] = part_g2_type
obj_gas = None
obj_new_stars = None
idstype = self.data['ids'].dtype
if ((npart_gas > 0) & (npart_gas != npart_tot)):
gas_ids = self.data['ids'][self.get_indexes_of_component('gas')]
dtype_gas.append(('ids', idstype))
obj_gas = np.zeros(npart_gas, dtype_gas)
obj_gas['ids'] = gas_ids
if ((npart_new_stars > 0) & (npart_new_stars != npart_tot)):
new_stars_ids = self.data['ids'][self.get_indexes_of_component('new stars')]
dtype_new_stars.append(('ids', idstype))
obj_new_stars = np.zeros(npart_new_stars, dtype_new_stars)
obj_new_stars['ids'] = new_stars_ids
for k, v in self.data.items():
if np.size(v) == npart_tot:
obj_tot[k] = v
elif ((np.size(v) == npart_gas) & (npart_gas != npart_tot)):
obj_gas[k] = v
elif ((np.size(v) == npart_new_stars) & (npart_new_stars != npart_tot)):
obj_new_stars[k] = v
sorted_obj_tot = np.sort(obj_tot, order='ids')
if ((npart_gas > 0) & (npart_gas != npart_tot)):
sorted_obj_gas = np.sort(obj_gas, order='ids')
else:
sorted_obj_gas = None
if ((npart_new_stars > 0) & (npart_new_stars != npart_tot)):
sorted_obj_new_stars = np.sort(obj_new_stars, order='ids')
else:
sorted_obj_new_stars = None
self.data['pos'] = np.vstack((sorted_obj_tot['xpos'], sorted_obj_tot['ypos'], sorted_obj_tot['zpos'])).T
self.data['vel'] = np.vstack((sorted_obj_tot['xvel'], sorted_obj_tot['yvel'], sorted_obj_tot['zvel'])).T
self.posini = np.vstack(
(sorted_obj_tot['xpos_ini'], sorted_obj_tot['ypos_ini'], sorted_obj_tot['zpos_ini'])).T
self.velini = np.vstack(
(sorted_obj_tot['xvel_ini'], sorted_obj_tot['yvel_ini'], sorted_obj_tot['zvel_ini'])).T
for k, v in self.data.items():
if np.size(v) == npart_tot:
self.data[k] = sorted_obj_tot[k]
elif ((np.size(v) == npart_gas) & (npart_gas != npart_tot)):
self.data[k] = sorted_obj_gas[k]
elif ((np.size(v) == npart_new_stars) & (npart_new_stars != npart_tot)):
self.data[k] = sorted_obj_new_stars[k]
self.data['g2_type'] = sorted_obj_tot['g2_type']
self._sorted_by_ids = True
[docs] def reset_pos_vel(self):
"""Can be used to reset positions and velocities of data dictionnary as read in snapshot, minus position and
velocity of centre of masse if centre_on_com is True.
Raises
------
ValueError
if store_pos_vel is False
"""
if self.posini is not None:
self.data['pos'] = self.posini
self.data['vel'] = self.velini
else:
raise ValueError(
"Set store_pos_vel=True when creating Snap object to be able to reset positions and velocities.")
[docs] def change_viewing_angle(self, azim, elev):
"""Projects positions and velocities on the right basis to obtain projection. This actually modifies
self.data['pos'] and self.data['vel']. Untouched positions and velocities are stored in
self.posini and self.velini.
Parameters
----------
azim, elev: float
See `~basefuncs.get_basis_los`
"""
if self.posini is None:
raise Warning(
'store_pos_vel must be set to True (default behaviour) to change viewing angle with this method.')
basis_los = get_basis_los(azim, elev)
self.data['pos'] = BaseTools.project_on_orthonormed_direct_basis(self.posini, basis_los)
self.data['vel'] = BaseTools.project_on_orthonormed_direct_basis(self.velini, basis_los)
[docs] def get_indexes_of_component(self, compstr):
"""Selects indexes of component type of self.listcomps.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
Returns
-------
ndarray
(N,) Indexes of the components selected by compstr.
Raises
------
ValueError
If comsptr is not in self.listcomps.
"""
if self.bool_mask is None:
bool_mask = np.ones(np.sum(self.data['npart']), dtype='bool')
else:
bool_mask = self.bool_mask
if self._sorted_by_ids:
if (compstr == 'total'):
inds = np.ones(np.sum(self.data['npart']), dtype='bool')
elif (compstr == 'gas'):
inds = self.data['g2_type'] == 0
elif (compstr == 'DM'):
inds = self.data['g2_type'] == 1
elif (compstr == 'disc stars'):
inds = self.data['g2_type'] == 2
elif (compstr == 'bulge stars'):
inds = self.data['g2_type'] == 3
elif (compstr == 'new stars'):
inds = self.data['g2_type'] == 4
elif (compstr == 'old stars'):
inds = ((self.data['g2_type'] == 2) | (self.data['g2_type'] == 3))
elif (compstr == 'all stars'):
inds = (self.data['g2_type'] == 2) | (self.data['g2_type'] == 3) | (self.data['g2_type'] == 4)
elif (compstr == 'baryons'):
inds = (self.data['g2_type'] == 0) | (self.data['g2_type'] == 2) | (self.data['g2_type'] == 3) | (
self.data['g2_type'] == 4)
else:
raise ValueError("Component is not present. 'total', 'gas', 'DM', 'disc stars', 'bulge stars', "
"'new stars', 'old stars', 'all stars' and 'baryons' can be tried.")
else:
inds = np.zeros(np.sum(self.data['npart']), dtype='bool')
if (compstr == 'total'):
inds[0:np.sum(self.data['npart'])] = True
elif (compstr == 'gas'):
inds[np.arange(0, np.sum(self.data['npart'][0:1]))] = True
elif (compstr == 'DM'):
inds[np.arange(int(np.sum(self.data['npart'][0:1])), np.sum(self.data['npart'][0:2]))] = True
elif (compstr == 'disc stars'):
inds[np.arange(int(np.sum(self.data['npart'][0:2])), np.sum(self.data['npart'][0:3]))] = True
elif (compstr == 'bulge stars'):
inds[np.arange(int(np.sum(self.data['npart'][0:3])), np.sum(self.data['npart'][0:4]))] = True
elif (compstr == 'new stars'):
inds[np.arange(int(np.sum(self.data['npart'][0:4])), np.sum(self.data['npart'][0:5]))] = True
elif (compstr == 'old stars'):
inds[np.arange(int(np.sum(self.data['npart'][0:2])), np.sum(self.data['npart'][0:4]))] = True
elif (compstr == 'all stars'):
inds[np.arange(int(np.sum(self.data['npart'][0:2])), np.sum(self.data['npart'][0:5]))] = True
elif (compstr == 'baryons'):
inds[np.arange(0, np.sum(self.data['npart'][0:1]))] = True
inds[np.arange(int(np.sum(self.data['npart'][0:2])), np.sum(self.data['npart'][0:5]))] = True
else:
raise ValueError("Component is not present. 'total', 'gas', 'DM', 'disc stars', 'bulge stars', "
"'new stars', 'old stars', 'all stars' and 'baryons' can be tried.")
inds = np.where(inds & bool_mask)[0]
return inds
[docs] def get_ids_from_indexes(self, inds):
"""Returns the ids (Gadget-2 identifiers integers) of the given indexes.
Parameters
----------
inds : ndarray
(N,) array of indexes.
Returns
-------
ndarray
(N,) array of particle identifiers corresponding to indexes.
"""
return self.data['ids'][inds]
[docs] def get_indexes_from_ids(self, ids):
"""Returns the indexes corresponding to the given ids (Gdget-2 integer identifiers). Does not keep the order of
input ids ndarray.
Parameters
----------
ids : ndarray
(N,) Particle identifiers.
Returns
-------
ndarray
(N,) Indexes corresponding to identifiers. Does not keep the order of input ids ndarray.
"""
return np.where(np.in1d(self.data['ids'], ids))[0]
[docs] def get_array_from_string(self, arrstr, inds, conv_type=None):
"""Associates array of data to the given string belonging to self.listquantities, for the given indexes.
Not defined as an attribute dictionary because some arrays are not created by default (e.g. cylindrical radius).
Parameters
----------
arrstr : str
Member of :attr:`~galaximview.snap.Snap.listquantities`.
inds : ndarray
(N,) indexes to select
conv_type : {None, 'astro', 'cgs'}
Type of units conversion.
Returns
-------
ndarray
Array obtained from string identification.
"""
if (self.listarrs_gas is not None):
if (arrstr in self.listarrs_gas):
indsgas = self.get_indexes_of_component('gas')
inds = np.searchsorted(indsgas, inds)
if (arrstr == 'xpos'):
arr = self.data['pos'][inds, 0]
elif (arrstr == 'ypos'):
arr = self.data['pos'][inds, 1]
elif (arrstr == 'zpos'):
arr = self.data['pos'][inds, 2]
elif (arrstr == 'xvel'):
arr = self.data['vel'][inds, 0]
elif (arrstr == 'yvel'):
arr = self.data['vel'][inds, 1]
elif (arrstr == 'zvel'):
arr = self.data['vel'][inds, 2]
elif (arrstr == '-zvel'):
arr = -self.data['vel'][inds, 2]
elif (arrstr == 'xpos_ini'):
arr = self.posini[inds, 0]
elif (arrstr == 'ypos_ini'):
arr = self.posini[inds, 1]
elif (arrstr == 'zpos_ini'):
arr = self.posini[inds, 2]
elif (arrstr == 'uint'):
arr = self.data['uint'][inds]
elif (arrstr == 'hsml'):
arr = self.data['hsml'][inds]
elif (arrstr == 'rho'):
arr = self.data['rho'][inds]
elif (arrstr == 'nH'):
arr = self.nhnucl(inds)
elif (arrstr == 'entropy'):
arr = self.entropy(inds)
elif (arrstr == 'pressure'):
arr = self.pressure(inds)
elif (arrstr == 'pot'):
arr = self.data['pot'][inds]
elif (arrstr == 'met'):
arr = self.data['met'][inds]
elif (arrstr == 'temperature'):
arr = self.get_temp_in_kelvin()[inds]
elif (arrstr == 'radcyl'):
arr = SpatialCoordsFromCart(self.data['pos'][inds, :], ctre_pos=self.ctre_pos).cylindrical_radius()
elif (arrstr == 'phi'):
arr = SpatialCoordsFromCart(self.data['pos'][inds, :], ctre_pos=self.ctre_pos).phi()
elif (arrstr == 'radspher'):
arr = SpatialCoordsFromCart(self.data['pos'][inds, :], ctre_pos=self.ctre_pos).spherical_radius()
elif (arrstr == 'thetaspher'):
arr = SpatialCoordsFromCart(self.data['pos'][inds, :], ctre_pos=self.ctre_pos).spherical_theta()
elif (arrstr == 'v_R'):
arr = VelCoordsFromCart(self.data['pos'][inds, :], self.data['vel'][inds, :],
ctre_pos=self.ctre_pos, ctre_vel=self.ctre_vel).cylindrical_radial_vel()
elif (arrstr == 'v_phi'):
arr = VelCoordsFromCart(self.data['pos'][inds, :], self.data['vel'][inds, :],
ctre_pos=self.ctre_pos, ctre_vel=self.ctre_vel).phi_vel()
elif (arrstr == 'velradspher'):
arr = VelCoordsFromCart(self.data['pos'][inds, :], self.data['vel'][inds, :],
ctre_pos=self.ctre_pos, ctre_vel=self.ctre_vel).spherical_radial_vel()
elif (arrstr == 'velthetaspher'):
arr = VelCoordsFromCart(self.data['pos'][inds, :], self.data['vel'][inds, :],
ctre_pos=self.ctre_pos, ctre_vel=self.ctre_vel).spherical_theta_vel()
elif (arrstr == 'jang_norm'): # norm of jang
arr = DynamicsFromCart(self.data['pos'][inds, :], self.data['vel'][inds, :],
self.data['mass'][inds], ctre_pos=self.ctre_pos,
ctre_vel=self.ctre_vel).angular_momentum_norm()
elif (arrstr == 'zjang'):
arr = DynamicsFromCart(self.data['pos'][inds, :], self.data['vel'][inds, :],
self.data['mass'][inds], ctre_pos=self.ctre_pos,
ctre_vel=self.ctre_vel).angular_momentum_3d()
arr = arr[:, 2]
elif (arrstr == 'spec_jang_norm'): # norm of jang
arr = DynamicsFromCart(self.data['pos'][inds, :], self.data['vel'][inds, :],
self.data['mass'][inds], ctre_pos=self.ctre_pos,
ctre_vel=self.ctre_vel).specific_angular_momentum_norm()
elif (arrstr == 'spec_zjang'):
arr = DynamicsFromCart(self.data['pos'][inds, :], self.data['vel'][inds, :],
self.data['mass'][inds], ctre_pos=self.ctre_pos,
ctre_vel=self.ctre_vel).specific_angular_momentum_3d()
arr = arr[:, 2]
elif (arrstr == 'zgas'):
arr = self.data['zgas'][inds]
elif (arrstr == 'zstars'):
arr = self.data['zstars'][inds]
else:
raise ValueError("Array string is not present, all values of self.listquantities can be tried.")
arr *= BaseUnitsDefs.get_conversion_factor_from_string(arrstr, self.dicunits, conv_type=conv_type)
return arr
[docs] def get_dic_of_components_masses(self, conv_type='astro'):
"""Returns masses of components in code units and also in sun masses if self.insim=True.
Parameters
----------
conv_type : {'astro', 'cgs', None}
Type of units conversion.
Returns
-------
dict
Dictionary of total masses of components in self.listcomps.
"""
dic = {}
for comptypestr in self.listcomps:
dic[comptypestr] = np.sum(self.data['mass'][self.get_indexes_of_component(comptypestr)])
if self._in_sim:
fac = BaseUnitsDefs.get_conversion_factor_from_string('mass', self.dicunits, conv_type=conv_type)
dicinmsun = {}
for k, v in dic.items():
dicinmsun[k + ' in m sun'] = v * fac
dic = {**dic, **dicinmsun}
return dic
[docs] def get_com(self, inds):
"""Returns the centre of mass of the particles having the given indexes.
Parameters
----------
inds : ndarray
(N,) indexes to select
Returns
-------
ndarray
(3,) Position of centre of mass.
"""
com = calc_com(self.data['pos'][inds, :], self.data['mass'][inds])
return com
[docs] def get_ang_mom(self, inds):
"""Returns angular momentum of the particles having the given indexes.
Parameters
----------
inds : ndarray
(N,) indexes to select
Returns
-------
ndarray
(3,) Angular momentum.
"""
dyn = DynamicsFromCart(self.data['pos'][inds, :], self.data['vel'][inds, :], self.data['mass'][inds])
return dyn.total_angular_momentum_3d()
[docs] def get_spin(self, inds):
"""Returns spin (angular momentum of the particles having the given indexes, normalized to unity).
Parameters
----------
inds : ndarray
(N,) indexes to select
Returns
-------
ndarray
(3,) Angular momentum.
"""
dyn = DynamicsFromCart(self.data['pos'][inds, :], self.data['vel'][inds, :], self.data['mass'][inds])
return dyn.spin()
[docs] def get_vel_of_com(self, inds):
"""Returns the velocity of the center of mass of the particles having the given indexes.
Parameters
----------
inds : ndarray
(N,) indexes to select
Returns
-------
ndarray
(3,) Velocity of centre of mass.
"""
com = calc_com(self.data['vel'][inds, :], self.data['mass'][inds])
return com
[docs] def select_gas_indexes_with_u_nh(self, ulimits, nhlimits, conv_type=None):
"""Select gas indexes that have internal energy and number density in given limits.
Parameters
----------
ulimits : (float, float)
limits of internal energy
nhlimits : (float, float)
limits of number density
conv_type : {None, 'astro', 'cgs'}
Type of units conversion.
Returns
-------
ndarray
(N,) Gas particles indexes that have internal energy in ulimits and number density in nhlimits.
"""
indsgasinall = self.get_indexes_of_component('gas')
nhnucl_in_units = self.nhnucl(self.get_indexes_of_component('gas'))
nhnucl_in_units *= BaseUnitsDefs.get_conversion_factor_from_string('nH', self.dicunits, conv_type=conv_type)
u_in_units = self.data['uint'] * BaseUnitsDefs.get_conversion_factor_from_string('uint', self.dicunits,
conv_type=conv_type)
inds = (np.log10(nhnucl_in_units) >= nhlimits[0]) & (np.log10(nhnucl_in_units) <= nhlimits[1]) & (
np.log10(u_in_units) >= ulimits[0]) & (np.log10(u_in_units) <= ulimits[1])
inds = np.where(inds)
return indsgasinall[inds]
[docs] def select_gas_indexes_with_temp_nh(self, ulimits, nhlimits, conv_type=None):
"""Select gas indexes that have internal energy and number density in given limits.
Parameters
----------
ulimits : (float, float)
limits of internal energy
nhlimits : (float, float)
limits of number density
conv_type : {None, 'astro', 'cgs'}
Type of units conversion.
Returns
-------
ndarray
(N,) Gas particles indexes that have internal energy in ulimits and number density in nhlimits.
"""
indsgasinall = self.get_indexes_of_component('gas')
nhnucl_in_units = self.nhnucl(self.get_indexes_of_component('gas'))
nhnucl_in_units *= BaseUnitsDefs.get_conversion_factor_from_string('nH', self.dicunits, conv_type=conv_type)
temp_in_units = self.get_temp_in_kelvin()
inds = (np.log10(nhnucl_in_units) >= nhlimits[0]) & (np.log10(nhnucl_in_units) <= nhlimits[1]) & (
np.log10(temp_in_units) >= ulimits[0]) & (np.log10(temp_in_units) <= ulimits[1])
inds = np.where(inds)
return indsgasinall[inds]
[docs] def select_new_stars_indexes_with_tbirth(self, tbirthlimits, conv_type='astro'):
"""Select indexes of new stars that have a time of birth between limits.
Parameters
----------
tbirthlimits : (float, float)
limits of birth time
conv_type : {'astro', 'cgs', None}
Type of units conversion.
Returns
-------
ndarray
(N,) Indexes of new stars that have a time of birth in tbirthlimits.
"""
indsnewstarsinall = self.get_indexes_of_component('new stars')
tbirth = self.data['stellage']
tbirth_in_units = tbirth * BaseUnitsDefs.get_conversion_factor_from_string('time', self.dicunits,
conv_type=conv_type)
inds = (tbirth_in_units >= tbirthlimits[0]) & (tbirth_in_units <= tbirthlimits[1])
inds = np.where(inds)
return indsnewstarsinall[inds]
[docs] def nhnucl(self, inds):
"""Returns the number of H nuclei per unit of volume, in internal units.
Parameters
----------
inds : ndarray
(Ngas,) indexes to select
Returns
-------
ndarray
(Ngas,) number of H nuclei per unit of volume, in internal units.
"""
indsgas = self.get_indexes_of_component('gas')
inds = np.searchsorted(indsgas, inds)
return self.data['rho'][inds] * BaseUnitsDefs.HYDROGEN_MASS_FRACTION * self.dicunits[
'mass_in_cgs'] / BaseUnitsDefs.CGS_PROTONMASS
[docs] def entropy(self, inds):
"""Returns the specific entropy, in internal units.
Parameters
----------
inds : ndarray
(Ngas,) indexes to select
Returns
-------
ndarray
(Ngas,) specific entropy, in internal units.
"""
indsgas = self.get_indexes_of_component('gas')
inds = np.searchsorted(indsgas, inds)
return (BaseUnitsDefs.GAMMA_MONOATOMIC_GAS - 1.) * self.data['uint'][inds] / self.data['rho'][inds] ** (
BaseUnitsDefs.GAMMA_MONOATOMIC_GAS - 1.)
[docs] def pressure(self, inds):
"""Returns the pressure, in internal units.
Parameters
----------
inds : ndarray
(Ngas,) indexes to select
Returns
-------
ndarray
(Ngas,) pressure, in internal units.
"""
indsgas = self.get_indexes_of_component('gas')
inds = np.searchsorted(indsgas, inds)
return (BaseUnitsDefs.GAMMA_MONOATOMIC_GAS - 1.) * self.data['rho'][inds] * self.data['uint'][inds]
[docs] def get_temp_in_kelvin(self):
"""Computes temperatrue in Kelvin, either from electron abundance, if present, or from some approximation."""
temp = np.zeros_like(self.data['uint'])
temp_on_mu = self.data['uint'] * (BaseUnitsDefs.GAMMA_MONOATOMIC_GAS - 1.) * self.dicunits[
'vel_in_cgs'] ** 2 * BaseUnitsDefs.CGS_PROTONMASS / BaseUnitsDefs.CGS_BOLTZMANN
if 'electron_ab' in self.data.keys():
# y = n_He/n_H
yhe = (1 - BaseUnitsDefs.HYDROGEN_MASS_FRACTION) / (4 * BaseUnitsDefs.HYDROGEN_MASS_FRACTION)
mean_molecular_weight = (1 + 4 * yhe) / (1 + yhe + self.data['electron_ab'])
temp = mean_molecular_weight * temp_on_mu
else: # very approximate treatment
# equivalent to (1+4y)/(1+y)
mean_molecular_weight_ionized = 4 / (
3 + 5 * BaseUnitsDefs.HYDROGEN_MASS_FRACTION) # assumes no H_2 and full ionization of H and He
# fh2 = self.data['h2mfrac'] / 2.
fh2 = 0
# equivalent to (1+4y)/(2+3y)
mean_molecular_weight_neutral = 4 / (
1 + 3 * BaseUnitsDefs.HYDROGEN_MASS_FRACTION - 4 * fh2 * BaseUnitsDefs.HYDROGEN_MASS_FRACTION)
temp_if_neutral = mean_molecular_weight_neutral * temp_on_mu
temp_if_ionized = mean_molecular_weight_ionized * temp_on_mu
temp_ionization = 1e4 # 1e4 K
neut = (temp_if_neutral <= temp_ionization)
inter = ((temp_if_neutral > temp_ionization) & (temp_if_ionized < temp_ionization))
ion = (temp_if_ionized >= temp_ionization)
if np.any(neut):
temp[neut] = temp_if_neutral[neut]
if np.any(inter):
temp[inter] = temp_ionization
if np.any(ion):
temp[ion] = temp_if_ionized[ion]
return temp
[docs] def plot_3d_for_subset(self, compstr, frac, xstr='xpos', ystr='ypos', zstr='zpos', xlimit=(-40, 40), ylimit=None,
zlimit=None, x_conv_type='astro', y_conv_type='astro', z_conv_type='astro', inds=None,
ids=None, seed=None, basis=None, zero_point=np.array([0, 0, 0]), ax=None, col=None,
leglab=None):
"""Plots a subset of particles of compstr (and returns the indexes of the subset). A seed can be provided for
identical reproduction.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
frac : float
Fraction of particles plotted.
xlimit : (float, float)
Limits of x axis in units chosen by x_conv_type.
ylimit : (float, float) or None, default = None
Limits of y axis in units chosen by y_conv_type. If None, set as same as xlimit.
zlimit : (float, float) or None, default = None
Limits of z axis in units chosen by z_conv_type. If None, set as same as xlimit.
xstr : str, default = 'xpos'
Member of :attr:`~galaximview.snap.Snap.listquantities`.
ystr : str, default = 'ypos'
Member of :attr:`~galaximview.snap.Snap.listquantities`.
zstr : str, default = 'zpos'
Member of :attr:`~galaximview.snap.Snap.listquantities`.
x_conv_type, y_conv_type, z_conv_type : {'astro', 'cgs', None}
Type of units conversion.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
seed : int, optional
Seed for random generator, can be set to specific natural integer for reproduction.
ax : matplotlib axis, optional
Specific matplotlib axis to plot on.
basis : ndarray or None, optional
(3,3) basis to project positions on.
zero_point : ndarray, optional
(3,) Removed from x, y, and z values for centering, after potential projection on basis. Defaults
to np.array([0, 0, 0]).
col, leglab :
See :func:`~galaximview.funcsplots.plot_3d`
Returns
-------
ndarray
(N,) Indexes of the subset.
Examples
--------
>>> from galaximview.snap import *
>>> snap = Snap('path/to/simulation', 'snapshot_name')
>>> indexes = snap.plot_3d_for_subset('all stars',0.1,xstr='xpos',ystr='ypos',zstr='zpos',xlimit=(-100,100))
"""
inds = self.get_indexes_from_user_input(compstr, inds, ids)
subsetindexes = BaseTools.subset_of_indexes(inds, frac, seed=seed)
xcoords, ycoords, zcoords = [self.get_array_from_string(x, subsetindexes, conv_type=coord_conv_type) for
x, coord_conv_type in
zip([xstr, ystr, zstr], [x_conv_type, y_conv_type, z_conv_type])]
if basis is not None:
pos = np.vstack((xcoords, ycoords, zcoords)).T
pos = BaseTools.project_on_orthonormed_direct_basis(pos, basis)
xcoords = pos[:, 0]
ycoords = pos[:, 1]
zcoords = pos[:, 2]
xlab, ylab, zlab = [self.get_label_from_string(x, conv_type=coord_conv_type) for x, coord_conv_type in
zip([xstr, ystr, zstr], [x_conv_type, y_conv_type, z_conv_type])]
if ylimit is None:
ylimit = xlimit
if zlimit is None:
zlimit = xlimit
plot_3d(xcoords - zero_point[0], ycoords - zero_point[1], zcoords - zero_point[2], xlimit, ylimit, zlimit,
xlab=xlab, ylab=ylab, zlab=zlab, linestyle=',', ax=ax, color=col, leglab=leglab)
return subsetindexes
[docs] def surface_density_of_radius(self, compstr, inds=None, ids=None, funcofx=identity_function, rad_conv_type='astro',
xrange=(None, None), nbins=128, edges=None, fill_range=False,
surfdens_conv_type='astro', show_plot=False, funcofyforpl=np.log10,
xlimit_pl=None, ylimit_pl=None, linestyle='-', linewidth=1, col=None, leglab=None):
"""Surface density as a function of cylindrical radius. Returns bins edges, middle-points and histogram values.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function of radius used for x axis. Defaults to no transformation.
rad_conv_type, surfdens_conv_type : {'astro', 'cgs', None}
Type of units conversion.
xrange, nbins, edges, fill_range : optional
See :class:`~galaximview.classeshists.HistOfX`
show_plot : bool, default = ` False`
`True` to show plot.
funcofyforpl : ufunc, default = np.log10
Function of surface density to plot.
xlimit_pl, ylimit_pl, linestyle, linewidth, col, leglab:
See :func:`~galaximview.funcsplots.plot_1d`
Returns
-------
ndarray, ndarray
(N,) Surface density of cylindrical radius histogram and (N+1,) cylindrical radius bins edges.
"""
inds = self.get_indexes_from_user_input(compstr, inds, ids)
rad = SpatialCoordsFromCart(self.data['pos'][inds, :], ctre_pos=self.ctre_pos).cylindrical_radius()
hofx = HistOfX(rad, xrange, nbins=nbins, funcofx=funcofx, edges=edges, fill_range=fill_range)
x, h = hofx.volumic_hist_of_f_of_x(weights=self.data['mass'][inds], volstr='annulus_surface')
xlab = self.get_label_from_string('radcyl', func_for_label=funcofx, conv_type=rad_conv_type)
ylab = self.get_label_from_string('surfdens', func_for_label=funcofyforpl, conv_type=surfdens_conv_type)
h *= BaseUnitsDefs.get_conversion_factor_from_string('surfdens', self.dicunits, conv_type=surfdens_conv_type)
if show_plot:
plot_1d(x, h, xlimit_pl=xlimit_pl, ylimit_pl=ylimit_pl, xlab=xlab, ylab=ylab, leglab=leglab,
funcofy=funcofyforpl, linestyle=linestyle, linewidth=linewidth, color=col)
return x, h
[docs] def volume_density_of_radius(self, compstr, inds=None, ids=None, funcofx=identity_function, rad_conv_type='astro',
xrange=(None, None), nbins=128, edges=None, fill_range=False,
voldens_conv_type='astro', show_plot=False, funcofyforpl=np.log10,
xlimit_pl=None, ylimit_pl=None, linestyle='-', linewidth=1, col=None, leglab=None):
"""Volume density as a function of spherical radius. Returns bins middleR-points and histogram values.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx : ufunc, default =
Function of radius used for x axis. Defaults to identity_function (no transfoormation).
rad_conv_type, voldens_conv_type : {'astro', 'cgs', None}
Type of units conversion.
xrange, nbins, edges, fill_range : optional
See :class:`~galaximview.classeshists.HistOfX`.
show_plot : bool, default = ` False`
`True` to show plot.
funcofyforpl : ufunc, default = np.log10
Function of volume density to plot.
xlimit_pl, ylimit_pl, linestyle, linewidth, col, leglab: optional
See :func:`~galaximview.funcsplots.plot_1d`
Returns
-------
ndarray, ndarray
(N,) Volume density of spherical radius histogram and (N+1,) spherical radius bins edges.
"""
inds = self.get_indexes_from_user_input(compstr, inds, ids)
rad = SpatialCoordsFromCart(self.data['pos'][inds, :], ctre_pos=self.ctre_pos).spherical_radius()
x, h = HistOfX(rad, xrange, nbins=nbins, funcofx=funcofx, edges=edges,
fill_range=fill_range).volumic_hist_of_f_of_x(weights=self.data['mass'][inds],
volstr='spherical_shell')
xlab = self.get_label_from_string('radspher', func_for_label=funcofx, conv_type=rad_conv_type)
ylab = self.get_label_from_string('voldens', func_for_label=funcofyforpl, conv_type=voldens_conv_type)
h *= BaseUnitsDefs.get_conversion_factor_from_string('voldens', self.dicunits, conv_type=voldens_conv_type)
if show_plot:
plot_1d(x, h, xlimit_pl=xlimit_pl, ylimit_pl=ylimit_pl, xlab=xlab, ylab=ylab, leglab=leglab,
funcofy=funcofyforpl, linestyle=linestyle, linewidth=linewidth, color=col)
return x, h
[docs] def mass_in_radius(self, compstr, radstr, inds=None, ids=None, funcofx=identity_function, rad_conv_type='astro',
xrange=(None, None), nbins=128, edges=None, fill_range=False, mass_conv_type='astro',
show_plot=False, funcofyforpl=identity_function, xlimit_pl=None, ylimit_pl=None, linestyle='-',
linewidth=1, col=None, leglab=None):
"""Mass inside radius (to be chosen between cylindrical or spherical). Returns bins middle-points and histogram
values.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
radstr : str
radcyl' for cylindrucal radius or 'radspher' for spherical radius.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function of radius used for x axis. Defaults to identity_function (no transfoormation).
rad_conv_type, mass_conv_type : {'astro', 'cgs', None}
Type of units conversion.
xrange, nbins, edges, fill_range : optional
See :class:`~galaximview.classeshists.HistOfX`
show_plot : bool, default = ` False`
`True` to show plot.
funcofyforpl : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function of mass in radius to plot. Defaults to no transformation.
xlimit_pl, ylimit_pl, linestyle, linewidth, col, leglab: optional
See :func:`~galaximview.funcsplots.plot_1d`
Returns
-------
ndarray, ndarray
(N,) Mass inside radius histogram and (N+1,) radius bins edges.
"""
inds = self.get_indexes_from_user_input(compstr, inds, ids)
rad = self.get_array_from_string(radstr, inds, conv_type=rad_conv_type)
histofr = HistOfX(rad, xrange, nbins=nbins, funcofx=funcofx, edges=edges, fill_range=fill_range)
x, h = histofr.hist_of_f_of_x(weights=self.data['mass'][inds])
massinrmin = np.sum(self.data['mass'][inds[0] + np.where(rad < histofr.edges[0])])
h0 = massinrmin
h = np.concatenate((np.array([h0]), h))
h = np.cumsum(h)
r = histofr.edges
xlab = self.get_label_from_string(radstr, func_for_label=funcofx, conv_type=rad_conv_type)
ylab = self.get_label_from_string('massinr', func_for_label=funcofyforpl, conv_type=mass_conv_type)
h *= BaseUnitsDefs.get_conversion_factor_from_string('massinr', self.dicunits, conv_type=mass_conv_type)
if show_plot:
plot_1d(r, h, xlimit_pl=xlimit_pl, ylimit_pl=ylimit_pl, xlab=xlab, ylab=ylab, leglab=leglab,
funcofy=funcofyforpl, linestyle=linestyle, linewidth=linewidth, color=col)
return r, h
[docs] def weighted_hist_of_array(self, compstr, xstr, inds=None, ids=None, funcofx=identity_function, x_conv_type='astro',
mass_conv_type='astro', xrange=(None, None), nbins=128, edges=None, fill_range=False,
weight_type='mass', show_plot=False, funcofyforpl=identity_function, xlimit_pl=None,
ylimit_pl=None, linestyle='-', linewidth=1, col=None, leglab=None):
"""Returns bins middle-points and histogram values.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
xstr : str
Member of :attr:`~galaximview.snap.Snap.listquantities`.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function of array. Defaults to no transfoormation.
x_conv_type, mass_conv_type : {'astro', 'cgs', None}
Type of units conversion.
xrange, nbins, edges, fill_range : optional
See :class:`~galaximview.classeshists.HistOfX`
weight_type : {'mass', 'number'}
Type of weighting for 2D histogramm: 'mass' for mass-weighting, 'number' for number-weighting.
show_plot : bool, default = ` False`
`True` to show plot.
funcofyforpl : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function of weighted histogram to plot. Defaults to no transformation.
xlimit_pl, ylimit_pl, linestyle, linewidth, col, leglab:
See :func:`~galaximview.funcsplots.plot_1d`
Returns
-------
ndarray, ndarray
(N,) Weighted histogram of quantity and (N+1,) bins edges.
"""
inds = self.get_indexes_from_user_input(compstr, inds, ids)
xcoords = self.get_array_from_string(xstr, inds, conv_type=x_conv_type)
hofx = HistOfX(xcoords, xrange, nbins=nbins, funcofx=funcofx, edges=edges, fill_range=fill_range)
if weight_type == 'mass':
mass = self.data['mass'][inds]
x, hv = hofx.hist_of_f_of_x(weights=mass)
hv *= BaseUnitsDefs.get_conversion_factor_from_string('mass', self.dicunits, conv_type=mass_conv_type)
elif weight_type == 'number':
x, hv = hofx.hist_of_f_of_x()
else:
raise ValueError('weight_type to be chosen between mass and number.')
xlab = self.get_label_from_string(xstr, func_for_label=funcofx, conv_type=x_conv_type)
if show_plot:
plot_1d(x, hv, xlimit_pl=xlimit_pl, ylimit_pl=ylimit_pl, xlab=xlab, leglab=leglab, funcofy=funcofyforpl,
linestyle=linestyle, linewidth=linewidth, color=col)
return x, hv
[docs] def binned_quantity_in_1d(self, compstr, xstr, ystr, inds=None, ids=None, funcofx=identity_function,
x_conv_type='astro',
y_conv_type='astro', xrange=(None, None), nbins=128, edges=None, fill_range=False,
weight_type='mass', show_plot=False, funcofyforpl=identity_function, xlimit_pl=None,
ylimit_pl=None, linestyle='-', linewidth=1, col=None, leglab=None):
"""Returns bins middle-points and histogram values.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
xstr, ystr : str
Members of :attr:`~galaximview.snap.Snap.listquantities` . Typically 'xpos', 'ypos' or 'zpos'.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function of array. Defaults to no transformation.
x_conv_type, y_conv_type : {'astro', 'cgs', None}
Type of units conversion for specified axis.
xrange, nbins, edges, fill_range : optional
See :class:`~galaximview.classeshists.HistOfX`
weight_type : {'mass', 'number'}
Type of weighting for 2D histogramm: 'mass' for mass-weighting, 'number' for number-weighting.
show_plot : bool, default = ` False`
`True` to show plot.
funcofyforpl : ufunc
Function of weighted histogram to plot. Defaults to no transformation.
xlimit_pl, ylimit_pl, linestyle, linewidth, col, leglab:
See :func:`~galaximview.funcsplots.plot_1d`
Returns
-------
ndarray, ndarray
(N,) Weighted histogram of y quantity at x and (N+1,) bins edges.
"""
inds = self.get_indexes_from_user_input(compstr, inds, ids)
xcoords, ycoords = [self.get_array_from_string(x, inds, conv_type=coord_conv_type) for x, coord_conv_type in
zip([xstr, ystr], [x_conv_type, y_conv_type])]
hofx = HistOfX(xcoords, xrange, nbins=nbins, funcofx=funcofx, edges=edges, fill_range=fill_range)
if weight_type == 'mass':
x, hv = hofx.binned_y_of_f_of_x(ycoords, weights=self.data['mass'][inds])
elif weight_type == 'number':
x, hv = hofx.binned_y_of_f_of_x(ycoords, weights=None)
else:
raise ValueError('weight_type to be chosen between mass and number.')
xlab = self.get_label_from_string(xstr, func_for_label=funcofx, conv_type=x_conv_type)
ylab = self.get_label_from_string(ystr, func_for_label=funcofyforpl, conv_type=y_conv_type)
if show_plot:
plot_1d(x, hv, xlimit_pl=xlimit_pl, ylimit_pl=ylimit_pl, xlab=xlab, ylab=ylab, leglab=leglab,
funcofy=funcofyforpl, linestyle=linestyle, linewidth=linewidth, color=col)
return x, hv
[docs] def velocity_dispersion_of_radius(self, compstr, xstr='radcyl', ystr='v_R', inds=None, ids=None,
funcofx=identity_function, x_conv_type='astro', y_conv_type='astro',
xrange=(None, None), nbins=128, edges=None, fill_range=False, show_plot=False,
funcofyforpl=np.log10, xlimit_pl=None, ylimit_pl=None, linestyle='-',
linewidth=1, col=None, leglab=None):
"""Default: radial velocity dispersion as a function of cylindrical radius. Returns bins middle-points and
histogram values.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
xstr, ystr : str
Members of :attr:`~galaximview.snap.Snap.listquantities` . Typically 'xpos', 'ypos' or 'zpos'.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function of array. Defaults to no transformation.
x_conv_type, y_conv_type : {'astro', 'cgs', None}
Type of units conversion for specified axis.
xrange, nbins, edges, fill_range : optional
See :class:`~galaximview.classeshists.HistOfX`
show_plot : bool, default = `False`
`True` to show plot.
funcofyforpl : ufunc
Function of surface density to plot. Defaults to np.log10.
xlimit_pl, ylimit_pl, linestyle, linewidth, col, leglab:
See :func:`~galaximview.funcsplots.plot_1d`
Returns
-------
ndarray, ndarray
(N,) Velocity dispersion and (N+1,) radius bins edges.
"""
inds = self.get_indexes_from_user_input(compstr, inds, ids)
xcoords = self.get_array_from_string(xstr, inds, conv_type=x_conv_type)
velcoords = self.get_array_from_string(ystr, inds, conv_type=y_conv_type)
hofx = HistOfX(xcoords, xrange, nbins=nbins, funcofx=funcofx, edges=edges, fill_range=fill_range)
x, h = hofx.binned_dispersion_of_y_of_f_of_x(velcoords, weights=self.data['mass'][inds])
xlab = self.get_label_from_string(xstr, func_for_label=funcofx, conv_type=x_conv_type)
ylab = self.get_label_from_string(ystr, func_for_label=funcofyforpl, label_of_dispersion=True,
conv_type=y_conv_type)
if show_plot:
plot_1d(x, h, xlimit_pl=xlimit_pl, ylimit_pl=ylimit_pl, xlab=xlab, ylab=ylab, leglab=leglab,
funcofy=funcofyforpl, linestyle=linestyle, linewidth=linewidth, color=col)
return x, h
[docs] def weighted_hist_2d(self, compstr, xstr, ystr, inds=None, ids=None, funcofx=np.log10, funcofy=np.log10,
x_conv_type='astro', y_conv_type='astro', mass_conv_type=None, xlimit=(-40., 40.), ylimit=None,
nbx=128, nby=None, weight_type='mass', cic_hist=False, show_plot=True, funcforpl=np.log10,
vmin_plot=None, vmax_plot=None, add_eps_in_plot=False, cmap=None, show_color_bar=True,
fig=None, cb_loc=None):
"""Computes and plot 2D mass-weighted hist of xstr and ystr value for component or particles from indexes.
Outputs/plots the mass in each bin.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
xstr, ystr : str
Members of :attr:`~galaximview.snap.Snap.listquantities` . Typically 'xpos', 'ypos' or 'zpos'.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx, funcofy : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function used to create bins on specified dimension. Defaults to no transformation.
x_conv_type, y_conv_type, mass_conv_type : {'astro', 'cgs', None}
Type of units conversion.
xlimit, ylimit, nbx, nby: optional
See :func:`~galaximview.snap.get_hist2d`
weight_type : {'mass', 'number'}
Type of weighting for 2D histogramm: 'mass' for mass-weighting, 'number' for number-weighting.
cic_hist : bool, default = `False`
If `True`, computes cloud in cloud histogram.
show_plot : bool, default = `True`
`True` to show plot.
funcforpl, vmin_plot, vmax_plot, add_eps_in_plot, cmap, show_color_bar, fig, cb_loc: optional
See :func:`~galaximview.funcsplots.plot_hist2d`
Returns
-------
ndarray, ndarray, ndarray
(nby, nbx) Weighted 2D histogram, (nby+1) y axis bins edges and (nbx+1) x axis bins edges.
"""
inds = self.get_indexes_from_user_input(compstr, inds, ids)
xcoords, ycoords = [self.get_array_from_string(x, inds, conv_type=coord_conv_type) for x, coord_conv_type in
zip([xstr, ystr], [x_conv_type, y_conv_type])]
mh2d = get_hist2d(xcoords, ycoords, xlimit, ylimit, nbx, nby, cic_hist=cic_hist, funcofx=funcofx,
funcofy=funcofy)
if weight_type == 'mass':
hm2d, edy, edx = mh2d.hist2d(weights=self.data['mass'][inds])
hm2d *= BaseUnitsDefs.get_conversion_factor_from_string('mass', self.dicunits, conv_type=mass_conv_type)
elif weight_type == 'number':
hm2d, edy, edx = mh2d.hist2d()
else:
raise ValueError('weight_type to be chosen between mass and number.')
xlab, ylab = [self.get_label_from_string(x, conv_type=coord_conv_type, func_for_label=func) for
x, coord_conv_type, func in
zip([xstr, ystr], [x_conv_type, y_conv_type], [funcofx, funcofy])]
if show_plot:
plot_hist2d(hm2d, edy, edx, funcforpl=funcforpl, vmin_plot=vmin_plot, vmax_plot=vmax_plot,
add_eps_in_plot=add_eps_in_plot, cmap=cmap, show_color_bar=show_color_bar, fig=fig,
cb_loc=cb_loc, xlab=xlab, ylab=ylab)
return hm2d, edy, edx
[docs] def weighted_hist_2d_with_marginals(self, compstr, xstr, ystr, inds=None, ids=None, funcofx=np.log10,
funcofy=np.log10, x_conv_type='astro', y_conv_type='astro', mass_conv_type=None,
xlimit=(None, None), ylimit=(None, None), nbx=128, nby=128, weight_type='mass',
funcforpl=np.log10, cmap=None, fig=None, ax=None):
"""
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
xstr, ystr : str
Members of :attr:`~galaximview.snap.Snap.listquantities` . Typically 'xpos', 'ypos' or 'zpos'.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx, funcofy : ufunc, default = np.log10
Function used to create bins on specified dimension.
x_conv_type, y_conv_type, mass_conv_type : {'astro', 'cgs', None}
Type of units conversion for specified axis.
xlimit, ylimit, nbx, nby: optional
See :func:`~galaximview.snap.get_hist2d`
weight_type : {'mass', 'number'}
Type of weighting for 2D histogramm: 'mass' for mass-weighting, 'number' for number-weighting.
funcforpl, cmap, fig, ax:
See :func:`~galaximview.funcsplots.plot_hist2d_with_pdfs`
Returns
-------
(ndarray, ndarray, ndarray, ndarray, ndarray, matplotlib axis, matplotlib axis, matplotlib axis)
(nby, nbx) Weighted 2D histogram, (nbx,) marginal x histogram, (nby,) marginal y histogram,
(nby+1) y axis bins edges, (nbx+1) x axis bins edges and matplolib axes.
"""
hm2d, edy, edx = self.weighted_hist_2d(compstr, xstr, ystr, inds=inds, ids=ids, funcofx=funcofx,
funcofy=funcofy, x_conv_type=x_conv_type, y_conv_type=y_conv_type,
mass_conv_type=mass_conv_type, weight_type=weight_type, xlimit=xlimit,
ylimit=ylimit, nbx=nbx, nby=nby, show_plot=False)
edxb, hmrho = self.weighted_hist_of_array(compstr, xstr, inds=inds, ids=ids, funcofx=funcofx, xrange=xlimit,
nbins=nbx, show_plot=False, x_conv_type=x_conv_type,
weight_type=weight_type)
edyb, hmu = self.weighted_hist_of_array(compstr, ystr, inds=inds, ids=ids, funcofx=funcofy, xrange=ylimit,
nbins=nby, show_plot=False, x_conv_type=y_conv_type,
weight_type=weight_type)
xlab, ylab = [self.get_label_from_string(x, conv_type=coord_conv_type, func_for_label=func) for
x, coord_conv_type, func in
zip([xstr, ystr], [x_conv_type, y_conv_type], [funcofx, funcofy])]
axy, axz, azy = plot_hist2d_with_pdfs(hm2d, hmrho, hmu, edx, edy, funcforpl=funcforpl, funcforplx=funcofx,
funcforply=funcofy, xlab=xlab, ylab=ylab, fig=fig, cmap=cmap, ax=ax)
return hm2d, hmrho, hmu, edy, edx, axy, axz, azy
[docs] def binned_quantity_in_2d_with_marginals(self, compstr, xstr, ystr, tobinstr, inds=None, ids=None,
funcofx=np.log10, funcofy=np.log10, x_conv_type='astro',
y_conv_type='astro', mass_conv_type=None, xlimit=(None, None),
ylimit=(None, None), nbx=128, nby=128, weight_type='mass',
funcforpl=np.log10, cmap=None, fig=None, ax=None):
"""
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
xstr, ystr : str
Members of :attr:`~galaximview.snap.Snap.listquantities` . Typically 'xpos', 'ypos' or 'zpos'.
tobinstr : str
Member of :attr:`~galaximview.snap.Snap.listquantities`.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx, funcofy : ufunc, default = np.log10
Function used to create bins on specified dimension.
x_conv_type, y_conv_type, mass_conv_type : {'astro', 'cgs', None}
Type of units conversion.
xlimit, ylimit, nbx, nby: optional
See :func:`~galaximview.snap.get_hist2d`
weight_type : {'mass', 'number'}
Type of weighting for 2D histogramm: 'mass' for mass-weighting, 'number' for number-weighting.
funcforpl, cmap, fig, ax:
See :func:`~galaximview.funcsplots.plot_hist2d_with_pdfs`
Returns
-------
ndarray, ndarray, ndarray, ndarray, ndarray, matplotlib axis, matplotlib axis, matplotlib axis
(nby, nbx) Weighted 2D histogram, (nbx,) marginal x histogram, (nby,) marginal y histogram,
(nby+1) y axis bins edges, (nbx+1) x axis bins edges and matplolib axes.
"""
hm2d, edy, edx = self.binned_quantity_in_2d(compstr, xstr, ystr, tobinstr, inds=inds, ids=ids,
funcofx=funcofx, funcofy=funcofy, x_conv_type=x_conv_type,
y_conv_type=y_conv_type, xlimit=xlimit, ylimit=ylimit, nbx=nbx,
nby=nby, weight_type=weight_type, quantity_conv_type=mass_conv_type,
show_plot=False)
edxb, hmrho = self.weighted_hist_of_array(compstr, xstr, inds=inds, ids=ids, funcofx=funcofx, xrange=xlimit,
nbins=nbx, show_plot=False, x_conv_type=x_conv_type,
weight_type=weight_type)
edyb, hmu = self.weighted_hist_of_array(compstr, ystr, inds=inds, ids=ids, funcofx=funcofy, xrange=ylimit,
nbins=nby, show_plot=False, x_conv_type=y_conv_type,
weight_type=weight_type)
xlab, ylab = [self.get_label_from_string(x, conv_type=coord_conv_type, func_for_label=func) for
x, coord_conv_type, func in
zip([xstr, ystr], [x_conv_type, y_conv_type], [funcofx, funcofy])]
axy, axz, azy = plot_hist2d_with_pdfs(hm2d, hmrho, hmu, edx, edy, funcforpl=funcforpl, funcforplx=funcofx,
funcforply=funcofy, xlab=xlab, ylab=ylab, fig=fig, cmap=cmap, ax=ax)
return hm2d, hmrho, hmu, edy, edx, axy, axz, azy
[docs] def sph_density_2d(self, compstr, xstr, ystr, inds=None, ids=None, funcofx=identity_function,
funcofy=identity_function, x_conv_type='astro', y_conv_type='astro', xlimit=(-40., 40.),
ylimit=None, nbx=128, nby=None, surfdens_conv_type='astro', nb_row_parts=4,
show_plot=True, funcforpl=np.log10, vmin_plot=None, vmax_plot=None, add_eps_in_plot=False,
cmap=None, show_color_bar=True, fig=None, cb_loc=None):
"""Computes SPH density.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
xstr, ystr : str
Members of :attr:`~galaximview.snap.Snap.listquantities` . Typically 'xpos', 'ypos' or 'zpos'.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx, funcofy : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function used to create bins on specified dimension. Defaults to no transformation.
x_conv_type, y_conv_type, surfdens_conv_type: {'astro', 'cgs', None}
Type of units conversion.
xlimit, ylimit, nbx, nby: optional
See :func:`~galaximview.snap.get_hist2d`
nb_row_parts : int
Number of spatial rows parts used for speeding up of SPH treatment. Defaults to 4.
show_plot : bool, default = `True`
`True` to show plot.
funcforpl, vmin_plot, vmax_plot, add_eps_in_plot, cmap, show_color_bar, fig, cb_loc: optional
See :func:`~galaximview.funcsplots.plot_hist2d`
Returns
-------
ndarray, ndarray, ndarray
(nby, nbx) Surface density histogram, (nby+1,) y axis bins edges and (nbx,) x axis bins edges.
"""
inds = self.get_indexes_from_user_input(compstr, inds, ids)
xcoords, ycoords = [self.get_array_from_string(x, inds, conv_type=coord_conv_type) for x, coord_conv_type in
zip([xstr, ystr], [x_conv_type, y_conv_type])]
mass = self.data['mass'][inds]
indsgas = self.get_indexes_of_component('gas')
inds = np.searchsorted(indsgas, inds)
hsml = self.data['hsml'][inds]
mh2d = get_hist2d(xcoords, ycoords, xlimit, ylimit, nbx, nby, funcofx=funcofx, funcofy=funcofy)
nbx = np.size(mh2d.xedges) - 1
nby = np.size(mh2d.yedges) - 1
xarr = mh2d.get_x_grid()
yarr = mh2d.get_y_grid()
xwidth = mh2d.get_x_binswidth()[0]
ywidth = mh2d.get_y_binswidth()[0]
xsemwidth = 0.5 * xwidth
ysemwidth = 0.5 * ywidth
xsidelength = (mh2d.xedges[-1] - mh2d.xedges[0])
xmidrow = np.zeros((nb_row_parts))
for i in range(0, nb_row_parts):
xmidrow[i] = mh2d.xedges[0] + (1. / (2 * nb_row_parts) + 1. / nb_row_parts * i) * xsidelength
hm2d = np.zeros((nby, nbx))
for iy in range(0, nby):
distytorow = ycoords - yarr[iy, 0]
indstorow = np.where(np.fabs(distytorow) <= (hsml + ysemwidth))[0]
ytorow = ycoords[indstorow]
xtorow = xcoords[indstorow]
hsmltorow = hsml[indstorow]
indstopartrow = []
xtomidpartrow = []
ytomidpartrow = []
hsmltomidpartrow = []
for i in range(0, nb_row_parts):
distxtomid = xtorow - xmidrow[i]
inds = np.where((np.fabs(distxtomid) <= (hsmltorow + 1. / (2 * nb_row_parts) * xsidelength)))[0]
indstopartrow.append(inds)
xtomidpartrow.append(xtorow[inds])
ytomidpartrow.append(ytorow[inds])
hsmltomidpartrow.append(hsmltorow[inds])
for ix in range(0, nbx):
if ((ix == 0) & (iy % 50 == 0)):
print(ix, iy)
ipartrow = int(ix // (nbx / nb_row_parts))
indstomid = indstopartrow[ipartrow]
xtomid = xtomidpartrow[ipartrow]
ytomid = ytomidpartrow[ipartrow]
hsmltomid = hsmltomidpartrow[ipartrow]
distxn = xtomid - xarr[iy, ix]
distyn = ytomid - yarr[iy, ix]
distrn = np.sqrt(distxn ** 2 + distyn ** 2)
indsn = np.where((distrn <= (hsmltomid)))[0]
indspic = \
np.where((np.fabs(distxn) < xsemwidth) & (np.fabs(distyn) < ysemwidth) & (distrn > hsmltomid))[
0]
indsnform = indstorow[indstomid[indsn]]
indspicform = indstorow[indstomid[indspic]]
hsmltomid = hsmltomid[indsn]
distrn = distrn[indsn]
wsph = np.zeros((np.size(indsn)))
indsclose = np.where(distrn <= 0.5 * hsmltomid)
indsfar = np.where((distrn > 0.5 * hsmltomid) & (distrn <= hsmltomid))
distrnclose = distrn[indsclose]
hsmlclose = hsmltomid[indsclose]
wsph[indsclose] = 1 - 6 * (distrnclose / hsmlclose) ** 2 + 6 * (distrnclose / hsmlclose) ** 3
wsph[indsfar] = 2 * (1 - distrn[indsfar] / hsmltomid[indsfar]) ** 3
wsph *= 8 / (np.pi * hsmltomid ** 3)
hm2d[iy, ix] = np.sum(mass[indsnform] * wsph) + np.sum(mass[indspicform])
indsinhist = np.where(
(xcoords <= mh2d.xedges[-1]) & (xcoords >= mh2d.xedges[0]) & (ycoords <= mh2d.yedges[-1]) & (
ycoords >= mh2d.yedges[0]))
hm2d *= np.sum(mass[indsinhist]) / np.sum(hm2d)
edx = mh2d.xedges
edy = mh2d.yedges
hm2d /= (xsemwidth * ysemwidth * 4)
hm2d *= BaseUnitsDefs.get_conversion_factor_from_string('surfdens', self.dicunits, conv_type=surfdens_conv_type)
xlab, ylab = [self.get_label_from_string(x, conv_type=coord_conv_type, func_for_label=func) for
x, coord_conv_type, func in
zip([xstr, ystr], [x_conv_type, y_conv_type], [funcofx, funcofy])]
cblab = self.get_label_from_string('surfdens', func_for_label=funcforpl, conv_type=surfdens_conv_type)
if show_plot:
plot_hist2d(hm2d, edy, edx, funcforpl=funcforpl, vmin_plot=vmin_plot, vmax_plot=vmax_plot,
add_eps_in_plot=add_eps_in_plot, cmap=cmap, show_color_bar=show_color_bar, fig=fig,
cb_loc=cb_loc, xlab=xlab, ylab=ylab, cblab=cblab)
return hm2d, edy, edx
[docs] def sph_binned_quantity_in_2d(self, compstr, xstr, ystr, tobinstr, inds=None, ids=None, funcofx=identity_function,
funcofy=identity_function, x_conv_type='astro', y_conv_type='astro',
xlimit=(-40., 40.),
ylimit=None, nbx=128, nby=None, quantity_conv_type='astro', nb_row_parts=4,
show_plot=True, funcforpl=np.log10, vmin_plot=None, vmax_plot=None,
add_eps_in_plot=False,
cmap=None, show_color_bar=True, fig=None, cb_loc=None):
"""Computes SPH density.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
xstr, ystr : str
Members of :attr:`~galaximview.snap.Snap.listquantities` . Typically 'xpos', 'ypos' or 'zpos'.
tobinstr : str
Member of :attr:`~galaximview.snap.Snap.listquantities`.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx, funcofy : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function used to create bins on specified dimension. Defaults to no transformation.
x_conv_type, y_conv_type, quantity_conv_type: {'astro', 'cgs', None}
Type of units conversion.
xlimit, ylimit, nbx, nby: optional
See :func:`~galaximview.snap.get_hist2d`
nb_row_parts : int
Number of spatial rows parts used for speeding up of SPH treatment. Defaults to 4.
show_plot : bool, default = `True`
`True` to show plot.
funcforpl, vmin_plot, vmax_plot, add_eps_in_plot, cmap, show_color_bar, fig, cb_loc: optional
See :func:`~galaximview.funcsplots.plot_hist2d`
Returns
-------
ndarray, ndarray, ndarray
(nby, nbx) Surface density histogram, (nby+1,) y axis bins edges and (nbx,) x axis bins edges.
"""
inds = self.get_indexes_from_user_input(compstr, inds, ids)
xcoords, ycoords = [self.get_array_from_string(x, inds, conv_type=coord_conv_type) for x, coord_conv_type in
zip([xstr, ystr], [x_conv_type, y_conv_type])]
tobin = self.get_array_from_string(tobinstr, inds, conv_type=None)
mass = self.data['mass'][inds]
indsgas = self.get_indexes_of_component('gas')
inds = np.searchsorted(indsgas, inds)
hsml = self.data['hsml'][inds]
mh2d = get_hist2d(xcoords, ycoords, xlimit, ylimit, nbx, nby, funcofx=funcofx, funcofy=funcofy)
nbx = np.size(mh2d.xedges) - 1
nby = np.size(mh2d.yedges) - 1
xarr = mh2d.get_x_grid()
yarr = mh2d.get_y_grid()
xwidth = mh2d.get_x_binswidth()[0]
ywidth = mh2d.get_y_binswidth()[0]
xsemwidth = 0.5 * xwidth
ysemwidth = 0.5 * ywidth
xsidelength = (mh2d.xedges[-1] - mh2d.xedges[0])
xmidrow = np.zeros((nb_row_parts))
for i in range(0, nb_row_parts):
xmidrow[i] = mh2d.xedges[0] + (1. / (2 * nb_row_parts) + 1. / nb_row_parts * i) * xsidelength
hm2d = np.zeros((nby, nbx))
hm2dtobin = np.zeros((nby, nbx))
for iy in range(0, nby):
distytorow = ycoords - yarr[iy, 0]
indstorow = np.where(np.fabs(distytorow) <= (hsml + ysemwidth))[0]
ytorow = ycoords[indstorow]
xtorow = xcoords[indstorow]
hsmltorow = hsml[indstorow]
indstopartrow = []
xtomidpartrow = []
ytomidpartrow = []
hsmltomidpartrow = []
for i in range(0, nb_row_parts):
distxtomid = xtorow - xmidrow[i]
inds = np.where((np.fabs(distxtomid) <= (hsmltorow + 1. / (2 * nb_row_parts) * xsidelength)))[0]
indstopartrow.append(inds)
xtomidpartrow.append(xtorow[inds])
ytomidpartrow.append(ytorow[inds])
hsmltomidpartrow.append(hsmltorow[inds])
for ix in range(0, nbx):
if ((ix == 0) & (iy % 50 == 0)):
print(ix, iy)
ipartrow = int(ix // (nbx / nb_row_parts))
indstomid = indstopartrow[ipartrow]
xtomid = xtomidpartrow[ipartrow]
ytomid = ytomidpartrow[ipartrow]
hsmltomid = hsmltomidpartrow[ipartrow]
distxn = xtomid - xarr[iy, ix]
distyn = ytomid - yarr[iy, ix]
distrn = np.sqrt(distxn ** 2 + distyn ** 2)
indsn = np.where((distrn <= (hsmltomid)))[0]
indspic = \
np.where((np.fabs(distxn) < xsemwidth) & (np.fabs(distyn) < ysemwidth) & (distrn > hsmltomid))[
0]
indsnform = indstorow[indstomid[indsn]]
indspicform = indstorow[indstomid[indspic]]
hsmltomid = hsmltomid[indsn]
distrn = distrn[indsn]
wsph = np.zeros((np.size(indsn)))
indsclose = np.where(distrn <= 0.5 * hsmltomid)
indsfar = np.where((distrn > 0.5 * hsmltomid) & (distrn <= hsmltomid))
distrnclose = distrn[indsclose]
hsmlclose = hsmltomid[indsclose]
wsph[indsclose] = 1 - 6 * (distrnclose / hsmlclose) ** 2 + 6 * (distrnclose / hsmlclose) ** 3
wsph[indsfar] = 2 * (1 - distrn[indsfar] / hsmltomid[indsfar]) ** 3
wsph *= 8 / (np.pi * hsmltomid ** 3)
hm2d[iy, ix] = np.sum(mass[indsnform] * wsph) + np.sum(mass[indspicform])
hm2dtobin[iy, ix] = np.sum(mass[indsnform] * tobin[indsnform] * wsph) + np.sum(mass[indspicform])
hm2dtobin /= hm2d
edx = mh2d.xedges
edy = mh2d.yedges
hm2dtobin *= BaseUnitsDefs.get_conversion_factor_from_string(tobinstr, self.dicunits,
conv_type=quantity_conv_type)
xlab, ylab = [self.get_label_from_string(x, conv_type=coord_conv_type, func_for_label=func) for
x, coord_conv_type, func in
zip([xstr, ystr], [x_conv_type, y_conv_type], [funcofx, funcofy])]
cblab = self.get_label_from_string(tobinstr, func_for_label=funcforpl, conv_type=quantity_conv_type)
if show_plot:
plot_hist2d(hm2dtobin, edy, edx, funcforpl=funcforpl, vmin_plot=vmin_plot, vmax_plot=vmax_plot,
add_eps_in_plot=add_eps_in_plot, cmap=cmap, show_color_bar=show_color_bar, fig=fig,
cb_loc=cb_loc, xlab=xlab, ylab=ylab, cblab=cblab)
return hm2d, edy, edx
[docs] def sph_binned_dispersion_of_quantity_in_2d(self, compstr, xstr, ystr, tobinstr, inds=None, ids=None,
funcofx=identity_function, funcofy=identity_function,
x_conv_type='astro', y_conv_type='astro', xlimit=(-40., 40.),
ylimit=None, nbx=128, nby=None, quantity_conv_type='astro',
nb_row_parts=4,
show_plot=True, funcforpl=np.log10, vmin_plot=None, vmax_plot=None,
add_eps_in_plot=False,
cmap=None, show_color_bar=True, fig=None, cb_loc=None):
"""Computes SPH density.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
xstr, ystr : str
Members of :attr:`~galaximview.snap.Snap.listquantities` . Typically 'xpos', 'ypos' or 'zpos'.
tobinstr : str
Member of :attr:`~galaximview.snap.Snap.listquantities`.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx, funcofy : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function used to create bins on specified dimension. Defaults to no transformation.
x_conv_type, y_conv_type, quantity_conv_type: {'astro', 'cgs', None}
Type of units conversion.
xlimit, ylimit, nbx, nby: optional
See :func:`~galaximview.snap.get_hist2d`
nb_row_parts : int
Number of spatial rows parts used for speeding up of SPH treatment. Defaults to 4.
show_plot : bool, default = `True`
`True` to show plot.
funcforpl, vmin_plot, vmax_plot, add_eps_in_plot, cmap, show_color_bar, fig, cb_loc: optional
See :func:`~galaximview.funcsplots.plot_hist2d`
Returns
-------
ndarray, ndarray, ndarray
(nby, nbx) Surface density histogram, (nby+1,) y axis bins edges and (nbx,) x axis bins edges.
"""
inds = self.get_indexes_from_user_input(compstr, inds, ids)
xcoords, ycoords = [self.get_array_from_string(x, inds, conv_type=coord_conv_type) for x, coord_conv_type in
zip([xstr, ystr], [x_conv_type, y_conv_type])]
tobin = self.get_array_from_string(tobinstr, inds, conv_type=None)
mass = self.data['mass'][inds]
indsgas = self.get_indexes_of_component('gas')
inds = np.searchsorted(indsgas, inds)
hsml = self.data['hsml'][inds]
mh2d = get_hist2d(xcoords, ycoords, xlimit, ylimit, nbx, nby, funcofx=funcofx, funcofy=funcofy)
nbx = np.size(mh2d.xedges) - 1
nby = np.size(mh2d.yedges) - 1
xarr = mh2d.get_x_grid()
yarr = mh2d.get_y_grid()
xwidth = mh2d.get_x_binswidth()[0]
ywidth = mh2d.get_y_binswidth()[0]
xsemwidth = 0.5 * xwidth
ysemwidth = 0.5 * ywidth
xsidelength = (mh2d.xedges[-1] - mh2d.xedges[0])
xmidrow = np.zeros((nb_row_parts))
for i in range(0, nb_row_parts):
xmidrow[i] = mh2d.xedges[0] + (1. / (2 * nb_row_parts) + 1. / nb_row_parts * i) * xsidelength
hm2d = np.zeros((nby, nbx))
hm2dtobin = np.zeros((nby, nbx))
hm2dtobinsq = np.zeros((nby, nbx))
for iy in range(0, nby):
distytorow = ycoords - yarr[iy, 0]
indstorow = np.where(np.fabs(distytorow) <= (hsml + ysemwidth))[0]
ytorow = ycoords[indstorow]
xtorow = xcoords[indstorow]
hsmltorow = hsml[indstorow]
indstopartrow = []
xtomidpartrow = []
ytomidpartrow = []
hsmltomidpartrow = []
for i in range(0, nb_row_parts):
distxtomid = xtorow - xmidrow[i]
inds = np.where((np.fabs(distxtomid) <= (hsmltorow + 1. / (2 * nb_row_parts) * xsidelength)))[0]
indstopartrow.append(inds)
xtomidpartrow.append(xtorow[inds])
ytomidpartrow.append(ytorow[inds])
hsmltomidpartrow.append(hsmltorow[inds])
for ix in range(0, nbx):
if ((ix == 0) & (iy % 50 == 0)):
print(ix, iy)
ipartrow = int(ix // (nbx / nb_row_parts))
indstomid = indstopartrow[ipartrow]
xtomid = xtomidpartrow[ipartrow]
ytomid = ytomidpartrow[ipartrow]
hsmltomid = hsmltomidpartrow[ipartrow]
distxn = xtomid - xarr[iy, ix]
distyn = ytomid - yarr[iy, ix]
distrn = np.sqrt(distxn ** 2 + distyn ** 2)
indsn = np.where((distrn <= (hsmltomid)))[0]
indspic = \
np.where((np.fabs(distxn) < xsemwidth) & (np.fabs(distyn) < ysemwidth) & (distrn > hsmltomid))[
0]
indsnform = indstorow[indstomid[indsn]]
indspicform = indstorow[indstomid[indspic]]
hsmltomid = hsmltomid[indsn]
distrn = distrn[indsn]
wsph = np.zeros((np.size(indsn)))
indsclose = np.where(distrn <= 0.5 * hsmltomid)
indsfar = np.where((distrn > 0.5 * hsmltomid) & (distrn <= hsmltomid))
distrnclose = distrn[indsclose]
hsmlclose = hsmltomid[indsclose]
wsph[indsclose] = 1 - 6 * (distrnclose / hsmlclose) ** 2 + 6 * (distrnclose / hsmlclose) ** 3
wsph[indsfar] = 2 * (1 - distrn[indsfar] / hsmltomid[indsfar]) ** 3
wsph *= 8 / (np.pi * hsmltomid ** 3)
hm2d[iy, ix] = np.sum(mass[indsnform] * wsph) + np.sum(mass[indspicform])
hm2dtobin[iy, ix] = np.sum(mass[indsnform] * tobin[indsnform] * wsph) + np.sum(mass[indspicform])
hm2dtobinsq[iy, ix] = np.sum(mass[indsnform] * tobin[indsnform] ** 2 * wsph) + np.sum(
mass[indspicform])
hm2dtobin /= hm2d
hm2dtobinsq /= hm2d
hm2dtobin = np.sqrt(hm2dtobinsq - hm2dtobin ** 2)
edx = mh2d.xedges
edy = mh2d.yedges
hm2dtobin *= BaseUnitsDefs.get_conversion_factor_from_string(tobinstr, self.dicunits,
conv_type=quantity_conv_type)
xlab, ylab = [self.get_label_from_string(x, conv_type=coord_conv_type, func_for_label=func) for
x, coord_conv_type, func in
zip([xstr, ystr], [x_conv_type, y_conv_type], [funcofx, funcofy])]
cblab = self.get_label_from_string(tobinstr, func_for_label=funcforpl, label_of_dispersion=True,
conv_type=quantity_conv_type)
if show_plot:
plot_hist2d(hm2dtobin, edy, edx, funcforpl=funcforpl, vmin_plot=vmin_plot, vmax_plot=vmax_plot,
add_eps_in_plot=add_eps_in_plot, cmap=cmap, show_color_bar=show_color_bar, fig=fig,
cb_loc=cb_loc, xlab=xlab, ylab=ylab, cblab=cblab)
return hm2d, edy, edx
[docs] def surface_density_2d(self, compstr, xstr, ystr, inds=None, ids=None, funcofx=identity_function,
funcofy=identity_function, x_conv_type='astro', y_conv_type='astro', xlimit=(-40., 40.),
ylimit=None, nbx=128, nby=None, volstrx='width', volstry='width', surfdens_conv_type='astro',
cic_hist=False, show_plot=True, funcforpl=np.log10, vmin_plot=None, vmax_plot=None,
add_eps_in_plot=False, cmap=None, show_color_bar=True, fig=None, cb_loc=None):
"""Computes and plots 2D surface density (mass weighted hist divided by surface of bins) of component or
particles from indexes.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
xstr, ystr : str
Members of :attr:`~galaximview.snap.Snap.listquantities` . Typically 'xpos', 'ypos' or 'zpos'.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx, funcofy : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function used to create bins on specified dimension. Defaults to no transformation.
x_conv_type, y_conv_type, surfdens_conv_type : {'astro', 'cgs', None}
Type of units conversion for specified axis.
xlimit, ylimit, nbx, nby: optional
See :func:`~galaximview.snap.get_hist2d`
volstrx, volstry : str, optional
Type of cell lengths used for 2D histogram. Set to 'width' (regular cells) if
cic_hist=True. Defaults to 'width' (regular cells).
cic_hist : bool, default = `False`
If `True`, computes cloud in cloud histogram.
show_plot : bool, default = `True`
`True` to show plot.
funcforpl, vmin_plot, vmax_plot, add_eps_in_plot, cmap, show_color_bar, fig, cb_loc: optional
See :func:`~galaximview.funcsplots.plot_hist2d`
Returns
-------
ndarray, ndarray, ndarray
(nby, nbx) Surface density histogram, (nby+1,) y axis bins edges and (nbx,) x axis bins edges.
"""
inds = self.get_indexes_from_user_input(compstr, inds, ids)
xcoords, ycoords = [self.get_array_from_string(x, inds, conv_type=coord_conv_type) for x, coord_conv_type in
zip([xstr, ystr], [x_conv_type, y_conv_type])]
mh2d = get_hist2d(xcoords, ycoords, xlimit, ylimit, nbx, nby, cic_hist=cic_hist, funcofx=funcofx,
funcofy=funcofy)
mass = self.data['mass'][inds]
hm2d, edy, edx = mh2d.volumic_hist2d(weights=mass, volstrx=volstrx, volstry=volstry)
hm2d *= BaseUnitsDefs.get_conversion_factor_from_string('surfdens', self.dicunits, conv_type=surfdens_conv_type)
xlab, ylab = [self.get_label_from_string(x, conv_type=coord_conv_type, func_for_label=func) for
x, coord_conv_type, func in
zip([xstr, ystr], [x_conv_type, y_conv_type], [funcofx, funcofy])]
cblab = self.get_label_from_string('surfdens', func_for_label=funcforpl, conv_type=surfdens_conv_type)
if show_plot:
plot_hist2d(hm2d, edy, edx, funcforpl=funcforpl, vmin_plot=vmin_plot, vmax_plot=vmax_plot,
add_eps_in_plot=add_eps_in_plot, cmap=cmap, show_color_bar=show_color_bar, fig=fig,
cb_loc=cb_loc, xlab=xlab, ylab=ylab, cblab=cblab)
return hm2d, edy, edx
[docs] def binned_quantity_in_2d(self, compstr, xstr, ystr, tobinstr, inds=None, ids=None,
funcofx=identity_function, funcofy=identity_function, x_conv_type='astro',
y_conv_type='astro', xlimit=(-40., 40.), ylimit=None, nbx=128, nby=None,
weight_type='mass', quantity_conv_type='astro', remove_average=False, cic_hist=False,
show_plot=True, funcforpl=np.log10, vmin_plot=None, vmax_plot=None, add_eps_in_plot=False,
cmap=None, show_color_bar=True, fig=None, cb_loc=None):
"""Computes and plot 2D weighted histogram of array for component or particles from indexes or identifiers.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
xstr, ystr : str
Members of :attr:`~galaximview.snap.Snap.listquantities` . Typically 'xpos', 'ypos' or 'zpos'.
tobinstr : str
Member of :attr:`~galaximview.snap.Snap.listquantities`.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx, funcofy : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function used to create bins on specified dimension. Defaults to no transformation.
x_conv_type, y_conv_type, quantity_conv_type : {'astro', 'cgs', None}
Type of units conversion for specified axis.
xlimit, ylimit, nbx, nby: optional
See :func:`~galaximview.snap.get_hist2d`
weight_type : {'mass', 'number'}
Type of weighting for 2D histogramm: 'mass' for mass-weighting, 'number' for number-weighting.
remove_average : bool, default = `False`
If `True`, removes mass-weighted average of histogram.
cic_hist : bool, default = `False`
If `True`, computes cloud in cloud histogram.
show_plot : bool, default = `True`
`True` to show plot.
funcforpl, vmin_plot, vmax_plot, add_eps_in_plot, cmap, show_color_bar, fig, cb_loc: optional
See :func:`~galaximview.funcsplots.plot_hist2d`
Returns
-------
ndarray, ndarray, ndarray
(nby, nbx) Weighted 2D histogram of array, (nby+1,) y axis bins edges and (nbx+1,) x axis bins edges.
"""
inds = self.get_indexes_from_user_input(compstr, inds, ids)
xcoords, ycoords = [self.get_array_from_string(x, inds, conv_type=coord_conv_type) for x, coord_conv_type in
zip([xstr, ystr], [x_conv_type, y_conv_type])]
mass = self.data['mass'][inds]
if (self.listarrs_gas is not None):
if tobinstr in self.listarrs_gas:
indsgas = self.get_indexes_of_component('gas')
inds = np.searchsorted(indsgas, inds)
tobin = self.get_array_from_string(tobinstr, inds, conv_type=None)
mh2d = get_hist2d(xcoords, ycoords, xlimit, ylimit, nbx, nby, cic_hist=cic_hist, funcofx=funcofx,
funcofy=funcofy)
if weight_type == 'mass':
hm2d, edy, edx = mh2d.binned_in_2d_array(tobin, weights=mass)
elif weight_type == 'number':
hm2d, edy, edx = mh2d.binned_in_2d_array(tobin)
else:
raise ValueError('weight_type to be chosen between mass and number.')
if remove_average:
mass_hist, edy, edx = mh2d.hist2d(weights=mass)
inds = np.where(~np.isnan(hm2d))
hist_av = np.average(hm2d[inds], weights=mass_hist[inds])
hm2d -= hist_av
hm2d *= BaseUnitsDefs.get_conversion_factor_from_string(tobinstr, self.dicunits, conv_type=quantity_conv_type)
xlab, ylab = [self.get_label_from_string(x, conv_type=coord_conv_type, func_for_label=func) for
x, coord_conv_type, func in
zip([xstr, ystr], [x_conv_type, y_conv_type], [funcofx, funcofy])]
cblab = self.get_label_from_string(tobinstr, func_for_label=funcforpl, conv_type=quantity_conv_type)
if show_plot:
plot_hist2d(hm2d, edy, edx, funcforpl=funcforpl, vmin_plot=vmin_plot, vmax_plot=vmax_plot,
add_eps_in_plot=add_eps_in_plot, cmap=cmap, show_color_bar=show_color_bar, fig=fig,
cb_loc=cb_loc, xlab=xlab, ylab=ylab, cblab=cblab)
return hm2d, edy, edx
[docs] def quiver(self, compstr, xstr, ystr, vxstr, vystr, inds=None, ids=None, x_conv_type='astro', y_conv_type='astro',
xlimit=(-40., 40.), ylimit=None, nbx=128, nby=None, quantity_conv_type='astro', remove_average=False,
color='k'):
"""Plots quiver plot (arrows), typically of 2D velocities with vxstr='xvel', vystr='yvel', xstr='xpos',
ystr='ypos'.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
xstr, ystr : str
Members of :attr:`~galaximview.snap.Snap.listquantities` . Typically 'xpos', 'ypos' or 'zpos'.
vxstr : str
Member of :attr:`~galaximview.snap.Snap.listquantities`. Typically 'vxpos' if xstr is 'xpos'.
vystr : str
Member of :attr:`~galaximview.snap.Snap.listquantities`. Typically 'vypos' if ystr is 'ypos'.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
x_conv_type, y_conv_type, quantity_conv_type : {'astro', 'cgs', None}
Type of units conversion for specified axis.
xlimit, ylimit, nbx, nby: optional
See :func:`~galaximview.snap.get_hist2d`
remove_average : bool, default = `False`
If `True`, removes mass-weighted average of histogram.
color : matplotlib color, default = 'k'
Matplotlib color of arrows. Defaults to black.
"""
hx, edy, edx = self.binned_quantity_in_2d(compstr, xstr, ystr, xstr, inds=inds, ids=ids,
x_conv_type=x_conv_type, y_conv_type=y_conv_type, xlimit=xlimit,
ylimit=ylimit, nbx=nbx, nby=nby, weight_type='mass',
quantity_conv_type=quantity_conv_type, show_plot=False)
hy, edy, edx = self.binned_quantity_in_2d(compstr, xstr, ystr, ystr, inds=inds, ids=ids,
x_conv_type=x_conv_type, y_conv_type=y_conv_type, xlimit=xlimit,
ylimit=ylimit, nbx=nbx, nby=nby, weight_type='mass',
quantity_conv_type=quantity_conv_type, show_plot=False)
hvx, edy, edx = self.binned_quantity_in_2d(compstr, xstr, ystr, vxstr, inds=inds, ids=ids,
x_conv_type=x_conv_type, y_conv_type=y_conv_type, xlimit=xlimit,
ylimit=ylimit, nbx=nbx, nby=nby, weight_type='mass',
quantity_conv_type=quantity_conv_type, remove_average=remove_average,
show_plot=False)
hvy, edy, edx = self.binned_quantity_in_2d(compstr, xstr, ystr, vystr, inds=inds, ids=ids,
x_conv_type=x_conv_type, y_conv_type=y_conv_type, xlimit=xlimit,
ylimit=ylimit, nbx=nbx, nby=nby, weight_type='mass',
quantity_conv_type=quantity_conv_type, remove_average=remove_average,
show_plot=False)
plt.quiver(hx, hy, hvx, hvy, color=color)
[docs] def binned_dispersion_of_quantity_in_2d(self, compstr, xstr, ystr, tobinstr, inds=None, ids=None,
funcofx=identity_function, funcofy=identity_function, x_conv_type='astro',
y_conv_type='astro', xlimit=(-40., 40.), ylimit=None, nbx=128, nby=None,
weight_type='mass', quantity_conv_type='astro', cic_hist=False,
show_plot=True,
funcforpl=np.log10, vmin_plot=None, vmax_plot=None, add_eps_in_plot=False,
cmap=None, show_color_bar=True, fig=None, cb_loc=None):
"""Computes and plot 2D weighted dispersion of array for component or particles from indexes or identifiers.
Parameters
----------
compstr : str
Member of :attr:`~galaximview.snap.Snap.listcomps`.
xstr, ystr : str
Member of :attr:`~galaximview.snap.Snap.listquantities`. Typically 'xpos', 'ypos' or 'zpos'.
tobinstr : str
Member of :attr:`~galaximview.snap.Snap.listquantities`.
inds, ids : optional
See :meth:`~galaximview.snap.Snap.get_indexes_from_user_input`
funcofx, funcofy : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function used to create bins on specified dimension. Defaults to no transformation.
x_conv_type, y_conv_type, quantity_conv_type : {'astro', 'cgs', None}
Type of units conversion for specified axis.
xlimit, ylimit, nbx, nby: optional
See :func:`~galaximview.snap.get_hist2d`
weight_type : {'mass', 'number'}
Type of weighting for 2D histogramm: 'mass' for mass-weighting, 'number' for number-weighting.
cic_hist : bool, default = `False`
If `True`, computes cloud in cloud histogram.
show_plot : bool, default = `True`
`True` to show plot.
funcforpl, vmin_plot, vmax_plot, add_eps_in_plot, cmap, show_color_bar, fig, cb_loc: optional
See :func:`~galaximview.funcsplots.plot_hist2d`
Returns
-------
ndarray, ndarray, ndarray
(nby, nbx) 2D weighted dispersion histogram of array, (nby+1,) y axis bins edges and (nbx+1,) x axis bins
edges.
"""
inds = self.get_indexes_from_user_input(compstr, inds, ids)
xcoords, ycoords = [self.get_array_from_string(x, inds, conv_type=coord_conv_type) for x, coord_conv_type in
zip([xstr, ystr], [x_conv_type, y_conv_type])]
mh2d = get_hist2d(xcoords, ycoords, xlimit, ylimit, nbx, nby, cic_hist=cic_hist, funcofx=funcofx,
funcofy=funcofy)
mass = self.data['mass'][inds]
if (self.listarrs_gas is not None):
if tobinstr in self.listarrs_gas:
indsgas = self.get_indexes_of_component('gas')
inds = np.searchsorted(indsgas, inds)
tobin = self.get_array_from_string(tobinstr, inds, conv_type=None)
if weight_type == 'mass':
hm2d, edy, edx = mh2d.binned_in_2d_dispersion_of_array(tobin, mass)
elif weight_type == 'number':
hm2d, edy, edx = mh2d.number_weighted_dispersion_2d_hist_of_array(tobin)
else:
raise ValueError('weight_type to be chosen between mass and number.')
hm2d *= BaseUnitsDefs.get_conversion_factor_from_string(tobinstr, self.dicunits, conv_type=quantity_conv_type)
xlab, ylab = [self.get_label_from_string(x, conv_type=coord_conv_type, func_for_label=func) for
x, coord_conv_type, func in
zip([xstr, ystr], [x_conv_type, y_conv_type], [funcofx, funcofy])]
cblab = self.get_label_from_string(tobinstr, func_for_label=funcforpl, label_of_dispersion=True,
conv_type=quantity_conv_type)
if show_plot:
plot_hist2d(hm2d, edy, edx, funcforpl=funcforpl, vmin_plot=vmin_plot, vmax_plot=vmax_plot,
add_eps_in_plot=add_eps_in_plot, cmap=cmap, show_color_bar=show_color_bar, fig=fig,
cb_loc=cb_loc, xlab=xlab, ylab=ylab, cblab=cblab)
return hm2d, edy, edx
[docs]class IterOnSim(SimCaracs):
"""Class for simulation from which snapshots can be obtained.
Parameters
----------
simu_directory : str
the directory containing the snapshots files.
snapshot_basis : str
the basis of the snapshots names, to which formated numbers will be added, see
:meth:`~galaximview.snap.IterOnSim.get_snap`
input_type : {'G1/G2', 'G2_str', 'G4_hdf5'}
Snapshot data type, 'G1/G2' for usual Gadget-2 format.
gadget2_header, ntypes:
Only for Gadget1/Gadget2 snapshots types used in Gadget4. See :func:`~galaximview.readgadget.read_gadget`
centre_on_com : bool, default = `False`
`True` to centre positions and velocities on the centre of mass
store_pos_vel : bool, default = `True`
`True` to store positions and velocities before any rotation, useful for viewer
"""
def __init__(self, simu_directory, snapshot_basis='snapshot_', input_type='G1/G2', gadget2_header=True, ntypes=6,
multiple_snaps=False, centre_on_com=False, store_pos_vel=True):
SimCaracs.__init__(self, simu_directory)
self.direc = simu_directory + '/'
self.input_type = input_type
self.gadget2_header = gadget2_header
self.ntypes = ntypes
self.multiple_snaps = multiple_snaps
self.centre_on_com = centre_on_com
self.snapnamebasis = snapshot_basis
self.storeposvel = store_pos_vel
self.isnapmin, self.isnapmax = self.check_snaps_present()
print('isnapmin, isnapmax = ', self.isnapmin, self.isnapmax)
self.snapt0 = None
[docs] def check_snaps_present(self):
"""Checks that snapshots are numbered as [min_nb, min_nb + 1, min_nb + 2, ..., min_nb + N], if not,
prints that they are not numbered this way. Returns the number of the first snapshot min_nb and the last one.
Returns
-------
int, int
Numbers of the first snapshot and of the last one
"""
snaplist = glob.glob(self.direc + self.snapnamebasis + '*')
if len(snaplist) == 0:
raise ValueError(
'There are no files with the right snapshot name basis in the specified simulation directory.')
isnaplist = []
for snapname in snaplist:
end_str = snapname.split(self.snapnamebasis)[-1]
if '.' in end_str:
isnap = int(end_str.split('.')[0])
else:
isnap = int(end_str)
isnaplist += [isnap]
isnapmin = min(isnaplist)
isnapmax = max(isnaplist)
nb_snaps = len(isnaplist)
if ((isnapmax - isnapmin) != (nb_snaps - 1)):
print('All snapshots are not present?'),
return isnapmin, isnapmax
[docs] def get_snap(self, isnap):
"""Returns object of class Snap.
Parameters
----------
isnap : int
snapshot number
Returns
-------
:class:`~galaximview.snap.Snap` object.
Snapshot object.
"""
fname = self.snapnamebasis + "{0:03d}".format(isnap)
return Snap(self.direc, fname, input_type=self.input_type, gadget2_header=self.gadget2_header,
multiple_snaps=self.multiple_snaps,
ntypes=self.ntypes, in_sim=True, centre_on_com=self.centre_on_com, store_pos_vel=self.storeposvel)
[docs]def get_hist2d(xcoords, ycoords, xrange, yrange, nbx, nby, cic_hist=False, funcofx=identity_function,
funcofy=identity_function):
"""Returns the histogram object (class Hist2DOfXY or CicHist2DOfXY). Sets nby=nbx if nby not given.
Parameters
----------
xcoords : ndarray
(N,) Array to bin for the 'x direction' of the 2D histogram.
ycoords : ndarray
(N,) Array to bin for the 'y direction' of the 2D histogram.
xrange : (float, float)
Limits of x axis.
yrange : (float, float) or (None, None)
Limits of y axis. If (None, None), set as equal to xlimit.
nbx : int
Number of bins for x axis.
nby : int or None
Number of bins for y axis. If is None, computed so that y bins have same length as x bins.
cic_hist : bool, default = `False`
If `True`, computes cloud in cloud histogram.
funcofx : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function used for x axis. Defaults to identity_function (no transfoormation).
funcofy : ufunc, default = :func:`~galaximview.basefuncs.identity_function`
Function used for y axis. Defaults to identity_function (no transfoormation).
Returns
-------
:class:`~galaximview.classeshists.Hist2DOfXY` object.
Hist2DOfXY object.
"""
if yrange is None:
yrange = xrange
if nbx is None:
raise ValueError("nbx should be given/")
if nby is None:
nby = int((yrange[1] - yrange[0]) / (xrange[1] - xrange[0]) * nbx)
if cic_hist:
mh2d = CicHist2DOfXY(xcoords, ycoords, nbx=nbx, nby=nby, xrange=xrange, yrange=yrange)
else:
mh2d = Hist2DOfXY(xcoords, ycoords, xrange, yrange, nbx=nbx, nby=nby, funcofx=funcofx, funcofy=funcofy)
return mh2d