I would like to write a FITS file out of an off-axis projection along some given direction of a field ("i_nu") which I create as a derived field. However, the parameters I set in the main Python program (pl_inu.py) using set_field_parameter() are not passed down to the derived field i_nu. Here is the main macro I am using:
""" pl_inu.py """
import numpy as np import yt from yt.utilities.cosmology import Cosmology from yt import derived_field import matplotlib.pyplot as plt from yt.fields.api import ValidateParameter
from bckfl_fns import b_p,_e_dens, e_dens, _chi_nu,_kappa, f_kappa, _i_nu, f_chi_nu, los_unit_vector, proj_dist
from yt import YTArray, YTQuantity
float_formatter = "{:.3f}".format
cosmo = Cosmology() LFLY = YTArray(1.0/cosmo.hubble_constant.value,"Mpc").in_cgs() TFLY = (2./3./cosmo.hubble_constant).in_cgs() MFLY = YTArray(5.23e12/cosmo.hubble_constant.value,"Msun").in_cgs() uo = {"length_unit":LFLY, "time_unit":TFLY, "mass_unit":MFLY}
fname = 'virgoS12_str_hdf5_plt_cnt_0015' ds=yt.load(fname,units_override=uo) dd=ds.all_data()
ds.add_field(("gas", "n_e"), function=_e_dens, units="cm**(-3)", sampling_type='cell', display_name="Electron density", force_override=True) ds.add_field(("gas", "kappa"), function=_kappa, units="kpc**(-1)", sampling_type='cell', display_name="N(E) norm coeff", force_override=True) ds.add_field(("gas", "chi_nu"), function=_chi_nu, units="kpc**(-1)", sampling_type='cell', display_name="synchr. abs. coeff.", force_override=True) ds.add_field(("gas", "i_nu"), function=_i_nu, units="W/(m**2)/Hz", sampling_type='cell', display_name="Synchr. emission", force_override=True, validators=[ValidateParameter(["lmin","l.o.s.","lmin origin","Background I_nu","Magnetic field","Electron p-law index","b(p)","nu (Ghz)"])])
# Output time (in Myr) # co = Cosmology() time = ds.parameters.get('time')*co.t_from_z(0.).in_units("Myr")
p = 5./2. E_min = 1.5 E_max = 2.e1 bm4 = 0.01
nu = 1.5 dd.set_field_parameter("Magnetic field", bm4) dd.set_field_parameter("Electron p-law index",p) dd.set_field_parameter("1-p", (1-p)) dd.set_field_parameter("b(p)",b_p(p)) dd.set_field_parameter('E_min', E_min) dd.set_field_parameter('E_max', E_max) dd.set_field_parameter("nu (Ghz)", nu)
eul_angles_deg=np.array([4.,22.]) n_dir = los_unit_vector(eul_angles_deg) pos_0 = YTArray([0.,0,-130],'kpc') grid_pos=YTArray([dd['gas','x'],dd['gas','y'],dd['gas','z']]).in_units('kpc') lmin = np.amin(proj_dist(n_dir,grid_pos,pos_0))
i_nu_lmin = YTQuantity(1.e-27,'W/(m**2)/Hz')
dd.set_field_parameter('Background I_nu', i_nu_lmin) dd.set_field_parameter("lmin", lmin) dd.set_field_parameter("l.o.s.", n_dir) dd.set_field_parameter('lmin origin', pos_0) dd.set_field_parameter('grid pos', grid_pos)
synchr_fits=yt.FITSOffAxisSlice(ds,n_dir,("gas", "i_nu"), width=((50,'kpc'))) fits_name = 'i_nu_r0000_15.fits' synchr_fits.writeto(fits_name, overwrite='True') ====
and this is the collection of functions I use to create the field "i_nu": """ bckfl_fns.py """
import sys import yt import numpy as np
import global_vars
from yt import derived_field from yt import YTArray, YTQuantity from yt.fields.api import ValidateParameter
""" Derived fields in yt-FLASH """
def _kappa(field, data): import math from yt import YTQuantity
if data.has_field_parameter("E_min"): E_min = data.get_field_parameter("E_min") else: E_min = 0.
if data.has_field_parameter("E_max"): E_max = data.get_field_parameter("E_max") else: E_max = 0.
ne_tot = data["gas", "n_e"].v
if data.has_field_parameter("Electron p-law index"): p = data.get_field_parameter("Electron p-law index") else: p = 0.
if data.has_field_parameter("E_min") and data.has_field_parameter("E_max") \ and data.has_field_parameter("Electron p-law index"): k=YTQuantity(1.0,"kpc**(-1)")*(1-p)*ne_tot/(pow(E_max,(1-p)) - pow(E_min,(1-p))) print("Inside kappa: p-law index is",p, E_min, E_max) else: k=YTQuantity(1.0,"kpc**(-1)")
if np.any(k < 0): print("k negative!", k)
return k
def _chi_nu(field, data): if data.has_field_parameter("Electron p-law index") and \ data.has_field_parameter("b(p)") and \ data.has_field_parameter("Magnetic field") and \ data.has_field_parameter("nu (Ghz)"): p = data.get_field_parameter("Electron p-law index") bp = data.get_field_parameter("b(p)") bm4 = data.get_field_parameter("Magnetic field") nu = data.get_field_parameter("nu (Ghz)") ne_tot = data[("gas","n_e")].v kappa = data[("gas","kappa")] a = 6.445e-6 chi_nu=a*np.power(10,(p/2))*np.power(5.67,p)*kappa*bp*pow(bm4,(0.5*p+1))/np.power(nu,(0.5*p+2)) else: chi_nu=YTQuantity(0.0,"kpc**(-1)")
return chi_nu
def _i_nu(field, data): from unyt import kpc from yt import YTArray, YTQuantity from yt.units import mass_hydrogen # if data.has_field_parameter("lmin"): lmin = data.get_field_parameter("lmin") print("Inside _i_nu, lmin=",lmin) else: lmin = YTQuantity(0.0,"kpc")
if data.has_field_parameter("l.o.s."): los = data.get_field_parameter("l.o.s.") print("Inside _i_nu, l.o.s.=",los) else: los = [0, 0, 0]
if data.has_field_parameter("lmin origin"): pos_0 = data.get_field_parameter("lmin origin") print("Inside _i_nu, pos_0=",pos_0) else: pos_0 = data.ds.arr(np.zeros(3), 'kpc')
if data.has_field_parameter("Background I_nu"): i_nulmin = data.get_field_parameter("Background I_nu") print("Inside _i_nu, i_nulmin=",i_nulmin) else: i_nu_lmin = 0.
if data.has_field_parameter("Magnetic field"): bm4 = data.get_field_parameter("Magnetic field") print("Inside _i_nu, bm4=",bm4) else: bm4 = 1.e-20
## i_nu_lmin=YTQuantity(i_nulmin,"W/(m**2)/Hz") # i_nu_lmin=YTQuantity(i_nulmin,"W/(m**2)/Hz").in_units("Jansky") if data.has_field_parameter("Electron p-law index"): p = data.get_field_parameter("Electron p-law index") print("Inside _i_nu, p=",p) fact_joc = np.power(10,-p) else: p=0. fact_joc = 0.
if data.has_field_parameter("b(p)"): bp = data.get_field_parameter("b(p)") print("Inside _i_nu, b(p)=",bp) else: bp = 0.
if data.has_field_parameter("nu (Ghz)"): nu = data.get_field_parameter("nu (Ghz)") print("Inside _i_nu, nu=",nu) else: nu = 0
ne_tot = data["gas","density"].v/mass_hydrogen.v kappa = data["kappa"] a = YTQuantity(4.97766e-17,'W/(m**2)/Hz') chi_nu = f_chi_nu(nu,bm4,p,bp,kappa) j_over_chi = a*fact_joc*(nu**2.5)/np.power(bm4,0.5) grid_pos=YTArray([data['gas','x'],data['gas','y'],data['gas','z']]).in_units('kpc') r = proj_dist(los,grid_pos,pos_0) ad = r-lmin print("Inside _i_nu", i_nulmin, j_over_chi, a) b = np.exp(-chi_nu*(r-lmin)) i_nu = i_nulmin*b + j_over_chi*(1-b) return i_nu
""" General purpose functions """
def _e_dens(field, data): """ Computes electron density from gas density.""" from yt.units import mass_hydrogen return data[("gas","density")]/mass_hydrogen
def e_dens(dens): """ Computes electron density from gas density.""" from yt.units import mass_hydrogen return dens.v/mass_hydrogen.v
def f_kappa(ntot,emin,emax,p): import sys, math fk = (1-p)*ntot/(np.power(emax,(1-p))-np.power(emin,(1-p))) return fk
def a_p(p): from scipy import interpolate as int_erp p=[1,1.5,2,2.5,3,3.5,4,4.5,5] bp=[2.056,0.909,0.529,0.359,0.269,0.217,0.186,0.167,0.157] f_interp=int_erp.interp1d(p,bp) return f_interp(p_inp)
def b_p(p_inp): from scipy import interpolate as int_erp p=[1,1.5,2,2.5,3,3.5,4,4.5,5] bp=[0.397,0.314,0.269,0.244,0.233,0.230,0.236,0.248,0.268] f_interp=int_erp.interp1d(p,bp) return f_interp(p_inp)
def f_chi_nu(nu,B,p,bp,kappa): import sys, math a1 = 5.67**p a2 = pow(10,p/2) c1 = 0.5*p + 1 c2 = 0.5*p + 2 abs_coeff = YTArray([6.445e-6],'kpc**(-1)') return abs_coeff
def los_unit_vector(inp_angles): from numpy import radians from math import sin, cos eul_angles_rad=radians(inp_angles) v_dir = np.array([sin(eul_angles_rad[0])*cos(eul_angles_rad[1]),sin(eul_angles_rad[0])*sin(eul_angles_rad[1]),cos(eul_angles_rad[0])]) return v_dir
def proj_dist(los,xyz,origin): import numpy as np d = ((xyz.T - origin)*los) dist = np.sqrt(np.sum(d*d, axis=1)) return dist =======
Fr some reasons I do not quite understand, the field parameters I define in the main program are not passed down to the "i_nu" field through the _i_nu function. Here the error message:
python3 pl_synchr_test.py
yt : [WARNING ] 2021-10-16 17:45:34,946 Overriding code units: Use this option only if you know that the dataset doesn't define the units correctly or at all. yt : [INFO ] 2021-10-16 17:45:34,946 Overriding length_unit: 4.346024761918768e+24 cm. yt : [INFO ] 2021-10-16 17:45:34,947 Overriding time_unit: 2.897349841279178e+17 s. yt : [INFO ] 2021-10-16 17:45:34,947 Overriding mass_unit: 1.4647063306760563e+46 g. yt : [INFO ] 2021-10-16 17:45:35,024 Parameters: current_time = 0.0016112618581139734 yt : [INFO ] 2021-10-16 17:45:35,024 Parameters: domain_dimensions = [8 8 8] yt : [INFO ] 2021-10-16 17:45:35,024 Parameters: domain_left_edge = [-0.7 -0.7 -0.7] yt : [INFO ] 2021-10-16 17:45:35,024 Parameters: domain_right_edge = [0.7 0.7 0.7] yt : [INFO ] 2021-10-16 17:45:35,024 Parameters: cosmological_simulation = 0 Inside _i_nu, lmin= 0.0 Inside _i_nu, l.o.s.= 0.0 Inside _i_nu, pos_0= 0.0 Inside _i_nu, i_nulmin= 0.0 Inside _i_nu, bm4= 0.0 Inside _i_nu, p= 0.0 Inside _i_nu, b(p)= 0.0 Inside _i_nu, nu= 0.0 /g100/home/userexternal/vantonuc/.local/lib/python3.8/site-packages/unyt/array.py:1778: RuntimeWarning: invalid value encountered in true_divide out_arr = func( Inside _i_nu 0.0 nan W/(Hz*m**2) 4.97766e-17 W/(Hz*m**2) yt : [INFO ] 2021-10-16 17:45:41,439 Making a fixed resolution buffer of (('gas', 'i_nu')) 512 by 512 Traceback (most recent call last): File "pl_synchr_test.py", line 64, in <module> synchr_fits=yt.FITSOffAxisSlice(ds,n_dir,("gas", "i_nu"), width=((50,'kpc'))) File "/g100/home/userexternal/vantonuc/.local/lib/python3.8/site-packages/yt/visualization/fits_image.py", line 1275, in __init__ super().__init__(frb, fields=fields, length_unit=lunit, wcs=w) File "/g100/home/userexternal/vantonuc/.local/lib/python3.8/site-packages/yt/visualization/fits_image.py", line 244, in __init__ this_img = img_data[field] File "/g100/home/userexternal/vantonuc/.local/lib/python3.8/site-packages/yt/visualization/fixed_resolution.py", line 139, in __getitem__ buff = self.ds.coordinates.pixelize( File "/g100/home/userexternal/vantonuc/.local/lib/python3.8/site-packages/yt/geometry/coordinates/cartesian_coordinates.py", line 221, in pixelize return self._oblique_pixelize(data_source, field, bounds, size, antialias) File "/g100/home/userexternal/vantonuc/.local/lib/python3.8/site-packages/yt/geometry/coordinates/cartesian_coordinates.py", line 570, in _oblique_pixelize data_source[field], File "/g100/home/userexternal/vantonuc/.local/lib/python3.8/site-packages/yt/data_objects/data_containers.py", line 258, in __getitem__ self.get_data(f) File "/g100/home/userexternal/vantonuc/.local/lib/python3.8/site-packages/yt/data_objects/selection_objects/data_selection_objects.py", line 169, in get_data finfo.check_available(self) File "/g100/home/userexternal/vantonuc/.local/lib/python3.8/site-packages/yt/fields/derived_field.py", line 244, in check_available validator(data) File "/g100/home/userexternal/vantonuc/.local/lib/python3.8/site-packages/yt/fields/derived_field.py", line 503, in __call__ raise NeedsParameter(doesnt_have) yt.fields.field_exceptions.NeedsParameter: (['lmin', 'l.o.s.', 'lmin origin', 'Background I_nu', 'Magnetic field', 'Electron p-law index', 'b(p)', 'nu (Ghz)']) ========
Many thanks in advance for any help/suggestions.