Download
Note: currently this script only works for ESO FLAMES or VIMOS or Gemini GMOS data (but could easily be adapted to work with any other IFU data).
Installation instructions
- Unbundle the tar file in some directory
- In arshift.com there are two links to IRAF libraries. These couple the Fortran code to IRAF, and enable (1) the FITS data and header parameters to be read and written, and (2) access to the iraf way of interpolating (as inimshift, etc), meaning you will be able to use your favourite option for the resampling. You need to find libiraf77.a and ibiminterp.a on your machine and then change the links in arshift.com to point to those files. You may have to consult your systems person but usually the libraries are somewhere accessible.
- Compile the code with >sh arshift.com. This makes arshift.e in the directory.
- In the directory in which you usually start IRAF there is a file called login.cl. Around where you have lines like
#
# Default USER package; extend or modify as you wish. Note that this can
# be used to call FORTRAN programs from IRAF.
#
package user
#
add a line
task diffatmref = /directory/arshft.e
where the directory is the location where you compiled the code.
When you next log in to IRAF you should find a new task called diffatmref and you can enter parameters with epar as any other task.
The method is described in Walsh, J. R., Roy, J.-R., 1990 “Area spectroscopy and correction for differential atmospheric refraction” in Proc. 2nd ESO/ST-ECF Data Analysis Workshop, eds. Baade, D.,Grosbol, P. J., ESO, Garching, p.95.
Help
X and Ypixin should be the spaxel sizes for your IFU data.
X and Ypixout you can set yourself and even make smaller than the input values, but beware. The value of FidWav has to be in the range of the spectrum.
The following header parameters need to be present and filled.
INSTRUME — valid values are 'FLAMES', 'VIMOS', or 'GMOS-S'
CRPIX3 — Start pixel for wavelength array (usually = 1.0)
CRVAL3 — Wavelength (${\rm \AA}$) corresponding to start pixel
CDELT3 — Wavelength increment per pixel (${\rm \AA}$)
DEC — Declination (deg.)
AZIMUTH — Telescope azimuth (deg.) Called ‘HIERARCH ESO TEL AZ’ in ESO headers
AIRMASS
PA_START — Parallactic angle at start of observation Called ‘HIERARCH ESO TEL PARANG START’ in ESO headers
PA_END — Parallactic angle at end of observation (‘HIERARCH ESO TEL PARANG END’ in ESO headers)
ADA_POSANG — Position angle of axis of instrument IFU (‘HIERARCH ESO ADA POSANG’ in ESO headers)
The pressure and temperature are also stored in the ESO hierarchical header under ‘HIERARCH ESO TEL AMBI PRES START’ and ‘HIERARCH ESO TEL AMBI TEMP’. As with things that have start and end values you should take the middle value. How you do this depends on whether you have single images or many combined.
Pressure should be input in units of hPa (1 Torr = 1mmHg, 760 Torr = 1013.25 mbar = 1013.25 hPa)
spline3 is a conservative choice of interpolation function. Don't be afraid to try out different interpolation methods to find which works best for your data.
Example script for running diffatmref (on a Gemini GMOS-S cube):
#first we must copy out the data extension from the cube (MEF file)
#but note problems occur with long filenames -- shorten first
string infile,outfile
real temp,pres
infile = "dcsqtexrgS20080303S0039.fits"
outfile = "pos4_39_DAR.fits"
#first copy out the relevant extension
imdel cube
imcopy(((infile)//"[1]"),"cube")
imgets(((infile)//"[0]"),"TELESCOP")
hedit images="cube.fits" fields="TELESCOP" value=(imgets.value) add+ addonly+ delete- verify+ show+ update+
imgets(((infile)//"[0]"),"INSTRUME")
hedit images="cube.fits" fields="INSTRUME" value=(imgets.value) add+ addonly+ delete- verify+ show+ update+
imgets(((infile)//"[0]"),"DEC")
hedit images="cube.fits" fields="DEC" value=real(imgets.value) add+ addonly+ delete- verify+ show+ update+
imgets(((infile)//"[0]"),"HA")
hedit images="cube.fits" fields="HA" value=(15.*real(imgets.value)) add+ addonly+ delete- verify+ show+ update+
imgets(((infile)//"[0]"),"AZIMUTH")
hedit images="cube.fits" fields="AZIMUTH" value=real(imgets.value) add+ addonly+ delete- verify+ show+ update+
imgets(((infile)//"[0]"),"AIRMASS")
hedit images="cube.fits" fields="AIRMASS" value=real(imgets.value) add+ addonly+ delete- verify+ show+ update+
imgets(((infile)//"[0]"),"PA")
hedit images="cube.fits" fields="ADA_POSANG" value=real(imgets.value) add+ addonly+ delete- verify+ show+ update+
#Get the temp and pressure from the header
imgets(((infile)//"[0]"),"TAMBIENT")
temp=(273+real(imgets.value))
imgets(((infile)//"[0]"),"PRESSUR2")
pres=(real(imgets.value)/100) #pressure measured in hPa (hectoPascals = mbar) = Pa/100
## RUN DAR correction
diffatmref(Incube="cube.fits",Xpixin=0.20,Ypixin=0.20,Xpixout=0.20,Ypixout=0.20,FidWav=6000.0,Xoffset=0.0,Yoffset=0.0,interpol="spline3",Pressure=(pres),Temperature=(temp),Outcube=(outfile))
# if there is a disagreement in the two airmasses printed to screen, must have something to do with how/when they calculate it -- middle or beginning of exposure
# it's fine as long as they agree to within 0.1
Test whether the DAR shifting has worked using wavextract:
# Testing the DAR correction - extract images at the two wavelength extremes
wavextract cube 5100 5150 S cube_5100
wavextract cube 6400 6450 S cube_6400
wavextract cube_DAR 5100 5150 S cube_DAR_5100
wavextract cube_DAR 6400 6450 S cube_DAR_6400
imcntr input="cube_5100.fits" x_init=14.0 y_init=8.0 cboxsize=5
imcntr input="cube_6400.fits" x_init=12.5 y_init=8.0 cboxsize=5
imcntr input="cube_DAR_5100.fits" x_init=13.0 y_init=8.0 cboxsize=5
imcntr input="cube_DAR_6400.fits" x_init=13.0 y_init=8.0 cboxsize=5
#or you could just display them in DS9
Note: x and y offsets specified in the command will offset the resulting cube in the positive direction, i.e. an x-offset of 1 will offset the final cube by +1 in the x direction.
Many thanks to Jeremy Walsh for providing this script and associated help
