def limitedfluxdensity(simulation, instrumentname, wavelengths, cmask, filterobject, fwhm, fluxlimit):
# get the path for the data cube corresponding to this instrument
fitspaths = filter(lambda fn: ("_"+instrumentname+"_") in fn, simulation.totalfitspaths())
if len(fitspaths) != 1: return None
# get the data cube and convert it to per-wavelength units
cube = pyfits.getdata(arch.openbinary(fitspaths[0])).T
cube = simulation.convert(cube, to_unit='W/m3/sr', quantity='surfacebrightness', wavelength=wavelengths)
# convolve the data cube to a single frame, and convert back to per-frequency units
frame = filterobject.convolve(wavelengths[cmask], cube[:,:,cmask])
frame = simulation.convert(frame, from_unit='W/m3/sr', to_unit='MJy/sr', wavelength=filterobject.pivotwavelength())
# get information on the simulated pixels (assume all frames are identical and square)
sim_pixels = simulation.instrumentshape()[0]
sim_pixelarea = simulation.angularpixelarea() # in sr
sim_pixelwidth = np.sqrt(sim_pixelarea) * 648000 / np.pi # in arcsec
# convolve the frame with a Gaussian of the appropriate size
frame = gaussian_filter(frame, sigma=fwhm/sim_pixelwidth/2.35482, mode='constant')
# get information on the observational instrument's pixels (assume pixel width ~ fwhm/3)
obs_pixelwidth = fwhm/3 # in arcsec
# calculate the bin size to obtain simulated pixels similar to observed pixels
bin_pixels = find_nearest_divisor(sim_pixels, obs_pixelwidth/sim_pixelwidth)
bin_pixelarea = sim_pixelarea * bin_pixels**2
# rebin the frame
frame = frame.reshape((sim_pixels//bin_pixels,bin_pixels,sim_pixels//bin_pixels,bin_pixels)).mean(axis=3).mean(axis=1)
# integrate over the frame to obtain the total flux density and convert from MJy to Jy
fluxdensity = frame[frame>fluxlimit].sum() * bin_pixelarea * 1e6
return fluxdensity
# This function returns the integer divisor of the first (integer) argument that is nearest to the second (float) argument
评论列表
文章目录
