# Herschel Data Processing workshop 22 Feb 2012 # # Example script: Calculate a gaussian convolution kernel to compare two # PACS line maps # # Bart Vandenbussche, KU Leuven # Get two projected cubes out of the archive obs1 = getObservation(1342187843, useHsa=True) obs2 = getObservation(1342187845, useHsa=True) cube1 = obs1.refs["level2"].product.refs["HPS3DPB"].product.refs[0].product cube2 = obs2.refs["level2"].product.refs["HPS3DPR"].product.refs[0].product # Integrate the lines and produce integrated line flux maps [map1] = integrateMapFromCube(cube1, nbStack=1, startArray=Double1d([15]), endArray=Double1d([40])) [map2] = integrateMapFromCube(cube2, nbStack=1, startArray=Double1d([5]), endArray=Double1d([15])) map1.setUnit("Jy") map2.setUnit("Jy") lineRatioMap = imageDivide(image1=map1, ref=1, image2=map2) d = Display(lineRatioMap) d.zoomFit() d.setCutLevels(0.5,1.5) # Wavelengths of the two scans: print "Cube 1 Wavelength ", MEAN(cube1.wave) print "Cube 2 Wavelength ",MEAN(cube2.wave) # Download the two beams beam1 = fitsReader(file = '/lhome/bart/pacs/psf/beamcalfilesv2/PCalSpectrometer_Beam_B2A62_mod12_FM_v2.fits') beam2 = fitsReader(file = '/lhome/bart/pacs/psf/beamcalfilesv2/PCalSpectrometer_Beam_R1150_mod12_FM_v2.fits') beamFit1 = sourceFitting(elongated=False,image=beam1, minX=5, minY=5, width=12, height=12) beamFit2 = sourceFitting(elongated=False,image=beam2, minX=5, minY=5, width=12, height=12) sigma1 = beamFit1.sigmaPixels*beam1.wcs.cdelt1 sigma2 = beamFit2.sigmaPixels*beam2.wcs.cdelt1 kernelSigma = SQRT(sigma2**2 - sigma1**2) # Assess the accuracy of the convolution kernel by applying it to beam1 and comparing to beam2: # make a gaussian convolution kernel evaluated for convolution with beam1 kernelModel = herschel.ia.dataset.image.SourceFittingProduct(\ 1,\ beam1.image.dimensions[0]/2.,\ beam1.image.dimensions[1]/2.,\ kernelSigma/beam1.wcs.cdelt1,\ 0.,\ beam1.wcs,\ beam1.unit) kernelImage = kernelModel.evaluate(beam1) convBeam1 = imageConvolution(image=beam1, kernel=kernelImage) # Apply the convolution kernel on map1 kernelModel = herschel.ia.dataset.image.SourceFittingProduct(\ 1,\ map1.image.dimensions[0]/2.,\ map1.image.dimensions[1]/2.,\ kernelSigma/3600./map1.wcs.cdelt1,\ 0.,\ map1.wcs,\ map1.unit) kernelImage = kernelModel.evaluate(map1) convMap1 = imageConvolution(image=map1, kernel=kernelImage) lineRatioMapMatched = imageDivide(image1=convMap1, ref=1, image2=map2) d2 = Display(lineRatioMapMatched) d2.zoomFit() d2.setCutLevels(0.5,1.5) # Line ratios of a staring observation # As an example we take the rebinned cubes of the central pointings of the # rasters above sCube1 = obs1.refs["level2"].product.refs["HPS3DRB"].product.refs[0].product.refs[12].product sCube2 = obs2.refs["level2"].product.refs["HPS3DRR"].product.refs[0].product.refs[4].product [sMap1] = integrateMapFromCube(sCube1, nbStack=1, startArray=Double1d([15]), endArray=Double1d([40])) [sMap2] = integrateMapFromCube(sCube2, nbStack=1, startArray=Double1d([5]), endArray=Double1d([15])) sMap1.setUnit("Jy") sMap2.setUnit("Jy") sLineRatioMap = imageDivide(image1=sMap1, ref=1, image2=sMap2) d1 = Display(sMap1) d1.addLayer(sMap2) d1.zoomFit() d2 = Display(sLineRatioMap) d2.zoomFit() d2.setCutLevels(0.5,1.5) # To be able to compare the line fluxes in the different spaxels: # * Make an interpolated line flux map from the 25 line fluxes measured at # positions sCube1.ra, sCube2.dec # 2D interpolation with irregularly sampled input date not available in hipe. # external packages can be used e.g. using scipy python package: # # from scipy import interpolate # measuredRa = reshape(sCube1.ra) # measuredDec = reshape(sCube1.dec) # measuredFlux = RESHAPE(sMap1) # interpolator = interpolate.interp2d(x, y, z, kind='cubic') # # then interpolate at all ra,dec of every pixel in your interpolated map # # interpolatedFlux = interpolator(ra,dec) # # * Convolve this map with the same kernel as above