"""
pyathena.fields.fields
Registry of derived fields and their plotting metadata.
Each `set_derived_fields_*` function builds and returns six dictionaries:
- func[name] : callable(d, u) -> array-like
- field_dep[name] : set/list of primitive fields required for func[name]
- label[name] : label string (typically LaTeX) for colorbars/axes
- cmap[name] : colormap (string or matplotlib/cmasher colormap)
- vminmax[name] : (vmin, vmax) defaults for visualization
- take_log[name] : bool; if True use LogNorm, else Normalize
"""
import matplotlib as mpl
import numpy as np
import xarray as xr
import astropy.constants as ac
import astropy.units as au
import cmasher as cmr
from matplotlib.colors import Normalize, LogNorm
from ..plt_tools.cmap import cmap_apply_alpha, cmap_shift
from ..microphysics.cool import get_xe_mol
# Optional dependency: xray emissivity helper requires `yt`.
# Keep import soft so users without yt can still use non-Xray functionality.
try:
from .xray_emissivity import get_xray_emissivity
except ModuleNotFoundError:
get_xray_emissivity = None
def static_vars(**kwargs):
"""
Decorator that attaches static attributes to a function.
Used in this file to stash a reference point `x0` on functions that need it.
"""
def decorate(func):
for k in kwargs:
setattr(func, k, kwargs[k])
return func
return decorate
[docs]def set_derived_fields_def(par, x0):
"""
Function to define derived fields info, for example,
functions to calculate derived fields, dependency, label, colormap, etc.
May not work correctly for problems using different unit system.
Assume that density = nH, length unit = pc, etc.
This includes:
- density (rho, nH)
- pressure (pok, pok_trbz)
- distance (r)
- velocity (vmag, vr, vx/vy/vz)
- sound speeds (cs, csound)
- momentum density, ram pressure (Mr, Mr_abs, rhovr2ok)
Parameters
----------
par: dict
Dictionary containing simulation parameter information
x0: sequence of floats
Coordinate of the center with respect to which distance is measured
Returns
-------
Tuple of dictionaries containing derived fields info
"""
func = dict()
field_dep = dict()
label = dict()
cmap = dict()
vminmax = dict()
take_log = dict()
# rho [g cm^-3] - gas density
f = 'rho'
field_dep[f] = ['density']
def _rho(d, u):
return d['density']*(u.muH*u.mH).cgs.value
func[f] = _rho
label[f] = r'$\rho\;[{\rm g\,cm^{-3}}]$'
cmap[f] = 'Spectral_r'
vminmax[f] = (1e-27,1e-20)
take_log[f] = True
# nH [cm^-3] (assume that density = nH)
f = 'nH'
field_dep[f] = ['density']
def _nH(d, u):
return d['density']
func[f] = _nH
label[f] = r'$n_{\rm H}\;[{\rm cm^{-3}}]$'
cmap[f] = 'Spectral_r'
vminmax[f] = (1e-3,1e4)
take_log[f] = True
# P/kB [K cm^-3] - thermal pressure
f = 'pok'
field_dep[f] = set(['pressure'])
def _pok(d, u):
return d['pressure']*(u.energy_density/ac.k_B).cgs.value
func[f] = _pok
label[f] = r'$P\;[k_B\,{\rm cm^{-3}\,K}]$'
cmap[f] = 'inferno'
vminmax[f] = (1e2,1e7)
take_log[f] = True
# rho vz^2/kB [K cm^-3]
f = 'pok_trbz'
field_dep[f] = set(['density','velocity'])
def _pok_trbz(d, u):
return d['density']*d['velocity3']**2*(u.energy_density/ac.k_B).cgs.value
func[f] = _pok_trbz
label[f] = r'$P_{\rm turb,z}\;[k_B\,{\rm cm^{-3}\,K}]$'
cmap[f] = 'inferno'
vminmax[f] = (1e2,1e7)
take_log[f] = True
# Distance from x0 [pc]
f = 'r'
field_dep[f] = set(['density'])
@static_vars(x0=x0)
def _r(d, u):
# NOTE: this assumes d provides coordinates 'x','y','z'
z, y, x = np.meshgrid(d['z'], d['y'], d['x'], indexing='ij')
return xr.DataArray(
np.sqrt((x - _r.x0[0])**2 + (y - _r.x0[1])**2 + (z - _r.x0[2])**2),
dims=('z', 'y', 'x'),
name='r',
)
func[f] = _r
label[f] = r'$r\;[{\rm pc}]$'
cmap[f] = 'viridis'
# Set vmax to the largest box dimension
if "domain1" in par:
blockname = "domain1"
elif "mesh" in par:
blockname = "mesh"
else:
raise KeyError('Expected "domain1" or "mesh" in parameter dictionary.')
Lx = par[blockname]['x1max'] - par[blockname]['x1min']
Ly = par[blockname]['x2max'] - par[blockname]['x2min']
Lz = par[blockname]['x3max'] - par[blockname]['x3min']
Lmax = max(max(Lx, Ly), Lz)
vminmax[f] = (0, Lmax)
take_log[f] = False
# Velocity magnitude [km/s]
f = 'vmag'
field_dep[f] = set(['velocity'])
def _vmag(d, u):
return (d['velocity1']**2 + d['velocity2']**2 + d['velocity3']**2)**0.5*u.kms
func[f] = _vmag
label[f] = r'$|\mathbf{v}|\;[{\rm km\,s^{-1}}]$'
vminmax[f] = (0.1, 1000.0)
cmap[f] = cmr.guppy_r
take_log[f] = True
# Radial velocity w.r.t. x0 [km/s]
f = 'vr'
field_dep[f] = set(['velocity'])
@static_vars(x0=x0)
def _vr(d, u):
z, y, x = np.meshgrid(d['z'], d['y'], d['x'], indexing='ij')
r = xr.DataArray(np.sqrt((x - _vr.x0[0])**2 + (y - _vr.x0[1])**2 + (z - _vr.x0[2])**2),
dims=('z','y','x'), name='r')
return (x*d['velocity1'] + y*d['velocity2'] + z*d['velocity3'])/r*u.kms
func[f] = _vr
label[f] = r'$v_r\;[{\rm km\,s^{-1}}]$'
vminmax[f] = (-20, 100)
# Set cmap midpoint accordingly (midpoint=abs(vmin)/(abs(vmin)+abs(vmax))
cmap[f] = cmap_shift(mpl.cm.BrBG,
midpoint=abs(vminmax[f][0]) / \
(abs(vminmax[f][0]) + abs(vminmax[f][1])),
name='cmap_vr')
take_log[f] = False
# Velocity components [km/s]
def _make_vcomp(idx):
def _vcomp(d, u):
return d[f'velocity{idx}'] * u.kms
return _vcomp
for idx, axis in enumerate(('x', 'y', 'z'), start=1):
f = f'v{axis}'
field_dep[f] = set(['velocity'])
func[f] = _make_vcomp(idx)
label[f] = rf'$v_{axis}\;[{{\rm km\,s^{{-1}}}}]$'
cmap[f] = 'bwr'
vminmax[f] = (-100.0, 100.0)
take_log[f] = False
# -------------------------------------------------------------------------
# Sound speeds [km/s]
# -------------------------------------------------------------------------
# cs [km/s]
# JGK: rename to csound_iso?
f = 'cs'
field_dep[f] = set(['pressure','density'])
def _cs(d, u):
return np.sqrt(d['pressure']/d['density'])
func[f] = _cs
label[f] = r'$c_s\;[{\rm km}\,{\rm s}^{-1}]$'
cmap[f] = 'magma'
vminmax[f] = (0.1,1e3)
take_log[f] = True
# cs [km/s] - sound speed
f = 'csound'
field_dep[f] = set(['pressure','density'])
def _csound(d, u):
return np.sqrt(par['problem']['gamma']*d['pressure']/d['density'])
func[f] = _csound
label[f] = r'$c_s\;[{\rm km}\,{\rm s}^{-1}]$'
cmap[f] = 'magma'
vminmax[f] = (0.1,1e3)
take_log[f] = True
# Radial momentum w.r.t. x0 [km/s cm^-3]
f = 'Mr'
field_dep[f] = set(['density','velocity'])
@static_vars(x0=x0)
def _Mr(d, u):
z, y, x = np.meshgrid(d['z'], d['y'], d['x'], indexing='ij')
r = xr.DataArray(np.sqrt((x - _r.x0[0])**2 + (y - _r.x0[1])**2 + (z - _r.x0[2])**2),
dims=('z','y','x'), name='r')
return d['density']*(x*d['velocity1'] + y*d['velocity2'] + z*d['velocity3'])/r*u.kms
func[f] = _Mr
label[f] = r'$p_r\;[{\rm cm^{-3}\,km\,s^{-1}}]$'
vminmax[f] = (-1e5, 1e5)
# Set cmap midpoint accordingly (midpoint=abs(vmin)/(abs(vmin)+abs(vmax))
cmap[f] = cmap_shift(mpl.cm.BrBG,
midpoint=abs(vminmax[f][0]) / \
(abs(vminmax[f][0]) + abs(vminmax[f][1])),
name='cmap_pyathena_Mr')
take_log[f] = False
# Absolute value of radial momentum w.r.t. x0 [km/s cm^-3]
f = 'Mr_abs'
field_dep[f] = set(['density','velocity'])
@static_vars(x0=x0)
def _Mr_abs(d, u):
z, y, x = np.meshgrid(d['z'], d['y'], d['x'], indexing='ij')
r = xr.DataArray(np.sqrt((x - _r.x0[0])**2 + (y - _r.x0[1])**2 + (z - _r.x0[2])**2),
dims=('z','y','x'), name='r')
return np.abs(d['density']*(x*d['velocity1'] + y*d['velocity2'] + z*d['velocity3'])/r)*u.kms
func[f] = _Mr_abs
label[f] = r'$|p_r|\;[{\rm cm^{-3}\,km\,s^{-1}}]$'
vminmax[f] = (1e-2, 1e4)
# Set cmap midpoint accordingly (midpoint=abs(vmin)/(abs(vmin)+abs(vmax))
cmap[f] = cmap_shift(mpl.cm.BrBG,
midpoint=abs(vminmax[f][0]) / \
(abs(vminmax[f][0]) + abs(vminmax[f][1])),
name='cmap_pyathena_Mr_abs')
take_log[f] = True
# rho vr^2 / kB [cm^-3 K]
f = 'rhovr2ok'
field_dep[f] = set(['density','velocity'])
@static_vars(x0=x0)
def _rhovr2ok(d, u):
z, y, x = np.meshgrid(d['z'], d['y'], d['x'], indexing='ij')
r = xr.DataArray(np.sqrt((x - _r.x0[0])**2 + (y - _r.x0[1])**2 + (z - _r.x0[2])**2),
dims=('z','y','x'), name='r')
return (d['density']*u.density.cgs.value)*\
((x*d['velocity1'] + y*d['velocity2'] + z*d['velocity3'])/r*u.velocity.cgs.value)**2/ac.k_B.cgs.value
func[f] = _rhovr2ok
label[f] = r'$\rho v_r^2\;[{\rm cm}^{-3}\,{\rm K}]$'
vminmax[f] = (1e2, 1e7)
# Set cmap midpoint accordingly (midpoint=abs(vmin)/(abs(vmin)+abs(vmax))
cmap[f] = 'inferno'
take_log[f] = True
return func, field_dep, label, cmap, vminmax, take_log
[docs]def set_derived_fields_cooling(par, newcool):
"""
Function to define derived fields info, for example,
functions to calculate derived fields, dependency, label, colormap, etc.
May not work correctly for problems using different unit system.
Assume that density = nH, length unit = pc, etc.
Parameters
----------
par: dict
Dictionary containing simulation parameter information
newcool: bool
Is new cooling turned on?
Returns
-------
Tuple of dictionaries containing derived fields info
"""
func = dict()
field_dep = dict()
label = dict()
cmap = dict()
vminmax = dict()
take_log = dict()
# T [K] - gas temperature
f = 'T'
if newcool:
field_dep[f] = set(['density','pressure','xe','xH2'])
def _T(d, u):
return d['pressure']/(d['density']*(1.1 + d['xe'] - d['xH2']))/\
(ac.k_B/u.energy_density).cgs.value
else:
field_dep[f] = set(['temperature'])
def _T(d, u):
return d['temperature']
func[f] = _T
label[f] = r'$T\;[{\rm K}]$'
cmap[f] = cmap_shift(mpl.cm.RdYlBu_r, midpoint=3./7., name='cmap_pyathena_T')
vminmax[f] = (1e1,1e7)
take_log[f] = True
# Td [K] - dust temperature
if newcool:
f = 'Td'
field_dep[f] = set(['temperature_dust'])
def _Td(d, u):
return d['temperature_dust']
func[f] = _Td
label[f] = r'$T_{\rm d}\;[{\rm K}]$'
cmap[f] = cmap_shift(mpl.cm.RdYlBu_r, midpoint=3./7., name='cmap_pyathena_Td')
vminmax[f] = (1e0,1e2)
take_log[f] = True
# Cooling rate per volume [erg/s/cm^3] - nH^2*Lambda
f = 'cool_rate_cgs'
field_dep[f] = set(['cool_rate'])
def _cool_rate(d, u):
return d['cool_rate']*(u.energy_density/u.time).cgs.value
func[f] = _cool_rate
label[f] = r'$\mathcal{L}\;[{\rm erg}\,{\rm cm^{-3}}\,{\rm s}^{-1}]$'
cmap[f] = 'cubehelix_r'
vminmax[f] = (1e-30,1e-20)
take_log[f] = True
# Heating rate per volume [erg/s/cm^3] - nH*Gamma
f = 'heat_rate_cgs'
field_dep[f] = set(['heat_rate'])
def _heat_rate(d, u):
return d['heat_rate']*(u.energy_density/u.time).cgs.value
func[f] = _heat_rate
label[f] = r'$\mathcal{G}\;[{\rm erg}\,{\rm cm^{-3}}\,{\rm s}^{-1}]$'
cmap[f] = 'cubehelix_r'
vminmax[f] = (1e-30,1e-20)
take_log[f] = True
# Net cooling rate per volume [erg/s/cm^3] - nH^2*Lambda - nH*Gamma
# (averaged over dt_mhd)
f = 'net_cool_rate'
field_dep[f] = set(['net_cool_rate'])
def _net_cool_rate(d, u):
return d['net_cool_rate']
func[f] = _net_cool_rate
label[f] = r'$\mathcal{L}\;[{\rm erg}\,{\rm cm^{-3}}\,{\rm s}^{-1}]$'
cmap[f] = 'bwr_r'
vminmax[f] = (-1e-20,1e-20)
take_log[f] = False
# Cooling efficiency [erg*cm^3/s]
f = 'Lambda_cool'
field_dep[f] = set(['density','cool_rate'])
def _Lambda_cool(d, u):
return d['cool_rate']*(u.energy_density/u.time).cgs.value/d['density']**2
func[f] = _Lambda_cool
label[f] = r'$\Lambda\;[{\rm erg}\,{\rm cm^{3}}\,{\rm s}^{-1}]$'
cmap[f] = 'cubehelix_r'
vminmax[f] = (1e-30,1e-20)
take_log[f] = True
# Specific cooling rate [erg/s/H]
f = 'nHLambda_cool'
field_dep[f] = set(['density','cool_rate'])
def _nHLambda_cool(d, u):
return d['cool_rate']*(u.energy_density/u.time).cgs.value/d['density']
func[f] = _nHLambda_cool
label[f] = r'$n_{\rm H}\Lambda\;[{\rm erg}\,{\rm cm^{3}}\,{\rm s}^{-1}]$'
cmap[f] = 'cubehelix_r'
vminmax[f] = (1e-30,1e-20)
take_log[f] = True
# Specific net cooling rate [erg/s/H]
f = 'nHLambda_cool_net'
field_dep[f] = set(['density','cool_rate','heat_rate'])
def _nHLambda_cool_net(d, u):
return (d['cool_rate'] - d['heat_rate'])*(u.energy_density/u.time).cgs.value/d['density']
func[f] = _nHLambda_cool_net
label[f] = r'$n_{\rm H}\Lambda_{\rm net}\;[{\rm erg}\,{\rm cm^{3}}\,{\rm s}^{-1}]$'
cmap[f] = 'cubehelix_r'
vminmax[f] = (1e-30,1e-20)
take_log[f] = True
# Heating efficiency [erg/s/H]
f = 'Gamma_heat'
field_dep[f] = set(['density','heat_rate'])
def _Gamma_heat(d, u):
return d['heat_rate']*(u.energy_density/u.time).cgs.value/d['density']
func[f] = _Gamma_heat
label[f] = r'$\Gamma_{\rm heat}\;[{\rm erg}\,{\rm s}^{-1}]$'
cmap[f] = 'cubehelix_r'
vminmax[f] = (1e-30,1e-20)
take_log[f] = True
# Cooling time [Myr]
f = 't_cool'
field_dep[f] = set(['pressure', 'cool_rate'])
def _t_cool(d, u):
return d['pressure']/d['cool_rate']*u.Myr
func[f] = _t_cool
label[f] = r'$t_{\rm cool}\;[{\rm yr}]$'
cmap[f] = 'cubehelix_r'
vminmax[f] = (1e-4,1e2)
take_log[f] = True
return func, field_dep, label, cmap, vminmax, take_log
[docs]def set_derived_fields_mag(par, x0):
"""
Magnetic-field / MHD derived fields.
"""
func = dict()
field_dep = dict()
label = dict()
cmap = dict()
vminmax = dict()
take_log = dict()
# Magnitude of Alfven velocity [km/s]
f = 'vAmag'
field_dep[f] = set(['density', 'cell_centered_B'])
def _vAmag(d, u):
return (d['cell_centered_B1']**2 +
d['cell_centered_B2']**2 +
d['cell_centered_B3']**2)**0.5/np.sqrt(d['density'])*u.kms
func[f] = _vAmag
label[f] = r'$v_A\;[{\rm km}\,{\rm s}^{-1}]$'
vminmax[f] = (0.1, 1000.0)
cmap[f] = 'cividis'
take_log[f] = True
def _make_vA_component(idx):
"""Return Alfven velocity component vA{idx} in km/s."""
def _vA(d, u):
return d[f'cell_centered_B{idx}']*np.sqrt(d['density'])*u.kms
return _vA
def _make_B_uG_component(idx):
"""Return magnetic field component B{idx} in micro-Gauss."""
def _B(d, u):
return d[f'cell_centered_B{idx}']*np.sqrt(u.energy_density.cgs.value)*np.sqrt(4.0*np.pi)*1e6
return _B
for idx, axis in enumerate(('x', 'y', 'z'), start=1):
# vA components
f = f'vA{axis}'
field_dep[f] = set(['density', 'cell_centered_B'])
func[f] = _make_vA_component(idx)
label[f] = rf'$v_{{A,{axis}}}\;[{{\rm km/s}}]$'
cmap[f] = 'RdBu'
vminmax[f] = (-1e2, 1e2)
take_log[f] = False
# B components
f = f'B{axis}'
field_dep[f] = set(['cell_centered_B'])
func[f] = _make_B_uG_component(idx)
label[f] = rf'$B_{{{axis}}}\;[\mu{{\rm G}}]$'
cmap[f] = 'RdBu'
vminmax[f] = (-1e2, 1e2)
take_log[f] = False
# Magnetic fields magnitude [micro-Gauss]
f = 'Bmag'
field_dep[f] = set(['cell_centered_B'])
def _Bmag(d, u):
return (d['cell_centered_B1']**2 +
d['cell_centered_B2']**2 +
d['cell_centered_B3']**2)**0.5*np.sqrt(u.energy_density.cgs.value)\
*np.sqrt(4.0*np.pi)*1e6
func[f] = _Bmag
label[f] = r'$|\mathbf{B}|\;[\mu{\rm G}]$'
vminmax[f] = (1e-1, 1e2)
cmap[f] = 'cividis'
take_log[f] = True
# Magnetic pressure magnitude [Pmag/kB]
f = 'pok_mag'
field_dep[f] = set(['cell_centered_B'])
def _pok_mag(d, u):
return (d['cell_centered_B1']**2 +
d['cell_centered_B2']**2 +
d['cell_centered_B3']**2)**0.5*(u.energy_density/ac.k_B).cgs.value
func[f] = _pok_mag
label[f] = r'$P_{\rm mag}\;[k_B\,{\rm cm^{-3}\,K}]$'
cmap[f] = 'inferno'
vminmax[f] = (1e2,1e7)
take_log[f] = True
return func, field_dep, label, cmap, vminmax, take_log
[docs]def set_derived_fields_rad(par, x0):
"""Register radiation field derived fields.
Includes normalized FUV (PE and LW band), C ionizing, and H2
dissociating radiation field strengths in Draine ISRF units
(chi_PE, chi_LW, chi_FUV, chi_CI, chi_H2).
Parameters
----------
par : dict
Dictionary containing simulation parameter information.
x0 : sequence of floats
Reference coordinate (unused here, kept for interface consistency).
Returns
-------
tuple
Six dictionaries: (func, field_dep, label, cmap, vminmax, take_log).
"""
func = dict()
field_dep = dict()
label = dict()
cmap = dict()
vminmax = dict()
take_log = dict()
# Dust PE opacity for Z'=1
kappa_dust_PE_def = 418.7
try:
Erad_PE0 = par['cooling']['Erad_PE0']
Erad_LW0 = par['cooling']['Erad_LW0']
except KeyError:
Erad_PE0 = 7.613e-14
Erad_LW0 = 1.335e-14
# Normalized PE radiation field strength (Draine field unit)
f = 'chi_PE'
field_dep[f] = set(['rad_energy_density_PE'])
def _chi_PE(d, u):
return d['rad_energy_density_PE']*(u.energy_density.cgs.value/Erad_PE0)
func[f] = _chi_PE
label[f] = r'$\chi_{\rm PE}$'
cmap[f] = 'viridis'
vminmax[f] = (1e-4,1e4)
take_log[f] = True
# Normalized LW radiation field strength (Draine field unit)
f = 'chi_LW'
field_dep[f] = set(['rad_energy_density_LW'])
def _chi_LW(d, u):
return d['rad_energy_density_LW']*(u.energy_density.cgs.value/Erad_LW0)
func[f] = _chi_LW
label[f] = r'$\chi_{\rm LW}$'
cmap[f] = 'viridis'
vminmax[f] = (1e-4,1e4)
take_log[f] = True
# Normalized FUV radiation field strength (Draine field unit)
f = 'chi_FUV'
field_dep[f] = set(['rad_energy_density_PE','rad_energy_density_LW'])
def _chi_FUV(d, u):
return (d['rad_energy_density_PE'] + d['rad_energy_density_LW'])*\
(u.energy_density.cgs.value/(Erad_PE0 + Erad_LW0))
func[f] = _chi_FUV
label[f] = r'$\chi_{\rm FUV}$'
cmap[f] = 'viridis'
vminmax[f] = (1e-4,1e4)
take_log[f] = True
# Normalized LW radiation field
f = 'chi_CI'
field_dep[f] = set(['rad_energy_density_CI'])
def _chi_CI(d, u):
return d['rad_energy_density_CI']*(u.energy_density.cgs.value/Erad_LW0)
func[f] = _chi_CI
label[f] = r'$\chi_{\rm C^0}$'
cmap[f] = 'viridis'
vminmax[f] = (1e-8,1e4)
take_log[f] = True
# Normalized LW radiation field strength (Draine ISRF)
f = 'chi_H2'
field_dep[f] = set(['rad_energy_density_LW_diss'])
def _chi_H2(d, u):
return d['rad_energy_density_LW_diss']*(u.energy_density.cgs.value/Erad_LW0)
func[f] = _chi_H2
label[f] = r'$\chi_{\rm H_2}$'
cmap[f] = 'viridis'
vminmax[f] = (1e-8,1e4)
take_log[f] = True
return func, field_dep, label, cmap, vminmax, take_log
[docs]def set_derived_fields_newcool(par, x0):
"""Register new-cooling (NCR) chemistry derived fields.
Includes temperature, molecular/atomic hydrogen fractions and densities
(T, xH2, nH2, xHI, nHI, xe), carbon/oxygen species (xCII, nCII, xCI,
nCI, xCO, nCO, xOI), cosmic-ray ionization rate, and cooling/heating rates.
Parameters
----------
par : dict
Dictionary containing simulation parameter information.
x0 : sequence of floats
Reference coordinate (unused here, kept for interface consistency).
Returns
-------
tuple
Six dictionaries: (func, field_dep, label, cmap, vminmax, take_log).
"""
func = dict()
field_dep = dict()
label = dict()
cmap = dict()
vminmax = dict()
take_log = dict()
# set metallicity and standard abundances
try:
Zgas = par['problem']['Z_gas']
Zdust = par['problem']['Z_dust']
except KeyError:
Zgas = 1.0
Zdust = 1.0
try:
xCstd = par['cooling']['xCstd']
except KeyError:
xCstd = 1.6e-4
print('xCstd not found. Use {:.2g}.'.format(xCstd))
xCtot = Zgas*xCstd
try:
xOstd = Zgas*par['cooling']['xOstd']
except KeyError:
xOstd = 3.2e-4
print('xOstd not found. Use {:.2g}.'.format(xOstd))
xOtot = Zgas*xOstd
# set xi_cr0
try:
if "photchem_ncr" in par:
xi_CR0= par["photchem_ncr"]["xi_cr0"]
else:
xi_CR0 = par['problem']['xi_CR0']
except KeyError:
xi_CR0 = 2.e-16
print('xi_CR0 not found. Use {:.2g}.'.format(xi_CR0))
# temperature
f = "T"
field_dep[f] = set(['density','pressure','xe','xH2'])
def _T(d, u):
return d['pressure']/(d['density']*(1.1 + d['xe'] - d['xH2']))/\
(ac.k_B/u.energy_density).cgs.value
func[f] = _T
label[f] = r'$T\;[{\rm K}]$'
cmap[f] = cmap_shift(mpl.cm.RdYlBu_r, midpoint=3./7., name='cmap_pyathena_T')
vminmax[f] = (1e1,1e7)
take_log[f] = True
# nH2 [cm^-3] (assume d=nH)
f = 'nH2'
field_dep[f] = set(['density', 'xH2'])
def _nH2(d, u):
return d['density']*d['xH2']
func[f] = _nH2
label[f] = r'$n_{\rm H_2}\;[{\rm cm^{-3}}]$'
cmap[f] = cmap_apply_alpha('Greens')
vminmax[f] = (1e0,1e4)
take_log[f] = True
# 2nH2 [cm^-3] (assume d=nH)
f = '2nH2'
field_dep[f] = set(['density', 'xH2'])
def _2nH2(d, u):
return 2.0*d['density']*d['xH2']
func[f] = _2nH2
label[f] = r'$2n_{\rm H_2}\;[{\rm cm^{-3}}]$'
cmap[f] = cmap_apply_alpha('Greens')
vminmax[f] = (1e0,1e4)
take_log[f] = True
# xH2 [cm^-3] (assume d=nH)
f = 'xH2'
field_dep[f] = set(['xH2'])
def _xH2(d, u):
return d['xH2']
func[f] = _xH2
label[f] = r'$x_{\rm H_2}$'
cmap[f] = 'viridis'
vminmax[f] = (0,0.5)
take_log[f] = False
# 2xH2 [cm^-3] (assume d=nH)
f = '2xH2'
field_dep[f] = set(['xH2'])
def _2xH2(d, u):
return 2*d['xH2']
func[f] = _2xH2
label[f] = r'$2x_{\rm H_2}$'
cmap[f] = 'viridis'
vminmax[f] = (0,1.0)
take_log[f] = False
# nHI [cm^-3]
f = 'nHI'
field_dep[f] = set(['density', 'xHI'])
def _nHI(d, u):
return d['density']*d['xHI']
func[f] = _nHI
label[f] = r'$n_{\rm H^0}\;[{\rm cm^{-3}}]$'
cmap[f] = cmap_apply_alpha('Blues')
vminmax[f] = (1e-3,1e4)
take_log[f] = True
# xHI
f = 'xHI'
field_dep[f] = set(['xHI'])
def _xHI(d, u):
return d['xHI']
func[f] = _xHI
label[f] = r'$x_{\rm H^0}$'
cmap[f] = 'viridis'
vminmax[f] = (0.0, 1.0)
take_log[f] = False
# nHII [cm^-3]
f = 'nHII'
field_dep[f] = set(['density', 'xHI', 'xH2'])
def _nHII(d, u):
return d['density']*(1.0 - d['xHI'] - 2.0*d['xH2'])
func[f] = _nHII
label[f] = r'$n_{\rm H^+}\;[{\rm cm^{-3}}]$'
cmap[f] = cmap_apply_alpha('Oranges')
vminmax[f] = (1e-3,1e4)
take_log[f] = True
# xHII
f = 'xHII'
field_dep[f] = set(['xH2','xHI'])
def _xHII(d, u):
return 1.0 - d['xHI'] - 2.0*d['xH2']
func[f] = _xHII
label[f] = r'$x_{\rm H^+}$'
cmap[f] = 'viridis'
vminmax[f] = (0.0, 1.0)
take_log[f] = False
# nHn [cm^-3]
f = 'nHn'
field_dep[f] = set(['density', 'xHI', 'xH2'])
def _nHn(d, u):
return d['density']*(d['xHI'] + 2.0*d['xH2'])
func[f] = _nHn
label[f] = r'$n_{\rm H^0} + 2n_{\rm H_2}\;[{\rm cm^{-3}}]$'
cmap[f] = cmap_apply_alpha('Blues')
vminmax[f] = (1e-3,1e4)
take_log[f] = True
# xn [cm^-3]
f = 'xn'
field_dep[f] = set(['xHI', 'xH2'])
def _xn(d, u):
return d['xHI'] + 2.0*d['xH2']
func[f] = _xn
label[f] = r'$x_{\rm n}$'
cmap[f] = cmap_apply_alpha('YlGn')
vminmax[f] = (0,1)
take_log[f] = False
# ne
f = 'ne'
field_dep[f] = set(['density','xe'])
def _ne(d, u):
return d['density']*d['xe']
func[f] = _ne
label[f] = r'$n_{\rm e}$'
cmap[f] = 'viridis'
vminmax[f] = (1e-4, 1e3)
take_log[f] = True
# nesq
f = 'nesq'
field_dep[f] = set(['density','xe'])
def _nesq(d, u):
return (d['density']*d['xe'])**2
func[f] = _nesq
label[f] = r'$n_{\rm e}^2$'
cmap[f] = 'plasma'
vminmax[f] = (1e-8, 1e6)
take_log[f] = True
# xe
f = 'xe'
field_dep[f] = set(['xe'])
def _xe(d, u):
return d['xe']
func[f] = _xe
label[f] = r'$x_{\rm e}$'
cmap[f] = 'viridis'
vminmax[f] = (1e-6, 1.208)
take_log[f] = True
# xCI - atomic neutral carbon
f = 'xCI'
field_dep[f] = set(['xCI_over_xCtot'])
def _xCI(d, u):
# Apply floor and ceiling
return np.maximum(0.0,np.minimum(xCtot,d['xCI_over_xCtot']*xCtot))
func[f] = _xCI
label[f] = r'$x_{\rm CI}$'
cmap[f] = 'viridis'
vminmax[f] = (0,xCtot)
take_log[f] = False
# nCI
f = 'nCI'
field_dep[f] = set(['density','xCI_over_xCtot'])
def _nCI(d, u):
return d['density']*xCtot*np.maximum(0.0,np.minimum(1.0,d['xCI_over_xCtot']))
func[f] = _nCI
label[f] = r'$x_{\rm CI}$'
cmap[f] = 'viridis'
vminmax[f] = (1e2*xCtot,1e4*xCtot)
take_log[f] = True
# xOII - single ionized oxygen
f = 'xOII'
field_dep[f] = set(['xH2','xHI'])
def _xOII(d, u):
return xOtot*(1.0 - d['xHI'] - 2.0*d['xH2'])
func[f] = _xOII
label[f] = r'$x_{\rm OII}$'
cmap[f] = 'viridis'
vminmax[f] = (0,xOtot)
take_log[f] = False
# xCII - singly ionized carbon
# Use with caution.
# (Do not apply to hot gas and depend on cooling implementation)
f = 'xCII'
field_dep[f] = set(['xe','xH2','xHI','pressure','density','CR_ionization_rate'])
def _xCII(d, u):
d['T'] = d['pressure']/(d['density']*(1.1 + d['xe'] - d['xH2']))/\
(ac.k_B/u.energy_density).cgs.value
xe_mol = get_xe_mol(d['density'],d['xH2'],d['xe'],d['T'],d['CR_ionization_rate'],
Zgas, Zdust)
# Apply floor and ceiling
return np.maximum(0.0,np.minimum(xCtot,
d['xe'] - (1.0 - d['xHI'] - 2.0*d['xH2'])*(1.0 + xOtot) - xe_mol))
func[f] = _xCII
label[f] = r'$x_{\rm CII}$'
cmap[f] = 'viridis'
vminmax[f] = (0,xCtot)
take_log[f] = False
# xCII - singly ionized carbon (use when CR ionization is uniform everywhere)
# Use with caution.
# (Do not apply to hot gas and depend on cooling implementation)
f = 'xCII_alt'
field_dep[f] = set(['xe','xH2','xHI','pressure','density'])
def _xCII_alt(d, u):
d['T'] = d['pressure']/(d['density']*(1.1 + d['xe'] - d['xH2']))/\
(ac.k_B/u.energy_density).cgs.value
xe_mol = get_xe_mol(d['density'],d['xH2'],d['xe'],d['T'],xi_CR0,
Zgas, Zdust)
# Apply floor and ceiling
return np.maximum(0.0,np.minimum(xCtot,
d['xe'] - (1.0 - d['xHI'] - 2.0*d['xH2'])*(1.0 + xOtot) - xe_mol))
func[f] = _xCII_alt
label[f] = r'$x_{\rm CII}$'
cmap[f] = 'viridis'
vminmax[f] = (0,xCtot)
take_log[f] = False
# xi_CR
f = 'xi_CR'
field_dep[f] = set(['CR_ionization_rate'])
def _xi_CR(d, u):
return d['CR_ionization_rate']
func[f] = _xi_CR
label[f] = r'$\xi_{\rm CR}\;[{\rm s}^{-1}]$'
cmap[f] = 'viridis'
vminmax[f] = (1e-17,1e-15)
take_log[f] = True
# T_alt [K]
f = 'T_alt'
field_dep[f] = set(['pressure','density','xe','xH2'])
def _T_alt(d, u):
return d['pressure']*(u.energy_density/ac.k_B).cgs.value/(d['density']*(1.1 + d['xe'] - d['xH2']))
func[f] = _T_alt
label[f] = r'$T_{\rm alt}\;[{\rm K}]$'
cmap[f] = cmap_shift(mpl.cm.RdYlBu_r, midpoint=3./7., name='cmap_pyathena_T')
vminmax[f] = (1e1,1e7)
take_log[f] = True
return func, field_dep, label, cmap, vminmax, take_log
[docs]def set_derived_fields_sixray(par, x0):
"""Register six-ray external radiation field derived fields.
Includes external (attenuated) PE and LW radiation field strengths in
Draine ISRF units (chi_PE_ext, chi_LW_ext, chi_H2_ext).
Parameters
----------
par : dict
Dictionary containing simulation parameter information.
x0 : sequence of floats
Reference coordinate (unused here, kept for interface consistency).
Returns
-------
tuple
Six dictionaries: (func, field_dep, label, cmap, vminmax, take_log).
"""
func = dict()
field_dep = dict()
label = dict()
cmap = dict()
vminmax = dict()
take_log = dict()
try:
Erad_PE0 = par['cooling']['Erad_PE0']
Erad_LW0 = par['cooling']['Erad_LW0']
except KeyError:
Erad_PE0 = 7.613e-14
Erad_LW0 = 1.335e-14
# Normalized FUV radiation field strength (Draine ISRF)
f = 'chi_PE_ext'
field_dep[f] = set(['rad_energy_density_PE_ext'])
def _chi_PE_ext(d, u):
return d['rad_energy_density_PE_ext']*(u.energy_density.cgs.value/Erad_PE0)
func[f] = _chi_PE_ext
label[f] = r'$\chi_{\rm PE,ext}$'
cmap[f] = 'viridis'
vminmax[f] = (1e-4,1e4)
take_log[f] = True
# Normalized LW radiation field
f = 'chi_LW_ext'
field_dep[f] = set(['rad_energy_density_LW_ext'])
def _chi_LW_ext(d, u):
return d['rad_energy_density_LW_ext']*(u.energy_density.cgs.value/Erad_LW0)
func[f] = _chi_LW_ext
label[f] = r'$\chi_{\rm LW,ext}$'
cmap[f] = 'viridis'
vminmax[f] = (1e-4,1e4)
take_log[f] = True
# Normalized FUV radiation field strength (Draine field unit)
f = 'chi_FUV_ext'
field_dep[f] = set(['rad_energy_density_PE_ext','rad_energy_density_LW_ext'])
def _chi_FUV_ext(d, u):
return (d['rad_energy_density_PE_ext'] + d['rad_energy_density_LW_ext'])*\
(u.energy_density.cgs.value/(Erad_PE0 + Erad_LW0))
func[f] = _chi_FUV_ext
label[f] = r'$\chi_{\rm FUV,ext}$'
cmap[f] = 'viridis'
vminmax[f] = (1e-4,1e4)
take_log[f] = True
# Normalized LW radiation field
f = 'chi_CI_ext'
field_dep[f] = set(['rad_energy_density_CI_ext'])
def _chi_CI_ext(d, u):
return d['rad_energy_density_CI_ext']*(u.energy_density.cgs.value/Erad_LW0)
func[f] = _chi_CI_ext
label[f] = r'$\chi_{\rm CI,ext}$'
cmap[f] = 'viridis'
vminmax[f] = (1e-4,1e4)
take_log[f] = True
# Normalized LW radiation field strength (Draine ISRF)
f = 'chi_H2_ext'
field_dep[f] = set(['rad_energy_density_LW_diss_ext'])
def _chi_H2_ext(d, u):
return d['rad_energy_density_LW_diss_ext']*(u.energy_density.cgs.value/Erad_LW0)
func[f] = _chi_H2_ext
label[f] = r'$\chi_{\rm H2,ext}$'
cmap[f] = 'viridis'
vminmax[f] = (1e-8,1e2)
take_log[f] = True
return func, field_dep, label, cmap, vminmax, take_log
[docs]def set_derived_fields_xray(par, x0, newcool):
"""Register X-ray emissivity derived fields.
Registers the X-ray emissivity field (j_X) in the 0.5–7 keV band,
computed using APEC emissivity tables via the ``yt`` library.
.. note::
Requires the optional dependency ``yt``. Raises
:exc:`ModuleNotFoundError` if ``yt`` is not installed.
Parameters
----------
par : dict
Dictionary containing simulation parameter information.
x0 : sequence of floats
Reference coordinate (unused here, kept for interface consistency).
newcool : bool
If ``True``, use new-cooling chemistry fields for the emissivity
calculation; otherwise use classic fields.
Returns
-------
tuple
Six dictionaries: (func, field_dep, label, cmap, vminmax, take_log).
"""
# Fail fast with a clear message if yt-dependent helper is missing.
if get_xray_emissivity is None:
raise ModuleNotFoundError(
"X-ray derived fields require the optional dependency `yt`.\n"
"Install it with:\n"
" python -m pip install yt\n"
"or:\n"
" conda install -c conda-forge yt"
)
func = dict()
field_dep = dict()
label = dict()
cmap = dict()
vminmax = dict()
take_log = dict()
# TODO-JKIM: Need Z_gas parameter in the problem block
# Or metallicity field
# set metallicity and standard abundances
try:
Z_gas = par['problem']['Z_gas']
except KeyError:
Z_gas = 1.0
emin = 0.5 # keV
emax = 7.0 # keV
energy = True # If set to False, returns photon emissivity [#/s/cm^3]
# Normalized FUV radiation field strength (Draine field unit)
f = 'j_X'
if newcool:
field_dep[f] = set(['density','pressure','xe','xH2'])
else:
field_dep[f] = set(['density','temperature'])
# Frequency integrated volume emissivity
def _j_Xray(d, u):
if newcool:
d['temperature'] = d['pressure']/(d['density']*(1.1 + d['xe'] - d['xH2']))/\
(ac.k_B/u.energy_density).cgs.value
em = get_xray_emissivity(d['temperature'].data, Z_gas,
emin, emax, energy=energy)
return d['density']**2*em
func[f] = _j_Xray
label[f] = r'$j_{\rm X}$ (0.5-7.0 keV) $[{\rm erg}\,{\rm s}^{-1}\,{\rm cm}^{-3}]$'
cmap[f] = 'plasma'
vminmax[f] = (1e-34,1e-25)
take_log[f] = True
return func, field_dep, label, cmap, vminmax, take_log
[docs]def set_derived_fields_cosmic_ray(par):
"""Register cosmic-ray derived fields.
Includes CR parallel diffusion coefficient (sigma_para), CR pressure
(pok_cr), CR energy density (Ec), and CR-related velocity fields.
Parameters
----------
par : dict
Dictionary containing simulation parameter information. Must include
a ``cr`` block with at least ``vmax`` (maximum CR streaming speed
in cm/s).
Returns
-------
tuple
Six dictionaries: (func, field_dep, label, cmap, vminmax, take_log).
"""
func = dict()
field_dep = dict()
label = dict()
cmap = dict()
vminmax = dict()
take_log = dict()
gamma_cr = 4./3.
vmax = par["cr"]["vmax"] # in cm/s
# cr diffusion coefficient
f = 'sigma_para'
field_dep[f] = set(['0-Sigma_diff1'])
def _sigma_para(d, u):
units = (u.length**2/u.time).cgs.value
vmax_ = vmax/(u.cm/u.s)
return d['0-Sigma_diff1']/vmax_/units
func[f] = _sigma_para
label[f] = r'$\sigma_\parallel\;[{\rm cm^{-2}\,s}]$'
cmap[f] = "viridis"
vminmax[f] = [6e-30,1.5e-27]
take_log[f] = True
# cr pressure
f = 'pok_cr'
field_dep[f] = set(['0-Ec'])
def _pok_cr(d, u):
return d['0-Ec']*(gamma_cr-1)*(u.energy_density/ac.k_B).cgs.value
func[f] = _pok_cr
label[f] = r'$P_{\rm cr}\;[k_B\,{\rm cm^{-3}\,K}]$'
cmap[f] = 'YlOrRd_r'
vminmax[f] = (1,5e4)
take_log[f] = True
# cr fluxes
f = 'Fcr1'
field_dep[f] = set(['0-Fc1'])
def _Fcr1(d, u):
vmax_ = vmax/(u.cm/u.s)
return d['0-Fc1']*vmax_*u.erg/u.kpc**2/(1.e6*u.Myr)
func[f] = _Fcr1
label[f] = r'$F_{\rm cr,x}\;[{\rm erg\,kpc^{-2}\,yr^{-1}}]$'
cmap[f] = 'pink'
vminmax[f] = (5e43,2.5e46)
take_log[f] = True
f = 'Fcr2'
field_dep[f] = set(['0-Fc2'])
def _Fcr2(d, u):
vmax_ = vmax/(u.cm/u.s)
return d['0-Fc2']*vmax_*u.erg/u.kpc**2/(1.e6*u.Myr)
func[f] = _Fcr2
label[f] = r'$F_{\rm cr,y}\;[{\rm erg\,kpc^{-2}\,yr^{-1}}]$'
cmap[f] = 'pink'
vminmax[f] = (5e43,2.5e46)
take_log[f] = True
f = 'Fcr3'
field_dep[f] = set(['0-Fc3'])
def _Fcr3(d, u):
vmax_ = vmax/(u.cm/u.s)
return d['0-Fc3']*vmax_*u.erg/u.kpc**2/(1.e6*u.Myr)
func[f] = _Fcr3
label[f] = r'$F_{\rm cr,z}\;[{\rm erg\,kpc^{-2}\,yr^{-1}}]$'
cmap[f] = 'pink'
vminmax[f] = (5e43,2.5e46)
take_log[f] = True
f = 'Fcr_mag'
field_dep[f] = set(['0-Fc1','0-Fc2','0-Fc3'])
def _Fcr_mag(d, u):
vmax_ = vmax/(u.cm/u.s)
Fmag = np.sqrt(d['0-Fc1']**2+d['0-Fc2']**2+d['0-Fc3']**2)
return Fmag*vmax_*u.erg/u.kpc**2/(1.e6*u.Myr)
func[f] = _Fcr_mag
label[f] = r'$F_{\rm cr}\;[{\rm erg\,kpc^{-2}\,yr^{-1}}]$'
cmap[f] = 'pink'
vminmax[f] = (5e43,2.5e46)
take_log[f] = True
# Ion Alfven speed
f = 'VAi_mag'
field_dep[f] = set(['0-Vs1','0-Vs2','0-Vs3'])
def _VAi_mag(d, u):
vmag = np.sqrt(d['0-Vs1']**2+d['0-Vs2']**2+d['0-Vs3']**2)
return vmag*u.kms
func[f] = _VAi_mag
label[f] = r'$|v_{\rm A,i}|\;[{\rm km\,s^{-1}}]$'
cmap[f] = 'managua_r'
vminmax[f] = (2,250)
take_log[f] = True
f = 'VAi1'
field_dep[f] = set(['0-Vs1','0-Vs2','0-Vs3'])
def _VAi1(d, u):
return d["0-Vs1"]*u.kms
func[f] = _VAi1
label[f] = r'$v_{\rm A,i,x}\;[{\rm km\,s^{-1}}]$'
cmap[f] = 'bwr'
vminmax[f] = (-100.0,100.0)
take_log[f] = False
f = 'VAi2'
field_dep[f] = set(['0-Vs1','0-Vs2','0-Vs3'])
def _VAi2(d, u):
return d["0-Vs2"]*u.kms
func[f] = _VAi2
label[f] = r'$v_{\rm A,i,y}\;[{\rm km\,s^{-1}}]$'
cmap[f] = 'bwr'
vminmax[f] = (-100.0,100.0)
take_log[f] = False
f = 'VAi3'
field_dep[f] = set(['0-Vs1','0-Vs2','0-Vs3'])
def _VAi3(d, u):
return d["0-Vs3"]*u.kms
func[f] = _VAi3
label[f] = r'$v_{\rm A,i,z}\;[{\rm km\,s^{-1}}]$'
cmap[f] = 'bwr'
vminmax[f] = (-100.0,100.0)
take_log[f] = False
# Vcr
f = 'Vcr_mag'
field_dep[f] = set(['0-Fc1','0-Fc2','0-Fc3','0-Ec'])
def _Vcr(d, u):
vmax_ = vmax/(u.cm/u.s)
Fmag = np.sqrt(d['0-Fc1']**2+d['0-Fc2']**2+d['0-Fc3']**2)
return Fmag*vmax_/(4/3.*d['0-Ec'])*u.kms
func[f] = _Vcr
label[f] = r'$|v_{\rm eff}|\;[{\rm km\,s^{-1}}]$'
cmap[f] = 'cool'
vminmax[f] = (2,250)
take_log[f] = True
# CR scalar
f = "rCR"
field_dep[f] = set(["rCR"])
def _rCR(d, u):
return d["rCR"]
func[f] = _rCR
label[f] = r'$r_{\rm CR}$'
cmap[f] = cmr.ghostlight
vminmax[f] = (0, 1)
take_log[f] = False
# CR injection
f = "pok_cr_inj"
field_dep[f] = set(["eCRinj"])
def _pok_CR_inj(d, u):
return d["eCRinj"]*(gamma_cr-1)*u.pok
func[f] = _pok_CR_inj
label[f] = r'$P_{\rm CR,inj}\;[k_B\,{\rm cm^{-3}\,K}]$'
cmap[f] = 'YlOrRd_r'
vminmax[f] = (1,5e4)
take_log[f] = True
return func, field_dep, label, cmap, vminmax, take_log
[docs]def set_derived_fields_feedback_scalars(par):
"""Register feedback scalar derived fields.
Includes gas metallicity (Zgas) and SN ejecta mass fraction (rSN)
derived from passive scalar fields tracked during stellar feedback.
Parameters
----------
par : dict
Dictionary containing simulation parameter information.
Returns
-------
tuple
Six dictionaries: (func, field_dep, label, cmap, vminmax, take_log).
"""
func = dict()
field_dep = dict()
label = dict()
cmap = dict()
vminmax = dict()
take_log = dict()
Zsolar=0.02 # need to get from parameter
# metallicity
f = 'Zgas'
field_dep[f] = set(['rmetal'])
def _Zgas(d, u):
return d['rmetal']/Zsolar
func[f] = _Zgas
label[f] = r'$Z/Z_\odot$'
cmap[f] = cmr.pepper_r
vminmax[f] = [1,10]
take_log[f] = True
# SN ejecta mass fraction
f = 'rSN'
field_dep[f] = set(['rSN'])
def _rSN(d, u):
return d['rSN']
func[f] = _rSN
label[f] = r'$f_{\rm SN}$'
cmap[f] = cmr.bubblegum
vminmax[f] = [0,1]
take_log[f] = False
# mass return fraction
f = 'rret'
field_dep[f] = set(['rret'])
def _rret(d, u):
return d['rret']
func[f] = _rret
label[f] = r'$f_{\rm ret}$'
cmap[f] = cmr.bubblegum
vminmax[f] = [0,1]
take_log[f] = False
return func, field_dep, label, cmap, vminmax, take_log
[docs]class DerivedFields(object):
"""Registry of derived fields for an Athena/Athena++ simulation.
Reads simulation parameters and builds dictionaries of derived field
metadata (compute function, primitive field dependencies, plot label,
colormap, value range, and log-scale flag). The active set of fields
is determined automatically from the ``par`` configuration (e.g. MHD,
cooling, radiation, cosmic rays).
Parameters
----------
par : dict
Simulation parameter dictionary as returned by
:func:`pyathena.read_athinput`.
x0 : array-like of float, optional
Reference position ``[x, y, z]`` used for distance-based fields
such as radial velocity. Default is ``[0.0, 0.0, 0.0]``.
Attributes
----------
func : dict
Maps field name → callable ``f(d, u)`` that computes the field.
field_dep : dict
Maps field name → set of required primitive field names.
label : dict
Maps field name → LaTeX label string for colorbars/axes.
cmap : dict
Maps field name → colormap.
vminmax : dict
Maps field name → ``(vmin, vmax)`` defaults for visualization.
take_log : dict
Maps field name → bool; if ``True`` use log normalization.
norm : dict
Maps field name → :class:`matplotlib.colors.Normalize` instance.
imshow_args : dict
Maps field name → dict of kwargs suitable for imshow/colorbar calls.
dfi : dict
Combined per-field metadata dict (union of all the above).
Examples
--------
>>> from pyathena import read_athinput
>>> from pyathena.fields.fields import DerivedFields
>>> par = read_athinput("/path/to/athinput")
>>> df = DerivedFields(par)
>>> print(list(df.func.keys())[:5])
"""
def __init__(self, par, x0=np.array([0.0,0.0,0.0])):
# Create a dictionary containing all information about derived fields
self.dfi = dict()
athenapp = "mesh" in par
# default configuration
cooling = False
newcool = False
mhd = False
radps = False
sixray = False
wind = False
xray = False
cosmic_ray = False
feedback_scalars = False
if athenapp:
# Athena++ configuration
try:
mhd = par["configure"]["Magnetic_fields"] == "ON"
except KeyError:
pass
try:
newcool = par["photchem"]["mode"] == "ncr"
except KeyError:
pass
try:
cooling = par["cooling"]["cooling"] != "none"
except KeyError:
pass
try:
cosmic_ray = par["configure"]["Cosmic_Ray_Transport"] == "Multigroups"
except KeyError:
pass
try:
feedback_scalars = par["configure"]["Number_of_feedback_scalar"] > 0
except KeyError:
pass
else:
# Athena configuration
try:
newcool = par['configure']['new_cooling'] == 'ON'
except KeyError:
pass
try:
mhd = par['configure']['gas'] == 'mhd'
except KeyError:
pass
try:
cooling = par['configure']['cooling'] == 'ON'
except KeyError:
pass
try:
radps = par['configure']['radps'] == 'ON'
except KeyError:
pass
try:
sixray = par['configure']['sixray'] == 'ON'
except KeyError:
pass
try:
wind = par['feedback']['iWind'] != 0
except KeyError:
pass
try:
xray = ((par['feedback']['iSN'] > 0) or
(par['feedback']['iWind'] > 0) or
(par['feedback']['iEarly'] > 0))
except KeyError:
pass
self.func, self.field_dep, \
self.label, self.cmap, \
self.vminmax, self.take_log = set_derived_fields_def(par, x0)
dicts = (self.func, self.field_dep, self.label, self.cmap, \
self.vminmax, self.take_log)
if mhd:
dicts_ = set_derived_fields_mag(par, x0)
for d, d_ in zip(dicts, dicts_):
d.update(d_)
if cooling or newcool:
dicts_ = set_derived_fields_cooling(par, newcool)
for d, d_ in zip(dicts, dicts_):
d.update(d_)
if newcool:
dicts_ = set_derived_fields_newcool(par, x0)
for d, d_ in zip(dicts, dicts_):
d.update(d_)
if radps:
dicts_ = set_derived_fields_rad(par, x0)
for d, d_ in zip(dicts, dicts_):
d.update(d_)
if sixray:
dicts_ = set_derived_fields_sixray(par, x0)
for d, d_ in zip(dicts, dicts_):
d.update(d_)
# if wind:
# dicts_ = set_derived_fields_wind(par, x0)
# for d, d_ in zip(dicts, dicts_):
# d.update(d_)
if xray:
dicts_ = set_derived_fields_xray(par, x0, newcool)
for d, d_ in zip(dicts, dicts_):
d.update(d_)
# TIGRESS++ specific
if cosmic_ray:
dicts_ = set_derived_fields_cosmic_ray(par)
for d, d_ in zip(dicts, dicts_):
d.update(d_)
if feedback_scalars:
dicts_ = set_derived_fields_feedback_scalars(par)
for d, d_ in zip(dicts, dicts_):
d.update(d_)
self.derived_field_list = self.func
# Set colormap normalization and scale
self.norm = dict()
self.scale = dict()
self.imshow_args = dict()
for f in self.derived_field_list:
if self.take_log[f]:
self.norm[f] = LogNorm(*self.vminmax[f])
self.scale[f] = 'log'
else:
self.norm[f] = Normalize(*self.vminmax[f])
self.scale[f] = 'linear'
self.imshow_args[f] = dict(norm=self.norm[f], cmap=self.cmap[f],
cbar_kwargs=dict(label=self.label[f]))
for f in self.derived_field_list:
self.dfi[f] = dict(field_dep=self.field_dep[f],
func=self.func[f],
label=self.label[f],
norm=self.norm[f],
vminmax=self.vminmax[f],
cmap=self.cmap[f],
scale=self.scale[f],
take_log=self.take_log[f],
imshow_args=self.imshow_args[f])
