# 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