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.
