Dear Nathan, Suppose I have this function (given below), where ne, and Temperature is the value corresponding to the grid, but this will be a yt-field, number density, and temperature of the respective grid. def max_jut_dist(p, ne, Temperature) : func = ne.d / (4.0*np.pi * (m_electron.d * speed_of_light.d)**3) theta = boltzmann.d * Temperature.d / (m_electron.d * speed_of_light.d**2) gamma = np.sqrt(1.0 + (p / (m_electron.d * speed_of_light.d))**2) func /= theta func /= syspc.kn(2, 1.0/theta) func *= np.exp( -gamma/theta) return func if I call this function inside the field to compute normalization of this function on each grid over p variable: using ---->>> scipy.integrate.quad( max_jut_dist, g1, g2, args=( data["electron_density"] , data["Temperature"]) ) but this will give an error, shape () doesn't match the broadcast shape (16, 16, 16) I hope this will explain what I want to do. Thanks -Prateek