Subaru Data Reduction Cookbook:
Grism Spectroscopic Observations with IRCS
— Version. 2.1.3e (January 5, 2010) —Based on the textbook in Japanese written by M. Imanishi
for the Subaru Data Reduction School held in May, 2008 (Version 2.0)
Current Editor of English Version: R. S. Furuya,
together with the combined eﬀort of the past and current staﬀ at Subaru Telescope11

Foreward

This Cookbook presents the typical procedure used to analyze infrared spectra. It is our expectation that a novice at grism spectroscopic observations using Infrared Camera and Spectrograph(IRCS) at Subaru Telescope will be able to produce a ﬁnal infrared spectrum by him/herself. The
focus of this Cookbook is to describe how to analyze low-wavelength resolution infrared spectra
having a moderate signal-to-noise (S/N) ratio. The targets of infrared low-resolution spectroscopy,
like Subaru IRCS grism spectroscopic observations, are generally faint where the S/N ratio withrespect to the continuum may be at most 10–30. On the other hand, high-dispersion spectroscopic
observations are often carried out on bright objects, with an aim of achieving very high S/N ratios
exceeding e.g., ∼100. Analysis of such high dispersion and high S/N ratio spectra requires a few
more procedures where special care is needed in addition to the procedures described in the Cookbook. As for the data handling of such spectra, we comment brieﬂy on such additional procedures.
Last but not least, we appreciate any feedback from the readers to help improve this Cookbook.

2008 October 30
Masatoshi IMANISHI
author of the original version in Japanese
E-mail: masa.imanishi "at" nao.ac.jp
Ray S. FURUYA
current editor of English version
E-mail: rsf "at" subaru.naoj.org

2

Revision History
Version
2.1.0eDate2008 October 30

2.1.1e2.1.2e2008 November 30
2008 December 10

2.1.3e2010 January 5First release based on the Japanese version [v.2.0 (2008 April 10)]by R.S.F. & M.I.
Language correction for the Subaru Asia Winter School by R.S.F. & A.H.
Minor corrections after the Winter School & changed the directory structureof the sample data by R.S.F.
Minor cosmetic corrections by R.S.F.

3

Contents
1 Foreward

2

2 Introduction to Spectroscopic Observations in the Near-Infrared

5

3 IRCS

6

4 A Standard Procedure of Infrared Spectroscopic Observation

7

5 A Standard Procedure of Infrared Spectroscopic Data

9

5.1

First fact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

5.2

Data Reduction: the Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

5.2.1

Getting Started: Setting up Your IRAF . . . . . . . . . . . . . . . . . . . . . . . . .

10

5.2.2

The Sample Data and Initial Data Inspection . . . . . . . . . . . . . . . . . . . . . .

12

5.2.3

Subtraction: A – B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

5.2.4

Dividing with Spectroscopic Flat . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

5.2.5

Eliminating Bad Pixels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

5.2.6

Extracting Signals and Combining Data taken at Diﬀerent Slit Positions . . . . . . .

26

5.2.7

Correction for Non-linearity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

5.2.8

Wavelength Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

5.2.9

Airmass Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

5.2.10 Binning Spectrum Channels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

5.2.11 Dividing the Object Spectrum by Standard Star Spectrum . . . . . . . . . . . . . .

42

5.2.12 Multiplying Standard Stars Blackbody Spectrum . . . . . . . . . . . . . . . . . . . .

43

5.2.13 Flux Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

A Error Estimates

48

4

2

Introduction to Spectroscopic Observations in the NearInfrared

The basic concepts of near-infrared (NIR) spectroscopic observations are essentially the same as
those at optical wavelengths. However, there are a few diﬀerences that characterize NIR spectroscopy, as described below:
1. Since CCDs used in the optical observations do not have sensitivity at the wavelength longer
than 1µm, detectors with mercury cadmium telluride (HgCdTe) or indium antimonide (InSb)
are generally used.
2. Compared to the optical wavelengths, thermal emission from the Earth’s atmosphere andtelescope itself is much higher (Fig. 1), making observations (e.g., calibration) diﬃcult. Thisrequires additional data handling processes1 . In addition, the resultant signal-to-noise ratio(S/N) is degraded by Poisson noises from background emission. Furthermore, since infrared
detectors cannot count numbers of photons exceeding a certain value (i.e., the full-well level,or simply full-well), one has to acquire data with a short integration time to prevent detector
saturation. Acquiring data with such a short exposure time results in low observing eﬃciency,
and produces a signiﬁcant number of data frames (which occupies a large amount of hard disk
space).
3. In general, calibration accuracy produced by infrared data reduction is lower than that for the
optical regime. However, in the case of low-resolution infrared spectroscopy, our experiencesuggests that the quality of the ﬁnal spectra is determined by the background noise rather
than the calibration uncertainties.
Astronomical signal can be observed at all the optical wavelengths, whereas there exist someinfrared wavelength ranges where the Earth’s atmosphere is not transparent, as shown in Figure
2. Therefore, ground-based infrared observations are limited to the wavelength ranges where the1

The major noise sources associated with ground-based observations can be categorized into the following threesources: (1) signal from the astronomical object(s) itself becomes a major noise source, i.e., signal-limited. This eﬀect
is obvious especially in the high-dispersion observations toward very bright sourcesin the optical wavelength, (2) readout noise of the detector(s) generally becomes the major noise source in the case of the high-dispersion observations at
infrared (readout-limited), and (3) Poisson noise of the background emission is generally the dominant noise source
for the case of the low-dispersion observations in the infrared (background-limited). The degree of contributions
from these noise sources depends on various parameters such as, e.g., source brightness, instrumental performance,
wavelength resolution and/or atmospheric conditions. Therefore, we suggest that the user should always pay attention
to which noise source most aﬀects the quality of the data in order to select the best observing parameters.

5

Figure 1: A typical spectrum of the background emission between the optical and the infrared
wavelengths where one ﬁnds that the noise level will increase signiﬁcantly λ & 2µm. Here, AE
denotes thermal radiation from Earth’s atmosphere (220 – 273 K), GBT : thermal emission fromtelescope, OH : OH airglow, i.e., OH vibration-rotation bands, ZSL : scattered light from the Sun
due to dust in zodiacal plane, ZE : thermal emission from such dust. The plot is from Astronomy
Lecture Notebook (in Japanese) by Prof. Iwamuro at Kyoto University (http://www.kusastro.kyotou.ac.jp/∼ iwamuro/LECTURE/OBS/).

atmospheric transparency is at a reasonable level (ideally close to unity). Such a wavelength range
is often referred to as a band.

3

IRCS

IRCS (Infrared Camera and Spectrograph) at Subaru Telescope oﬀers both imaging and spectroscopic observational capability in the infrared. For further information, please visit the IRCS web
page, http://www.naoj.org/Observing/Instruments/IRCS/index.html. Spectroscopic observations using IRCS can be categorized into low-dispersion (i.e., low wavelength-resolution) single-orderobservations and the high resolution cross-dispersion (such as the HDS in the case of the Subaru
6Figure 2: A plot of Earth’s atmospheric transmission, i.e, infrared wavelength window(http://www.jach.hawaii.edu/UKIRT/astronomy/utils/atmos-index.html).
Ground-based
infrared observations are limited to the (wavelength) bands where the Earth’s atmosphere is transparent. These bands are referred to as J, H, K, L, M, N, and Q. The Earth’s atmospheric opacity
is highly time variable with several conditions. In particular, it is known to be sensitive to the total
column density of water vapor above the telescope site. The illustrated transmission curve is thegood conditions case (i.e., the water vapor level is low). As water vapor level rises, the transmission
level degrades at all wavelength. Such an eﬀect is especially prominent at wavelength whose transmission is intrinsically low. See Figure 11 for a magniﬁed view of the L band transmission curve.

Telescope). However, for the novice, we will focus on describing on how to reduce the former. If you
are interested in reducing data from cross-dispersion spectroscopy, we suggest reading ”A User’s
Guide to Reducing Echelle Spectra with IRAF” (Daryl Willmarth and Jeanette Barnes, 1994 May;
http://iraf.noao.edu/docs/spectra.html).

4

A Standard Procedure of Infrared Spectroscopic Observation

In this section, we describe a standard procedure for data reduction of infrared spectroscopic observations. For simplicity’s sake, we limit our discussion only to observations of a point-like source.

7

Atmosphere

Figure 3: A sketch illustrating a concept of airmass. Airmass is deﬁned by sec(z) = 1/ cos(z) wherez is the zenith angle. As is clear from the deﬁnition, it takes a minimum (= 1.0) at zenith. The
degree of the attenuation of the signal is proportional to the distance in the atmosphere that thephoton has traveled, i.e., proportional to the airmass. One should plan to observe target sources at
an airmass as small as possible: less than ∼ 2.0 (image courtesy: Prof. Iwamuro as the same URL
for Figure 1).
1. Observing an object at two positions (hereafter named A and B) along a slit — Conﬁgure nondestructive readout (NDR)2 , of the detector and the numbers of exposure at each position(COADD) accordingly. Keep in mind that each exposure should be completed within a reasonable timescale comparable to (or less than) the characteristic time scale of the Earth’s
atmospheric variation.
2. Taking several sets of data with an ABBA dithering sequence
3. Observing bright standard star (a main-sequence G star in many cases) before and after observing
the science target(s) — We suggest observing such a standard star at an airmass diﬀerence(Fig.3) of less than 0.1 with respect to the target object. Recall that atmospheric transmission
in the infrared will vary signiﬁcantly with weather conditions and observing wavelength(s).
One can correct for the variation of the Earth’s atmospheric transmission through dividing
the target spectrum by that of the standard star whose intrinsic spectrum (i.e., the spectrum
free from the absorption by the Earth’s atmosphere) is well-known.
Given the above observing procedure, we will describe a typical procedure to reduce lowdispersion infrared spectroscopic data.
2

Contrary to the CCDs used in the optical, infrared detectors can read out the accumulated electric charge manytimes without destruction (Non Destructive Readout; NDR), allowing us to reduce read-out noise. The maximum
number of NDRs is given by the read-out speed and the exposure time. The default setting at IRCS is conﬁgured
to maximize the NDR numbers. In the case of L-band (λ = 2.8 − 4.2 µm) low-resolution spectroscopy such as the
sample data used as an example in this Cookbook, Poisson noise from the background would overcome the read-out
noise. Therefore, instead of using the default setting at IRCS, you have to wisely optimize the numbers of NDRs,
because unnecessarily large number of NDRs reduce the observing eﬃciency.

85

A Standard Procedure of Infrared Spectroscopic Data

5.1

First fact

Listed below are the typical steps required for the analysis of infrared spectroscopic data.
1. Subtracting data taken at a slit position from those at another, i.e., A-B
2. Correcting for bad pixels
3. Dividing by a spectroscopic ﬂat frame
4. Extracting the signal, then combining A and B
5. Correction for non-linearity (we will skip this procedure for the sample data.)
6. Wavelength calibration
7. Airmass correction (we will skip this procedure for the sample data.)
8. Binning the obtained spectra (if necessary)
9. Dividing the object spectrum by a standard star spectrum
10. Multiplying the intrinsic spectrum of the standard star
11. Flux calibration
12. Obtaining the ﬁnal spectrum
Due to diﬃculty in the long-wavelength infrared observations described above, the number ofsuccessful observations is clearly smaller than those at the short wavelength. The smaller numberof observation yields less feedbacks from the data reduction process as compared to shorter wavelengths. We, however, believe that it will not be so diﬃcult to learn data reduction of the shorter
wavelength observations once you have familiarized yourself with relatively long wavelength data.
In this Cookbook, we will use the already published L-band (λ = 2.8 − 4.2 µm) spectroscopic
data toward the ultra-luminous infrared galaxy IRAS 00188−0856 (Imanishi et al. 2006 ApJ 637
114).

9

5.2

Data Reduction: the Details

We will utilize the widely used multi-purpose data reduction package — the Image Reduction andAnalysis Facility (IRAF) distributed and maintained by the National Optical Astronomy Observatory (NOAO) of the United States. See http://iraf.noao.edu/ for getting further information.
At http://iraf.noao.edu/docs/spectra.html, you will ﬁnd a helpful guide book ”A Beginner’s
Guide to Using IRAF (IRAF Version 2.10) by Jeanette Barnes (1993 August).
If you are looking for a document focusing on single-order spectroscopy, we suggest reading ”AUser’s Guide to Reducing Slit Spectra with IRAF” by Phil Massey, Frank Valdes, and Jeanette
Barnes (1992 April; http://iraf.noao.edu/docs/spectra.html).
In this Cookbook, we describe a method to reduce spectroscopic data on the basis of therecipes already available in these documents. We believe that the basic concept of data handling isessentially the same as those used in many references, although some optimizations will be requiredin order to reduce the sample data.

5.2.1

Getting Started: Setting up Your IRAF

We assume that you have a Linux/Unix machine where IRAF is available. If not, please consult
with your local experts. The PC-IRAF web page
http://iraf.noao.edu/projects/pc-iraf/pc-iraf.html would also help you when you havetrouble 3 . To start IRAF, go to any directory where you will work with IRAF, e.g., ”ircs sp”
in the case of the sample data described in §5.2.2, type ”mkiraf”. You may have to answer somequestions with y and xgterm.

ana03 /home/imanshms/ircs_sp> mkirafInitialize uparm? (y|n): y-- initializing uparm
Terminal types: xgterm,xterm,gterm,vt640,vt100,etc.
Enter terminal type: xgterm
A new LOGIN.CL file has been created in the current directory.
You may wish to review and edit this file to change the defaults.

3

If some IRAF commands do not work on your computer, your version of IRAF may not be compatible withthose commands.

10This will create an ASCII text ﬁle named login.cl which deﬁnes your IRAF environment.
Open the ﬁle with your favorite editor (e.g., vi, emacs, etc) and ﬁnd the line:#set

imtype= "imh"We suggest modifying this line to (don’t forget to remove # at the beginning),
set

imtype= "fits"

This tells IRAF that your default image format is FITS rather than the IRAF one.
Type ”xgterm &” from any Linux/Unix command line, then type cl in the xgterm window.
You will get an IRAF environment with the following message.

NOAO Sun/IRAF Revision 2.12.1-EXPORT Wed Jul 10 13:33:04 MST 2002
This is the EXPORT version of Sun/IRAF V2.12 for SunOS 4 and Solaris 2.8
Welcome to IRAF. To list the available commands, type ? or ??. To get
detailed information about a command, type ‘help command’.
To run acommand or load a package, type its name.
Type ‘bye’ to exit a
package, or ‘logout’ to get out of the CL.
Type ‘news’ to find out
what is new in the version of the system you are using.
The following
commands or packages are currently defined:
apropos
color.
dataio.
dbms.

guiapps.
guidemo.
images.
language.

lists.
mscred.
noao.
obsolete.

plot.
proto.
softools.
stsdas.

system.
tables.
utilities.

cl>Here cl> is the IRAF prompt. To work on your data, you also need to start an image viewer
ds9 (or saoimage). Type ”ds9 &” to start it.

11

5.2.2The Sample Data and Initial Data Inspection

We assume that you were able to download the archive of sample data as well as some scripts andprograms, ircs sp.tar.gz, fromhttp://optik2.mtk.nao.ac.jp/subaru red/IRCS/files/ircs sp.tar.gz. To extract the thetar archive, typetar zxvf ircs sp.tar.gz
You have a directory named ircs sp which has the four sub-directories as shown below. Although name of the directory is arbitrary, we suggest that the users follow our nomenclature usedin this Cookbook. Both the ana and spec directories are empty at this stage because they are thedirectories where we will store the FITS and ASCII text format data created in the data reduction
steps.

ircs_spana/(scratch directory for analysis: fits data)

fits/(the raw data in fits format)

spec/(scratch directory for analysis: ASCII data)

tools/(scripts and programs)

You will ﬁnd a series of data ﬁles in FITS format under the fits directory. If you are
curious to see more information, pick one of the FITS ﬁles to see its header by typing ”more
IRCA00070351.fits”, you will see the basic parameters of the data such as object name, observing
date, observing mode, and integration time.

OBJECT= ’HR72

DATE-OBS= ’2002-08-19’ /’ / UT date of Observation (yyyy-mm-dd)

Let’s go back to the working directory, ”ircs sp”, by typing ”cd ..” (from the fits/ directory) to see the raw data. In the case of sample data, sequential FITS ﬁles, IRCA00070351 –IRCA00070362, are the data of a standard star. To inspect the raw data, we use the task display.
Optimize several input parameters to show the image accordingly. See below.
12cl> epar display
I R A F
Image Reduction and Analysis Facility
PACKAGE = tv
TASK = display
image
=
frame
=
(bpmask =(bpdispl=(bpcolor=(overlay=(ocolors=(erase =(border_=(select_=(repeat =(fill=(zscale =(contras=(zrange =(zmask =(nsample=(xcenter=(ycenter=(xsize =(ysize =
(xmag
=
(ymag
=
(order =
(z1
=
(z2
=
(ztrans =(lutfile=

IRCA00070351.fits image to be displayed
1 frame to be written intoBPM) bad pixel masknone) bad pixel display (none|overlay|interpolate)red) bad pixel colors) overlay mask
green) overlay colorsyes) erase frameno) erase unfilled area of windowyes) display frame being loadedno) repeat previous display parametersno) scale image to fit display windowno) display range of greylevels near median0.25) contrast adjustment for zscale algorithmno) display full image intensity range) sample mask
1000) maximum number of sample pixels to use
0.5) display window horizontal center0.5) display window vertical center
1.) display window horizontal size
1.) display window vertical size0.5) display window horizontal magnification
0.5) display window vertical magnification0) spatial interpolator order (0=replicate, 1=linea) minimum greylevel to be displayed) maximum greylevel to be displayed
linear) greylevel transformation (linear|log|none|user)) file containing user defined look up table13

(mode=ql)

Here, we have modiﬁed the default parameters to xmag = 0.5, ymag = 0.5, zscale = no,and zrange = no. The ﬁrst two parameters display the whole area of the image using 1024×1024
pixels. The third and fourth parameters manually control the display range of the data, i.e., theintensity scale. You can leave all the remaining default parameters as they are. To modify theparameters, bring your cursor to the desired line and put the parameter value, and type ”:q” toquit the input interface. If you want to change (or erase) a variable, just hit the SPACE key onthat parameter. To initialize all the parameters to the default values, issue ”unlearn display” at
IRAF command line (i.e., cl> unlearn display).
After setting all the parameters, try the command,cl> display fits/IRCA00070351.fits 1 z1=0 z2=30000

You will see the standard star data (the very bright one) taken at the slit position A as shown in
the left panel of Figure 4 where the dispersed signal is seen in the vertical direction. On the other
hand, it would be diﬃcult to identify the signal from the targets as they are generally faint. In theﬁgure, you can recognize that the data consists of three parts split in the horizontal direction; they
represent data taken with the diﬀerent slit widths (each section corresponds to a slit width of 0.′′ 9,0.′′ 6, and 0.′′ 3 from the left). It should be noted that the data having the wider slit width have thehigher background signal level as the total power of the incident radiation is proportional to the slit
width. In the case of the sample data, the standard star can be recognized in the left hand portionof the data, i.e., on the 0.′′ 9 width slit.
In general, slit width should be selected on the basis of the stellar image size (i.e., seeing) duringobservations. Considering scientiﬁc goals as well, you may have to select an appropriate slit width
to achieve the highest S/N ratios. If you select a narrow slit width, you reduce the backgroundnoise level, but also you may lose signals under poor seeing conditions (because of slit loss). Sucha degradation may also occur if the tracking accuracy of telescope is poor. Needless to say, you
should be able to get an excellent data set with a narrow slit under excellent seeing conditions.

cl> display fits/IRCA00070352.fits 1 z1=0 z2=3000

14

Figure 4: An image of the standard star without any calibration (i.e., the raw data). The left andright panels, respectively, show the image taken at the slit positions A and B.
Check quality of the data taken at the position B using the above command. You will see animage shown in the right hand panel of Figure 4 where you should verify that the location of thesignal of the standard star has moved to the left.

5.2.3

Subtraction: A – B

The raw image, of course, contains not only the desired target signal but also the backgroundemission which is usually very strong in the L-band. To better visualize signals from the star, we
can simply subtract the background emission, namely, making an image through A-B. Issue thefollowing IRAF commands,

cl> imarith fits/IRCA00070351 - fits/IRCA00070352 ana/HR72_351_352
cl> display ana/HR72_351_352 1 z1=-1000 z2=1000

You will see an excellent spectroscopic (the resultant subtracted) image for the standard star as
shown in Figure 5 left. Here, ana/HR72 351 352 is the name of the output image which can bearbitrary. It would be wise of you to establish your own naming convention for ﬁle names. As seen
above, IRAF recognizes a ﬁle name without its extension (e.g., .FITS) as we have deﬁned in §5.2.1.
Also try to see how the images will be displayed with diﬀerent z1 and z2 values.
15

Figure 5: (Left) – Spectroscopic data of the standard star obtained from A-B. The mean value forthe source-emission-free region ideally should have zero count within the errors. (Right) – The sameas the left panel, but a sub-image containing only the source of interest; the sub-image has beenobtained by truncating the original.
In the following steps, we will work on the A-B image, but we will use solely the left hand portionof the image where the star of interest is visible. Therefore, it would be good idea to get rid of
the remaining two-third portion on the right hand side to save disk space. Another reason is thatthese region may contain ”jumped” pixel counts (i.e., showing an extraordinary high count) that
may not be handled correctly. For this purpose, let us truncate the left portion using the following
command:cl> imcopy ana/HR72_351_352[1:340,40:900] ana/HR72a

where the ﬁrst two values in the square brackets are the x-coordinates and the last two are they-coordinates for the image to be clipped. After clipping, make sure that you have selected thedesired portion of the image by,cl> display ana/HR72a 1 z1=-1000 z2=1000

We assume that you have an A-B image, like the right hand panel of Figure 5, that will be reducedin the following manner. In practice, almost all the observers repeat the observing sequence of A
16and B. Namely, we have more than one pair of A and B images that should be combined after making
A-B images from adjacent sequences. Since the standard star, HR 72, has been observed before andafter observing the target galaxy, we have 12 frames from IRCA00070351 through IRCA00070362(giving 6 sets of A-B) and 8 frames from IRCA00070413 through IRCA00070420 (giving 4 sets of
A-B). In this Cookbook, let us use the former set of data.
To make subtracted-images, execute the IRAF task imarith, as follows:cl> imarith fits/IRCA00070351 - fits/IRCA00070352 ana/HR72_351_352
cl> imarith fits/IRCA00070354 - fits/IRCA00070353 ana/HR72_354_353
cl> imarith fits/IRCA00070355 - fits/IRCA00070356 ana/HR72_355_356
cl> imarith fits/IRCA00070358 - fits/IRCA00070357 ana/HR72_358_357
cl> imarith fits/IRCA00070359 - fits/IRCA00070360 ana/HR72_359_360
cl> imarith fits/IRCA00070362 - fits/IRCA00070361 ana/HR72_362_361

Consequently all the six images should be combined into an image in order to improve the S/N.
For this purpose, we decided to sum up the images straightforwardly. To combine the images,
you can also take another approach that calculates the mean. Notice that both methods should
essentially work the same from our point of view; the resultant images created from summation andaveraging only have a diﬀerence in their scale factors.
Open a blank ASCII ﬁle typing ”emacs HR72 list &” on the Unix/Linux command line. The
ﬁle should contain the following lines:

ana/HR72_351_352.fits
ana/HR72_354_353.fits
ana/HR72_355_356.fits
ana/HR72_358_357.fits
ana/HR72_359_360.fits
ana/HR72_362_361.fits

Note that you may omit the ”.fits” extension. Use the ASCII text ﬁle as an input ﬁle to sum
up the 6 images,
17

cl> imsum @HR72_list ana/HR72_351_362

Did you get an improvement of the S/N ratio for the standard star image? If you are happy
with the results, it would be prudent to throw away the unnecessary region by issuing the following
command:cl> imcopy ana/HR72_351_362[1:340,40:900] ana/HR72_1

Just for your knowledge, we do not need to subtract dark-current images, as we are going to
reduce the A-B images (in fact, the dark current and bias level stay at a negligible level for the
majority of the L-band observations).
If you are familiar with IRAF scripts, you may want to skip some steps that are often repeatedsuch as ”subsum fits.cl”, which is found in ”ircs sp/tools”. To use the IRAF script, follow
below:

cl> task subsum_fits = tools/subsum_fits.clcl> subsum_fits IRCA0007 351 362 HR72
cl> imcopy ana/HR72_351_362[1:340,40:900] ana/HR72_1

Here, you may see a write permission warning for the output ﬁles. In this case, delete the ﬁlesbefore using the script:

cl> imdel ana/HR72_1

18

The script will help you to increase your work eﬃciently. We hope that you will be able to write
such scripts after familiarizing yourself with IRAF. However, if you consider yourself as being anovice, we suggest repeating each command and inspecting all the results by eye.
Subsequently, let us apply the same procedure4 on the target ultra-luminous infrared galaxy
IRAS 00188−0856. Although we have a full set of data set (ﬁle number = IRCA00070368 through
IRCA00070407), we are going to use the subsets of the ﬁles here i.e., between IRCA00070368
and IRCA00070387 for illustrating how to make a subtracted (= A-B) ﬁle that will be named
IR00188 368 387.fits5 After that, we will clip only the area in the same fashion as applied to
the HR 72 image (see Figure 6).

Figure 6: An example of spectroscopic raw data generated by A-B for the target source. The image
shown here is the clipped sub-image. Compared with Figure 5, the image toward the infrared galaxy
image is noisier than that of the standard star.

5.2.4

Dividing with Spectroscopic Flat

Ideally, both the standard star and the target should be observed at the identical slit positions.
However, in practice, it is diﬃcult to match exactly the position of the standard star and the4

The task subsum fits requires a rather large amount of memory, making the system unusually slow. When you
attempt to combine large numbers of images using a computer with a limited amount of memory and encounter sucha problem, try to reduce the number of images being combined.
5
Since the observations were carried out by repeating the sequence of ABBA, the FITS ﬁle numbers of IRCA00070368
and IRCA00070369 correspond to the slit positions A and B, respectively. The sub-data sets ends with IRCA00070387
which is taken at A.

19

source(s). Moreover, we cannot exclude the possibility that slit width varies with slit position
and/or that the transmission is not uniform over the entire IRCS optics system along the slit. In
order to eliminate these artifacts, we correct for the transmission along the slit using an image
obtained by observing spatially uniform incident light. We have to make a diﬀerential image bysubtracting the reference lamp ON (i.e., the light from the lamp plus thermal emission of the system)from OFF (the system thermal emission only). Hence, the subtraction should give an image that
contains the spatially uniform incident light emanated from the lamp.
Since we now have the ﬁve sets of lamp ON (IRCA00070227 – IRCA00070231) and OFF(IRCA00070247 – IRCA00070251) images, we should make a diﬀerential image for each pair.
Here, you may see an error when running imcombine command in IRAF. This error is due toimproper writing of the FITS header by IRCS (unfortunately this has not been ﬁxed yet). Pleasecheck if the header parameter CTYPE2 is LINEAR. If this is the case, an extra step to modify theheader information is needed.
Creat ﬁles named ”flaton” and ”flatoff”, using an editor. The ﬁle of flaton contains:fits/IRCA00070227.fits
fits/IRCA00070228.fits
fits/IRCA00070229.fits
fits/IRCA00070230.fits
fits/IRCA00070231.fits

The ﬁle of flatoff contains

fits/IRCA00070247.fits
fits/IRCA00070248.fits
fits/IRCA00070249.fits
fits/IRCA00070250.fits
fits/IRCA00070251.fits

Let’s see and modify the input parameters for the hedit task:
20

cl> epar hedit
imagesfields
value

=
=
=

@flatonCTYPE2
DEC---TAN

images to be edited
fields to be edited
value expression@flatoff
CTYPE2
DEC---TAN

images to be edited
fields to be edited
value expressioncl> heditcl> epar hedit
imagesfields
value

=
=
=

cl> hedit

The hedit task asks you many questions at each step. If this is annoying, you can select thesilent mode by changing verify to no. Parameter editing can be done with epar hedit.
Subsequently, we will subtract the OFF ﬁles from each corresponding ON ﬁle, then extract a
region that contains the target object (i.e., left portion of the detector). For combining the data,
our experience suggests adopting the so-called ”median” ﬁltering rather than averaging one (theoption ”average”). This is because, cosmic rays may hit during an exposure, causing pixels whosecounts are enormously high. The presence of such pixels would skew the resultant average value
improperly, while the ”median” ﬁltering avoids such cases.

cl> imcombine @flaton ana/flaton_227_231 comb=mediancl> imcombine @flatoff ana/flatoff_247_251 comb=mediancl> imarith ana/flaton_227_231 -ana/flatoff_247_251 ana/flat1

cl> imcopy ana/flat1[1:340,40:900] ana/flat2

21

After getting spectroscopic ﬂat data, it would be convenient to normalize its intensity scale,although this is not a strict requirement as the result of not normalizing is just an image that isscaled by a constant value.

cl> epar imstat
Image Reduction and Analysis Facility
PACKAGE = imutil
TASK = imstatistics
images =ana/flat2) List of input images(fields = image,npix,mean,midpt,stddev,min,max) Fields to be printed(lower =INDEF) Lower limit for pixel values(upper =INDEF) Upper limit for pixel values(nclip =0) Number of clipping iterations(lsigma =3.) Lower side clipping factor in sigma(usigma =3.) Upper side clipping factor in sigma(binwidt=
1.0000000000000E-5) Bin width of histogram in sigma(format =yes) Format output and print column labels ?
(cache =no) Cache image in memory ?
(mode
=
ql)

Add midpt to the fields parameter, and edit image,npix,mean,midpt,stddev,min,max asshown above.

cl> imstat ana/flat2
IMAGE
ana/flat2

NPIX
292740

MEAN
1440.

MIDPT
1520.

STDDEV
365.3

MIN
-58.

MAX
2459

Editing the parameters as shown above, you now are ready to normalize the frame with its
median value (= MIDPT).
22

cl> imarith ana/flat2 / 1520.0 ana/flat

Although you should learn by trial and error, we suggest that a small value of binwidt (i.e., thewidth for each bin of the histogram) will give a reasonable median value.
The next step is to divide the object and the standard star frames by the normalized spectroscopic ﬂat.

cl> imarith ana/HR72_1 / ana/flat ana/dHR72_1
cl> imarith ana/IR00188_1 / ana/flat ana/dIR00188_1

Here, the IR00188 1.fits should be created from the data set described in p.19 with the same
as for the standard star. In practice, ﬁnal IRCS spectra that have been divided by spectroscopic ﬂat
do not diﬀer signiﬁcantly from ﬁnal IRCS spectra that have not been divided by spectroscopic ﬂat.
We strongly encourage readers to verify the results of with and without spectroscopic ﬂat-ﬁelding.

5.2.5

Eliminating Bad Pixels

It is known that some of the pixels within the 1024×1024 IRCS detector arrays do not behave as
expected. Some of them permanently do not work. They always generate enormously high pixel
counts (i.e., hot pixel) or no-signal (i.e., dead pixel). Contrary to these ”bad” pixels, the majorityof the working pixels — ”good” pixels — generate output signals whose levels are proportional to
the numbers of the incident photons. As we described in the previous section, even ”good” pixels
sometimes may return extremely high counts after being hit by cosmic rays. Clearly, data fromthese pixels should not be used in any data reduction process. In the following, we describe how to
identify and eliminate such bad pixels. Let’s go back to IRCA00070351.fits (see Figure 4 left) forthe purpose of inspecting pixel behavior.

cl> display fits/IRCA00070351 1 z1=0 z2=30000

23

You may recognize that there are a few compact white points around the image center. Moreover,
one can identify a distinct white point (i.e., showing a very high count value compared to its
surrounding pixels) both on the left and right portions of the image. Accordingly we should correctfor (i.e., replace) these bad pixels with the appropriate ones, by extrapolating from the surrounding
pixel values. For this purpose, a general approach is to replace the bad pixel values with either themean or the median of the adjacent 8 pixels. We believe that this approach works for many cases.
But, if you have a preference, try an alternative approach.
In principle, you can identify bad pixels by eye by typing the following command:cl> display ana/dHR72_1 1 z1=-1000 z2=1000

Ideally, bad pixel correction should be done on each individual frame. However, it will be tediousto inspect very large numbers of the frames because, in general, frame rates of infrared observations
are high. Here, we will apply bad pixel corrections onto the data set combined from several frames.
We will use the IRAF task, cosmicrays on dHR72 1.fits to identify and correct for the bad
pixels that have elevated values compared to their surrounding pixels. The task is implemented in
the crutil package under noao/imred:

cl> noaono> imred
im> crutil
cr> epar cosmicrays(interac=no) Examine parameters interactively?

This allows you to use the task cosmicrays6 .
6For your knowledge, once you have executed cosmicrays at cr>, you can issue the command from any prompt.
Namely, not only at the cr> prompt, but also at the e.g., cl> and no> prompts.

24cr> cosmicrays ana/dHR72_1 abc

Notice that cosmicrays command is designated solely to ﬁnd bad pixels having positive values.
You may want to ﬁnd the negative bad pixels as well, e.g., in the case of the A-B image. For thispurpose, just multiply −1 to the image before running cosmicrays, as shown below.

cr> imarith abc * -1 abc1
cr> cosmicrays abc1 abc2
cr> imarith abc2 * -1 ana/cHR72_1

Comparing the original dHR72 1.fits with the corrected cHR72 1.fits, you will see that thelatter image has smaller numbers of bad pixels than the former. The task cosmicrays has succeededin eliminating almost all the bad pixels. If you are not satisﬁed with the results, you can perform
further bad pixel corrections by hand using the task fixpix. The following shows how to accomplish
this:cr> imcopy ana/cHR72_1 ana/fHR72_1
cr> fixpix ana/fHR72_1 badHR72_1

The ﬁrst line just makes a copy of the original data as fixpix overwrites the original image.
The second operation corrects the count values of the pixels whose coordinates are given in the
badHR72 1 ﬁle, which should contain pixel coordinates:

18 18 97 97
23 24 67 69

25

This tells fixpix that count values for the pixel located at (x, y) = (18, 97) and (23 to 24, 67
to 69) should be replaced with those extrapolated from their surroundings.
If you want to give a speciﬁc value to replace the bad pixel, instead of using extrapolation, thetask epix allows you to do this:cl> epix ana/fHR72_1 18 9717
18
19
96 -39.0465 -7.45007 44.8567
97 2.20168 -51.8033 -30.5948
98 -17.9674 -55.1241 -52.4665
median -30.59483, mean -23.04381, sigma 32.69519, sample 9 pixels
new value for pixel (-51.8033): 5.0

The above example replaces the count at (X, Y) = (18, 97) with 5.0.
Last, don’t forget to apply the bad pixel corrections to the target data (i.e., dIR00188 1.fits
to make cIR00188 1.fits and ana/fIR00188 1.fits.

5.2.6

Extracting Signals and Combining Data taken at Diﬀerent Slit Positions

We assume that you have successfully completed all the calibration steps so far. Now we are readyto extract spectra of the standard star and the target galaxy. The extraction can be done with thetask apall in IRAF.
If you want to learn the detail on how apall works on single-order spectroscopic data, we suggest
reading ”A User’s Guide to Reducing Slit Spectra with IRAF” by Phil Massey, Frank Valdes, &Jeanette Barnes (1992 April; http://iraf.noao.edu/docs/spectra.html).
cl> noaono> twodspec
tw> apextract26

ap> apextract.dispaxis=1 (This is just to make sure that the vertical axis
represents the wavelength dispersion.)ap> epar apall
I R A F
Image Reduction and Analysis Facility
PACKAGE = apextract
TASK = apall
input=(output =(apertur=(format =(referen=(profile=ana/fHR72_1 List of input images) List of output spectra) Apertures
multispec) Extracted spectra format) List of aperture reference images) List of aperture profile images(interac=(find=(recente=(resize =(edit=(trace =(fittrac=(extract=(extras =(review =yes) Run task interactively?
yes) Find apertures?
yes) Recenter apertures?
yes) Resize apertures?
yes) Edit apertures?
yes) Trace apertures?
yes) Fit the traced points interactively?
yes) Extract spectra?
no) Extract sky, sigma, etc.?
yes) Review extractions?

(line
(nsum

400) Dispersion line20) Number of dispersion lines to sum or median==# DEFAULT APERTURE PARAMETERS(lower(upper==-5.) Lower aperture limit relative to center5.) Upper aperture limit relative to center
27

(apidtab=) Aperture ID table (optional)# DEFAULT BACKGROUND PARAMETERS(b_funct=(b_order=(b_sampl=(b_naver=(b_niter=(b_low_r=(b_high_=(b_grow =chebyshev) Background function
1) Background function order-40:-30,30:40) Background sample regions-3) Background average or median0) Background rejection iterations
3.) Background lower rejection sigma
3.) Background upper rejection sigma0.) Background rejection growing radius# APERTURE CENTERING PARAMETERS(width =(radius =(thresho=5.) Profile centering width
10.) Profile centering radius
0.) Detection threshold for profile centering# AUTOMATIC FINDING AND ORDERING PARAMETERS

nfind=(minsep =(maxsep =(order =1 Number of apertures to be found automatically5.) Minimum separation between spectra
1000.) Maximum separation between spectra
increasing) Order of apertures# RECENTERING PARAMETERS(aprecen=(npeaks =(shift =) Apertures for recentering calculation
INDEF) Select brightest peaksyes) Use average shift instead of recentering?
# RESIZING PARAMETERS(llimit =INDEF) Lower aperture limit relative to center
28

(ulimit =(ylevel =
(peak
=
(bkg
=
(r_grow =(avglimi=

INDEF) Upper aperture limit relative to center
0.1) Fraction of peak or intensity for automatic widtyes) Is ylevel a fraction of the peak?
yes) Subtract background in automatic width?
0.) Grow limits by this factorno) Average limits over all apertures?
# TRACING PARAMETERS(t_nsum =(t_step =(t_nlost=(t_funct=(t_order=(t_sampl=(t_naver=(t_niter=(t_low_r=(t_high_=(t_grow =4) Number of dispersion lines to sum4) Tracing step3) Number of consecutive times profile is lost befo
legendre) Trace fitting function
3) Trace fitting function order*) Trace sample regions
1) Trace average or median0) Trace rejection iterations
3.) Trace lower rejection sigma
3.) Trace upper rejection sigma0.) Trace rejection growing radius# EXTRACTION PARAMETERS(backgro=(skybox =(weights=(pfit=(clean =(saturat=(readnoi=(gain=(lsigma =(usigma =(nsubaps=(mode=fit) Background to subtract
1) Box car smoothing length for skynone) Extraction weights (none|variance)
fit1d) Profile fitting type (fit1d|fit2d)no) Detect and replace bad pixels?
INDEF) Saturation level0.) Read out noise sigma (photons)
1.) Photon gain (photons/data number)
4.) Lower rejection threshold4.) Upper rejection threshold1) Number of subapertures per apertureql)

29

apall has many parameters, but most of them are not critical. Bring your cursor to the ﬁeld you
want to change and type ”?”; you will get options that can be selected. Set the following ”critical”
parameters, explained below, as follows: line = 400, nsum = 20, t nsum = 4, t step = 4,t order = 3, b sample = -40:-30,30:40, backgro = fit, weights = none, clean = no, andextras = no. Here, extras = no lets the program extract only signals in the speciﬁc case of thesample data.
Assuming that the spatial proﬁle of the signal does not change with wavelength, we will checkit at a certain wavelength (i.e., at a y position where y is the vertical coordinate). line = 400 andnsum = 20 speciﬁes the location of the spatial proﬁle from the pixel value y = 400±10. t nsum =4 sums 4 pixels along the wavelength direction. t step = 4 sets the search for the peak position of
the signal in the x coordinate, by shifting 4 pixels along the wavelength axis (= the y coordinate).
t order is the order of the function that will be used to ﬁt the peak position of the signal (atthe x coordinate). In general, ﬁtting will be successful with t order = 3. Notice that t order =3 speciﬁes a quadratic equation, i.e., y = Ax2 + Bx + C. We selected the order set to 3 for the
standard star because the S/N of the sample data is very high, allowing us to trace the peak x value
with high accuracy. On the other hand, ﬁtting with a linear equation (t order = 2; y = Ax + B)
usually works well for scientiﬁc target as they are faint.
The counts of the A-B image, except for the signal counts, should ideally be scattered around thezero level within the noise. However, there often exist pixels whose values are non-zero (taking somepositive or negative values), and whose absolute values vary with wavelength. These non-zero values
indicate that the sky background level has changed during the exposure at the A and B positions
7. If you specify backgro = fit, you can further subtract this residual pattern. b sample is the
distance in the x-coordinate from the peak position that deﬁnes a source-emission-free area where
the sky-background level can be estimated. If you deﬁne the sky area on both sides of the peak
position in the x-coordinate, then the linear sky gradient in the x-direction can be corrected for.
weights = none and clean = no extract a spectrum by summing the signal using the moststraightforward method instead of using the optimal extraction. Here, optimal extraction is amethod to add the signals at each pixel by giving a weight based on the noise at each pixel. If
optimal extraction works well, you can expect to obtain an improved S/N ratio, compared to simple
summation (see http://archive.stsci.edu/imaps/expastro/node29.html). Optimal extraction
is known to work well on optical spectroscopic data. This is because the Earth’s atmospheric
7

At the wavelength where the Earth’s atmospheric transmission is intrinsically small, the background emission ishigh (Recall 1 - transmission = absorption = emission), where a small diﬀerence of the relative sky level can resultin high absolute value.

30

transmission is almost 100% over the band. Another factor, thanks to the above reason, is thatmost of the optical signals tend to show smooth proﬁles along the wavelength axis. On the contrary,
the Earth’s atmospheric transmission in the infrared strongly depends on wavelength. Hence theintensity of the signal varies considerably with wavelength. We have witnessed a couple of obvious
examples where optimal extraction produces ”strange” spectra accompanied by unusually large
scatter of the data points especially at the wavelengths where the Earth’s atmospheric transmission is
intrinsically low. Therefore, we suggest simply summing signal over the aperture without weighting.
If you plan to use optimal extraction, you should carefully read the original documents. Moreover,test with various combinations of the input parameters. Don’t forget to compare the results using
optimal extraction with those using the standard extraction (i.e., without weighting). Last, please
note that optimal extraction executs automatically with clean = yes even if you specify weight= none.
It is highly likely that the results will be essentially the same even if you try to change the otherparameters in apall that are mentioned above. Of course, there may be exceptions, e.g., in thecase of a very faint source where apall may fail to trace a spectrum unless you have optimized the”minor” parameters. If this is the case, try to ﬁnd the best parameters for your source by checking
the help page (type help apall to get the help page). Shown below is apall usage:ap> apall ana/fHR72_1 out=ana/fHR72_1posFind apertures for ana/fHR72_1? (yes):
Number of apertures to be found automatically (1):Resize apertures for ana/fHR72_1? (yes):
Edit apertures for ana/fHR72_1? (yes):

You will see that a Tektro window shown in Figure 7 pops up.
On the Tektro window, you can optimize graphically the range used for the spectral extraction.
Move your cursor onto the window, and type e.g., ”:low -8”, ”:upper 8”, and ”:center 215”.
Once you have provided a range to the program, type ”t” to execute ﬁtting. In Figure 7, the twohorizontal (short) lines (around x = 180 and 250 in this example) with very short vertical bars are
the ranges over which the sky-level will be calculated for its subtraction. Make sure the ranges are
OK, then, press ”Yes”. The left hand panel of Figure 8 shows how the ﬁtting works, and the righthand shows the extracted (i.e., one-dimensionalized) spectrum. If you are happy with the spectrum,31

Figure 7: The output window for apall. This window is referred to as the Tektro window.
return to IRAF command line (cl>) by typing ”q”. Use caution when you repeating apall on thesame input ﬁle, as apall overwrites the results to the same output ﬁle. You may have seen a queryfrom Clobber *** at a later stage of apall execution where the default setting of ”No” prohibits
overwriting of the results. If you want to save the second (or reﬁned) ﬁtting results after changing
the parameters, don’t forget to set the parameter to ”Yes”.
Recall that we have to extract the signal observed at slit position B in order to combine themto improve the S/N. Repeat the above procedures for the slit position B as follows:ap> imarith ana/fHR72_1 * -1 ana/fHR72_1n
ap> apall ana/fHR72_1n out=ana/fHR72_1neg
ap> imarith ana/fHR72_1pos + ana/fHR72_1neg ana/fHR72_1ana

After combining both the A and B spectra of the standard star, the next step is reducing thetarget galaxy. It is certainly prudent to set line=550 and nsum=40,t order=2 because the source
is faint. Trace the signal as shown in Figure 9, and make the resultant ﬁle fIR00188 1ana.fits.
32

Figure 8: An example of the working windows that will pop up when running apall. The left panel
shows how apall traces a spectrum. In this example, notice that the ﬁtting curve (the dashed line)
has obviously failed to trace the data at ”Line” (the x coordinate) & 600 (the ﬁtting curve is shiftedby up to ∼0.25 pixel). You may have a better result by increasing the order of the ﬁtting function.
The right panel indicate the one-dimensionalized, i.e., the extracted spectrum.

Figure 9: Results from apall of the target galaxy IRAS 00188. The left- and right hand panels,
respectively, show the spectral ﬁtting and the spectrum obtained from the target galaxy.
Now, we have obtained a one-dimensional spectrum of the target galaxy. You may elect tocontinue to work with IRAF for further analysis such as binning and calculating errors. However,
we must admit that IRAF is not necessarily the best tool for such detailed analysis. Because ofthis, many users seem to prefer transferring the extracted one-dimensional spectra into their own
programs or/and other software packages.
To write the extracted spectrum into an ASCII text ﬁle, use the IRAF task wspectext as follows:ap> noao

33

no> onedspecon> epar wspectext
I R A F
Image Reduction and Analysis Facility
PACKAGE = onedspec
TASK = wspectext
input=
output =(header =(wformat=(mode=

fHR72_1ana

Input list of image spectra
Output list of text spectrano) Include header?
) Wavelength format
ql)

Don’t forget to select ”no” for the header ﬁeld, otherwise you may get long, unnecessary, header
information.

on> wspectext ana/fHR72_1ana spec/fHR72_1.qdp

Here the extracted ﬁles are stored under ”spec” directory. The same procedure should be
repeated for the target galaxy as well. Here, the extension of the ﬁle, .qdp, is arbitrary8 . We prefer
to use it as it is convenient for the later steps.
If your source is too faint to be traced by apall, you can use data of a bright source (i.e.,
standard stars) as a reference to tell apall how to trace the data. Such a reference image can besupplied to apall with the ”referen” parameter. In the above example, you may use it as, e.g.,
apall ana/fIR00188 1 out=ana/fIR00188 1pos refe=ana/fHR72 1. Keep in mind that, in order
to apply the results from the bright source (standard star), both the bright and the dim (target)
sources must have been observed at the same slit position. If not, you have to shift either of theimages slightly along the x-direction, using the IRAF command ”imshift”.
8

To use qdp, add the following path in your login shell. If the case of bash, add the following to .bashrc:
PATH=$PATH:/usr/local/headas-6.5/x86 64-unknown-linux-gnu-libc2.5/bin,and execute:source/usr/local/headas-6.5/heagen/x86 64-unknown-linux-gnu-libc2.5/headas-init.sh.

34

5.2.7

Correction for Non-linearity

The ideal detector should generate a signal whose intensity is perfectly proportional to the numbersof the incident photons. However, infrared detectors (CCDs at the optical regime as well) in practice
do not show such a behavior. Below, we attempt to explain how the infrared detectors work by
returning to the principals, with the understanding that its mechanism cannot be so idealized. If
you are interested in the details, please refer to ”Electronic Imaging in Astronomy Detectors andInstrumentation” by Ian S. McLean in the WILEY series.
Detectors produce a depletion layer by setting an inverse-bias voltage on the pn-junction which
consists of p- and n-type semiconductors. A photon incident to the depletion layer will generate anelectron9 that causes an electron and a hole to apparently move in opposite directions from eachother. This reduces the inverse-bias voltage. If there are a large number of the incident photons andif the inverse voltage drops to near zero, the depletion layer disappears, causing a loss of sensitivity
to incident photons. Such a saturated situation corresponds to the so-called ”full-well” level of thedetector.
In other words, the above mechanism may be compared to how a capacitor with a capacity, C,
works. If charge Q is stored at both the sides of the polar plates, a voltage of V = Q/C will be
generated. The Q value will vary with the amount of the incident photons, i.e., the number of theelectrons generated by the photo-electric eﬀect. If the capacity, C, remains constant, the variationof the voltage should be proportional to the number of the generated electrons, unless the detector
is saturated. In summary, the detector counts the number of the incident photons by counting
those generated electrons, namely, the variation of the voltage. Variation of the voltage is an analog
quantity which is digitized by the AD converter. In other words, we are dealing with digitized countnumbers in ADU, as seen in the image viewer.
Figure 10 shows that the generated (output) voltage tends to decrease once the detector has
stored a large number of electrons. In other words, the output voltage does not increase proportionally to the numbers of the incident photons. In the above explanation, this would correspond to asituation where the capacity of the capacitor has increased as the number of the stored electrons has
increased. The diﬀerence between the actual and the expected output voltage increases abruptly
once the incident photon numbers exceeds a certain threshold (see Figure 10). Below this threshold,the diﬀerence stays at a negligible level. For this reason, one must select a proper exposure time
for the detector used in observations so as not reach the threshold value. In the speciﬁc case of9

If the detector outputs an electron for an input photon, its quantum eﬃciency (QE) is referred to as 100%. The
infrared detectors currently being used are considered to have an eﬃciency of 70 – 90%, which is high enough to becomparable to the CCDs used in the optical regime.

35

IRCS, the oﬀset values with respect to the ideal linear-response is known to be at most 1–2 %, if
the ADU count is conﬁgured to be less than 4000. If your data are low-resolution spectroscopy datawith a modest S/N (e.g., ∼ 10 − 30), the linear deviation is certainly smaller than the noise level.
Therefore, correction for the non-linearity should not be so critical.

Figure 10: A plot of the ideal (the dashed line) and the actual (the thick curve) detectors’ responses
as a function of the incident photon numbers.
With the following explanation, it is possible to understand the mechanism of why the capacity,
C, increases with the numbers of the stored charges. It is known that the capacity, C, is representedby ∝ S/d where S and d are the surface area and separation of the two polar plates, respectively.
The decrease of the depletion layer by the incident photons can be thought of as a decrease in theseparation of the dual polar plates, i.e., the decrease of d causes an increase of C, explaining theactual behavior of detectors.
Although you must verify the signal level before starting the data reduction, let us checkthe signal levels of the IRCS data. Notice that the count value of the IRCS data have been multiplied by the product of the NDR and COADD values, so we have to divide the data by this productbefore checking the count values (in ADU unit). Let’s take IRCA00070351.fits as an example.
The ﬁts header tells us that the data were taken with NDR = 2 and COADD = 10; therefore, weshould divide the image by a factor of 20 (= 2 × 10), as follows:cl> imarith fits/IRCA00070351 / 20.0 fits/IRCA00070351d
cl> display fits/IRCS00070351d 1 z1=0 z2=4000
36

If you move your cursor onto the ds9 window, you get the count value in ADU unit for eachpixel. If you want to have more detailed information, try this:cl> imexamine
cl> z
Have you obtained the count values for the region of interests? Type ”q” to go back to the cl>prompt. Carefully check whether or not all the pixels have counts below 4,000, except for those in
the uppermost (long wavelength) region of the image and some hot pixels10 . We also need to checkthis for IRAS 00188. Keep in mind that the object frames have NDR = 6 and COADD = 15.
cl> imarith fits/IRCA00070368 / 90.0 fits/IRCA00070368d
cl> display fits/IRCS00070368d 1 z1=0 z2=4000
Last but not least, we strongly recommend checking the count values for all the ﬁles.

5.2.8

Wavelength Calibration

In order to precisely associate the pixel coordinate with the absolute wavelength, it is ideal toobserve the light from a lamp that emits many lines whose wavelengths are well-known and line
widths are small. Unfortunately, there is no such an ideal lamp for wavelength calibrations in theinfrared L band (λ = 2.8 − 4.2 µm) spectroscopy with IRCS, unlike the optical regime and shortwavelength (λ . 2.5 µm) infrared. We therefore utilize the transmission spectrum of the Earth’s
atmosphere for wavelength calibration. This situation is similar at the longer wavelength M band
(λ = 4.5–5.0 µm) as well.
Figure 11 enlarges the wavelength dependency of the Earth’s atmospheric transmission (Figure2) in the L band (λ = 2.8 − 4.2 µm); one can utilize the strong (intrinsic) wavelength dependenttransmission for this purpose. Of course, this concept can be applied to the data taken at the shortwavelength, i.e., some observers use the very peaky OH airglow (vibration-rotation lines) emanatedfrom Earth’s atmosphere.
The pixel values showing counts higher than 4000 ADU are the data for λ ∼ 4.2 µm in the case of the sample
data. The data around this wavelength range are not our interest. We (M. Imanishi), therefore, have intentionally
set the signals at λ ∼ 4.2 µm above the 4000 ADU level (i.e., long exposure), to maximize the observing eﬃciency.
1037

Figure 11:
The Earth’s atmospheric transmission curve for the L band (λ = 2.8–
4.2 µm).
This plot is a enlarged version of Figure 2.
Image courtesy of the UKIRT(http://www.jach.hawaii.edu/UKIRT/astronomy/utils/atmos-index.html).
If you want to know the exact wavelengths of the peaks in Figure 11, you can get ASCII formateddata from http://www.jach.hawaii.edu/UKIRT/astronomy/utils/atmos-index.html and clickon ”Text version”.
For a simplicity, we are going to make a linear ﬁt to the data using the methane absorption line
at 3.315 µm and the emission at 3.9005 µm for many cases. You can try to ﬁt the data with a highorder functions using least-square ﬁtting. However, the resultant accuracy would not improve somuch.
Subsequently, we will compare the spectrum of the bright standard star with the intrinsic transmission spectrum in the ASCII text format. In the following steps, we will abandon IRAF, anduse software written by one of the authors (Imanishi) in gawk. Of course, you can use any existing
software you like, but we recommend to writing your own for a better understanding.
Since the ﬁrst column of the extracted ASCII data (fHR72 1.qdp) does not have line numbers
staring from 1, we have to correct it by:
gawk -f tools/pix fHR72_1.qdp > fHR72_1c.qdp
Here pix contains a simple sentence of{printf("%5.6f %5.6f\n",NR,$2)}38

where NR is the system variable implemented in gawk.
The next step is to verify the spectra by displaying with qdp (or your favorite, e.g., SuperMongo,
gnuplot, IDL, etc). Here, you can obtain qdp fromhttp://wwwastro.msfc.nasa.gov/qdp/. In this Cookbook, we present basic concepts of further
analysis on the extracted spectrum using qdp, as we are familiar with qdp.
qdp fHR72_1c.qdp/xs (or /xw)
which provides the plot shown in the left hand panel of Figure 12. Comparing the two panels
in Figure 12, the very narrow dip at x = 310 corresponds to 3.315 µm, and the peak aroundx = 680 with 3.9005 µm. Such a correspondence can be obtained by e.g., checking the ASCII
text ﬁle by opening it with an editor, or by magnifying the data around x = 300 to 400 with qdp
command ”r x 300 400”. If we ﬁt the data with a linear equation, the best-ﬁt is represented by λ= 1.5975×10−3 x + 2.8174. Using the results, let us associate the ﬁrst column from the row number
with the absolute wavelength. Here we will use a program written by us,

gawk -f tools/wavecal spec/fHR72_1c.qdp > spec/fHR72_1cw.qdp

where wavecal contains{printf("%5.6f %5.6f\n",$1*0.0015975+2.8174,$2)}

It is known that the Earth’s atmospheric transmission at wavelengths shorter than 2.8 µm
degrades signiﬁcantly. Moreover, the background emission level around λ ∼ 4.2 µm is too high forthe linear response range (Recall §5.2.7) 11 . In fact, there are some data points exceeding 4000 ADU
at λ = 4.10 − 4.15 µm in the object frame. We therefore delete the wavelength ranges where theEarth’s atmospheric transmission is intrinsically small, and the counts exceeds 4000 ADU (λ ∼ 4.2µm), as follows.
The lines of our interests exist in λ . 4.1 µm, hence the excess over the linearity range does not harm ourscientiﬁc discussion. If we want to obtain the galaxy’s spectrum up to ∼ 4.2 µm, we need to shorten the exposure
time. However, this makes the observing eﬃciency low.
1139

cp spec/fHR72_1cw.qdp fHR72_1ccw.qdp (Copying the file, just for the sake of saving)
emacs fHR72_1ccw.qdp & (delete unnecessary rows)

Later on it will be convenient, when spectral channel binning (§5.2.10), to have a pixel numbers
be multiples of 8.

qdp fHR72_1ccw.qdp/xs

We suppose that you obtained a spectrum whose x-axis has been converted into wavelength,
similar to the right hand panel of Figure 12.

Figure 12: Spectra of the standard star HR 72 before (left) and after (right) the wavelength calibration.
If you are satisﬁed with the standard stars spectrum, apply the same procedure to the target
galaxy to obtain the wavelength-calibrated spectrum shown in Figure 13.
For cases where you have performed spectroscopic observations at wavelengths shorter than 2.5µm using a calibration lamp, you may use the IRAF tasks, e.g., identify, refs and/or dispcor,
which are optimized for wavelength calibration. Please see the IRAF manual for more details.
40

Figure 13: The wavelength calibrated spectrum of the object galaxy, IRAS 00188 (z = 0.128). The
peak emission seen at 3.7 µm is considered to be the redshifted PAH emission whose wavelength in
the rest-frame is 3.3 µm.
5.2.9

Airmass Correction

The path-length of where the astronomical signal has traveled through the Earth’s atmospherevaries with airmass. However, as long as you have observed both the target(s) and standard stars
within an airmass diﬀerence of less than ≈ 0.1, you may not need to apply such corrections for
low-dispersion grism spectroscopic data taken with IRCS with an S/N . 10 − 30.
However, if your data have substantially good S/N, e.g., exceeding ∼ 30, you probably realize
a need for applying airmass corrections. Please see IRAF documents and the other resources fordetails.

5.2.10Binning Spectrum Channels

If you have observed with a slit width of 0.′′ 9, this corresponds to 16 pixels on the IRCS detector.
This deﬁnes an eﬀective wavelength resolution, R ∼150. On the other hand, you will realize, e.g.,
by opening the extracted ASCII-text ﬁle, that the extracted spectrum consists of a numerical pair of
wavelength and the signal intensity at EACH pixel of the detector. This spectral pixel (i.e., channel)
increment is 16 times larger than the eﬀective wavelength resolution, namely, oversampled. We,
therefore, perform binning of the data along the wavelength axis at 16-pixel intervals. According to
the Nyquist’s sampling theorem, information that is two times R is sensible (i.e., having independent
information). Let us apply spectral channel binning of each 8-pixel group. Although you can do41this with your favorite software, we show an example using the software that is written by us (M.
Imanishi).
cd spec (going to the spec directory)
tools/bin2_2col fHR72_1ccw.qdp > fHR72_1ccwb2.qdp
tools/bin2_2col fHR72_1ccwb2.qdp > fHR72_1ccwb4.qdp
tools/bin2_2col fHR72_1ccwb4.qdp > fHR72_1ccwb8.qdp
Since the software, bin2 2col, performs binning at every two pixels (i.e., channels), we haveto run it three times12 . In addition, the program assumes that the user will supply data having
an even number of pixels. So it would be the best to select your data with a multiple of 8 before
running the program.
If your source is faint, you may consider applying further smoothing. Generally speaking, many
L band spectra do not have good S/N ratio. In addition, they show a bumpy shape because theEarth’s atmospheric transmission is intrinsically low at that wavelength, especially at λ . 3.4 µm(recall Figure 11). Therefore, many observers prefer to perform channel binning particularly for the
short-wavelength regime.
For the sources whose accumulated signal photons are deﬁcient, you may need smoothing thespectrum as many times as possible to improve the S/N and/or reduce the bumpy features seen in
the proﬁle. However, if you smooth spectra too much, you lose the needed wavelength resolution
required for your scientiﬁc goal. As we have seen, there is no way to improve such a situation (unless
you integrate further using more telescope time) for the data whose amount of the signal photons
is intrinsically small.
If you continue to work with IRAF, try the task blkavg to bin a spectrum extracted by apall,see the IRAF help page for the detail.

5.2.11Dividing the Object Spectrum by Standard Star Spectrum

Now, we have wavelength-calibrated spectrum whose intensity represents (i.e., is convolved with) a
product of not only the intrinsic spectrum of the object but also the other factors. These factors
12

It is prudent to check line numbers of the output ﬁles using e.g., Unix/Linux command wc.

42

are the Earth’s atmospheric transmission curve, and the various wavelength-dependent response
functions of the instrument itself. Clearly, we have to eliminate contributions from these factors to
extract the intrinsic spectrum of the source. This correction can be done as follows,

paste fIR00188_1ccwb8.qdp fHR72_1ccwb8.qdp > tmp1
gawk -f tools/div2 tmp1 > fIR00188_fHR72_1b8ccw.qdp

The above operation by the gawk commands (written in div2 ﬁle) divides the target source databy that of the standard source. Here, paste is a Unix/Linux command (see its man page for usage).
For your knowledge, the div2 ﬁle consists of a very straightforward command of{printf("%5.5f %5.6f\n",$1,$2/$4)}

5.2.12Multiplying Standard Stars Blackbody Spectrum

Once you know the spectral type of the standard star, you have knowledge about its temperature.
For instance, the intrinsic spectra of G-type stars can be approximated by a (single) blackbody
spectrum given by a temperature. This will work in many cases for B-, A-, and F-type stars, except
for the narrow lines emanating from atomic hydrogen in the stellar atmosphere. On the other hand,
you cannot use late-type stars, such as K- and M-type stars, for spectroscopic calibrations because
their spectra cannot be treated as blackbody radiation owing to the signiﬁcant absorption by thestellar atmosphere.
Here, we know that our standard star of HR 72 is categorized as G0V star, yielding a T of
5930 K (see Allen’s Astrophysical Quantities, fourth edition, p. 151). We therefore simply multiply
the blackbody spectrum at 5930 K.
It should be noticed that, when you use B- or/and A-type star(s) as a standard star, you will
have spurious (i.e., not real) emission lines at the rest-frame wavelength of Pfγ (λ = 3.74µm) andBrα (4.05µm) in the L band due to their deep absorptions. In this case, you have to estimate ﬂux
density of the star by interpolating the stellar spectrum; you have to supply an appropriate value
by interpolating from the adjacent (or neighbor) values.

43

gawk -f tools/bb5930 fIR00188_fHR72_1b8ccw.qdp > fIR00188_fHR72_1b8ccwbb5930.qdp

Now we have the desired spectrum. In the above example, we have multiplied a (constant) scale
factor to the spectrum. The content of bb5930 is as follows.

{printf("%1.6f %lf\n",$1,$2*1.1906*10^5/($1^5)/(exp(14388.3/$1/5930)-1))}

5.2.13Flux Calibration

The very last step of data reduction is ﬂux calibration. If there have been photometric observations
taken with a small aperture at the same band, you can obtain the ﬂux scale from those observations.
If there are no such observations, you will have to derive the brightness of the target on the basis of
an intensity (measured over the slit width) ratio between the target and the standard star. We canestimate the L band magnitude of HR 72 to be 5.00 mag from the V band magnitude of 6.46 mag
and the color V − L = 1.46 for G0V (see Allen’s Astrophysical Quantities, fourth edition, p. 151).
Since the L band magnitude was measured over a wavelength range of λ = 3.547±0.285 µm (Allen’s
Astrophysical Quantities, fourth edition, p150), we should compare the integrated intensity (along∫ 3.547+0.285µm
the wavelength axis; i.e., S = 3.547−0.285µm Bλ dλ ) of the target source over the wavelength range
with that of the standard star. In addition, we have to consider the total integration time whencomparing them13 . In practice, we have to consider the following points,
1. The signal from HR 72 is approximately 80 times brighter than that of IRAS 00188 whenmeasured over the entire L band range (λ = 3.262 − 3.832 µm). If you only extract the
spectral data within the above wavelength range, a program named normalize written by us(Imanishi) will calculate the integrated intensity.
2. The target data have a larger product of EXPITIME×NDR×COADD as compared to thestandard source by a factor of 18.75. See their FITS headers to conﬁrm this; the standard
star and the object, respectively, have EXPTIME = 1.0 and 2.5 sec, NDR = 2 and 6, and
The total integration time for IRCS data can be calculated from EXPITIME (exposure time per frame) × NDR(number of NDRs) × COADD (numbers of the exposure at each slit position) × the frame numbers you combined(§5.2.3). If you had combined multiple frames with averaging (i.e., comb = average for task imcombine), instead ofsummation (i.e., imsum task or comb = sum for task imcombine), the last frame number term should be 1.0.
1344

COADD = 10 and 15. Recall that we have combined 12 and 20 frames for the standard starand the target, respectively (see §5.2.3).
3. HR72 has L = 5.00 mag
Using the above information, we can calculate that IRAS 00188 should have integrated intensity
over the slit of L = 12.9 magnitude 14 . Since L =0 mag is deﬁned by the ﬂux density of Fλ = 6.59 ×
10−11 (W m−2 µm−1 ) (see Allen’s Astrophysical Quantities, fourth edition, p.150), the signal of L =12.9 mag corresponds to Fλ = 4.6 × 10−16 (W m−2 µm−1 ) at λ ∼ 3.55 µm. If you want to present
the data in the ﬂux density scale [Fλ = 10−16 (W m−2 µm−1 )], you can get it by scaling by a factor of
about 2.0 because the spectrum before the ﬂux calibration (fIR00188 fHR72 1b8ccwbb5930.qdp)
has a mean signal of 2.3 over the wavelength range λ = 3.547 ± 0.285 µm.

gawk -f tools/multi20E0 fIR00188_fHR72_1b8ccwbb5930.qdp >fIR00188_fHR72_final.qdp

The above gives you the fully calibrated and desired FINAL! spectrum. The multi20E0 ﬁlecontains argument for gawk as follows,{printf("%5.6f %5.6f\n",$1,$2*2.0)}

Finally, the following command shows you the ﬁnal spectrum (shown in Figure 14):

qdp fIR00188_fHR72_final.qdp/xs

14¿From the deﬁnition, if the signal is a factor of A times dimmer, its magnitude increases by a factor of 2.5log10 (A)
mag.

45

Figure 14: The fully calibrated ﬁnal spectrum.

qdp> hardcopy
qdp> q
lpr pgplot.ps

The above command prints the spectrum. If you want to have a more elegant plot, e.g., writinglabels, changing scales etc, you can play with the qdp parameters by consulting the manual.

qdp fIR00188_fHR72_final.qdp/xsqdp> la x Wavelength (\gmm)
qdp> la y F\d\gl\u (10\u-16\d W m\u-2\d \gmm\u-1\d)Notice that the above brightness is that measured in the slit of IRCS. If you observe a compact
object with IRCS, the image size at the L-band (. 0.′′ 5) is usually smaller than the slit width (0.′′ 9).
46

The good tracking accuracy of Subaru Telescope produces very small slit loss. Even if the slit loss
stays around the 10–20% level, the eﬀect cancels out each other by dividing the object with thestandard star. In this case, the resultant ﬂux calibration uncertainty would be smaller than that ineach spectrum (10–20%).
In the case of the short wavelength observations at λ . 2.5 µm, the background emission is lessintense, hence, one can use a fainter standard star. In addition, we can utilize 2MASS data, allowing
us to compare directly the magnitudes of the standard star at the J, H, and K bands. Unfortunately,
there is no such widely used catalog for the L band. We therefore have to extrapolate the desired
ﬂux of a standard star to the L band from short wavelengths ﬂux using the stellar colors. In general,bright standard stars, which are suitable as standard stars for L-band spectroscopy, are saturated
in the 2MASS measurements and accompanied by large uncertainties. We usually estimate the Lband ﬂux of a bright standard star from the optical V band ﬂux. In conclusion, the users should pay
special attention to the fact that the photometric error of standard stars at L band is signiﬁcantly
larger than those for optical and short-wavelength infrared measurements.
Regardless of these drawbacks, we believe that the above approach will give reasonable resultsin the sense that the uncertainty level is at most 20–30 %. This is because the resultant ﬂuxes
for compact sources (e.g., AGNs) estimated in this way have a fairly good consistency with the
published photometry results taken with small apertures.

47

Appendix
A

Error Estimates
The readers should bear in mind that the method to evaluate errors
described below is one of the many possible approaches.

By dividing the spectrum of the target by that of the standard star, we obtained the ﬁnal
spectrum, but the spectrum does not have error bars. Even without error bars, the readers of yourpaper can guess the uncertainty from the scatter of the data points at the line-free wavelengths.
However, there is no doubt about the usefulness of providing uncertainties of the ﬂux levels for theindividual spectral elements. We realize that the method to evaluate uncertainties certainly is notunique and depends on researchers. For instance, one may adopt Poisson noise of the total signal,
including both the target and background emission, as an uncertainty.
One of the approaches, we describe here, is to give uncertainties from the standard deviation ofseveral (independent) measurements toward an object. Here, all the measurements must be donewith exactly the identical observing parameters. For this purpose, we will split the whole data
sets of the target and standard star into a few sub-sets, and obtain one-dimensional spectra for
the individual sub-datasets. Subsequently we normalize each spectrum to equalize the ﬂux levels,and then estimate uncertainties based on the scatter of each spectral element among the sub-sets.
Finally, we should obtain a ﬁnal spectrum with errors through dividing target spectrum by that ofstandard star (here, we suppose that the both have error bars), multiplying the blackbody spectrumof the standard star, and ﬂux-calibrating.
In §5.2.2, we reduced only the ﬁrst-half of the entire dataset for both the standard star HR 72
and the target IRAS 00188. By analyzing the remaining data in the same fashion, we can estimate
the uncertainties. Speciﬁcally, we perform the following operations using the software prepared byourselves. We create the text ﬁles, fHR72 1.qdp and fHR72 2.qdp, from the former and the latter
half of the HR 72 data. We also make the text ﬁles, fIR00188 1.qdp and fIR00188 2.qdp, from the
IRAS00188 sub-datasets. We simply repeat the previous procedures accordingly up to wavelength
calibration (§5.2.8), and then bin 8 spectral elements. After completing all the steps on the second
half of the data, do the following:48

cd spec
tools/normalize fHR72_1ccwb8.qdp > fHR72_1normccwb8.qdp
tools/normalize fHR72_2ccwb8.qdp > fHR72_2normccwb8.qdp
paste fHR72_1normccw.qdp fHR72_2normccwb8.qdp > fHR72_12normccwb8.tmp

Above, we normalized the spectra of the former and the latter halves of the standard stars data,
based on the signal levels, then combined them. In general, the slit losses do not stay at a constantvalue between sub-datasets, and varies slightly. We believe that we do not need to consider suchvariation of the absolute value of the signal as it causes unnecessarily large (i.e., physically no-sense)
error bars. Since our primary interest lies in the spectral shape, above procedures would be justiﬁed.
However, we encourage the readers to address this issue from various aspects.
emacs fHR72_12normccwb8.tmp
We should delete the ﬁrst three rows because the program appends unnecessary information,such as pixel numbers, sum of signals, and the normalization constant, all of which are not required
for the latter stages.
Next is:gawk -f tools/sigma fHR72_12normccwb8.tmp > fHR72_12normccwb8.qdp

We combine the former and latter half datasets after the normalization step. At the same time,
the uncertainties (standard deviation) are also derived.
The content of the ﬁle sigma is as follows:{printf("%1.5f %lf %lf\n",$1,1/2*($2+$4),sqrt((($2-($2+$4)/2)^2+($4-($2+$4)/2)^2)/2))}

49

When combining three sub-datasets, the following calculation is required.

{printf("%1.5f %lf %lf\n",$1,1/3*($2+$4+$6),sqrt((($2-($2+$4+$6)/3)^2+($3-($2+$4+$6)/3)^2+($4-($2+$4+$6)/3)^2)/3))}

The above procedures are for the standard star data. For the target data, the same proceduresmust be applied accordingly. The readers should perform the following:tools/normalize fIR00188_1ccwb8.qdp > fIR00188_1normccwb8.qdp
tools/normalize fIR00188_2ccwb8.qdp > fIR00188_2normccwb8.qdp
paste fIR00188_1normccw.qdp fIR00188_2normccwb8.qdp > fIR00188_12normccwb8.tmp
emacs fIR00188_12normccwb8.tmp (delete the first three unnecessary rows)
gawk -f tools/sigma fIR00188_12normccwb8.tmp > fIR00188_12normccwb8.qdp

At this stage, we have obtained spectra with error bars for both the standard star and thetarget. Next, we should divide (1) the target spectrum by that of the standard star, (2) multiply
the blackbody spectrum corresponding to the standard stars temperature, and (3) apply ﬂuxcalibration. After these operations, we hope that you have obtained a FINAL spectrum where theuncertainties are included. The following steps should guide you toward that goal:

paste fIR00188_12normccw.qdp fHR72_12normccwb8.qdp > fIR00188_fHR72_12normccwb8.tmp

Here, we have combined the text ﬁles of the standard star and the target using the Unix command”paste”, then,50

gawk -f tools/spectrum fIR00188_fHR72_12normccwb8.tmp >fIR00188_fHR72_12normccwb8.qdp

The ”gawk” scripts consists of the following,{printf("%1.6f %lf %lf\n",$1,$2/$5,sqrt($3*$3/$5/$5+$6*$6*$2*$2/$5/$5/$5/$5))}

It should be noted that we calculated an error propagation for the divided spectrum because
both the standard star and target spectra have their own uncertainties. If the target data arerepresented by X ±σX and those for the √standard star as Y ±σY , then the propagated uncertainties( σ )2 ( σ ·X )2
X
associated with X/Y can be written as+ YY 2 . Then, we multiply with the blackbody
Y
spectrum.
gawk -f tools/bb5930a fIR00188_fHR72_12normccwb8.qdp >IR00188_fHR72_12normccwb8bb5930.qdp
The content of the ﬁle ”bb5930a” is as follows:{printf("%1.6f %lf %lf\n",$1,$2*1.1906*10^5/($1^5)/(exp(14388.3/$1/5930)-1),\\$3*1.1906*10^5/($1^5)/(exp(14388.3/$1/5930)-1))}

Notice that we made a minor correction to the previously-used gawk ﬁle ”bb5930” because thecurrent data has error bars consisting of three columns, thengawk -f tools/multi530E-2a fIR00188_fHR72_12normccwb8bb5930.qdp >fIR00188_fHR72_12normccwb8bb5930flux.qdp

51

We subsequently make the ﬂux-calibration to the target spectrum based on the calculation
mentioned in §5.2.13; the gawk ﬁle ”multi530E-2a” consists of{printf("%5.6f %5.6f %5.6f\n",$1,$2*0.05300,$3*0.0530)}

Of course, the normalization factor is diﬀerent for data with and data without error bars.
Insert read serr 2 to the ﬁrst line of the ﬁle opened by the editor,

emacs fIR00188_fHR72_12normccwb8bb5930flux.qdp &qdp fIR00188_fHR72_12normccwb8bb5930flux.qdp/xs

The last two commands display the spectrum with error bars (Figure 15).

Figure 15: Obtained spectrum of IRAS00188 with error bars.
Inspecting the 8 pixel-bin’ed spectrum in Figure 15, one probably realize that the scatter of thespectrum is large at the short wavelength region, which is caused by the descent Earth’s atmospherictransmission. We thus decided to increase the numbers of the binning pixels in the wavelength range.

52

As described before, the number of binning pixels should be optimized on the basis of thespectral quality of the target. Notice that we dealt with the target and the standard star data
separately before dividing the former by the latter when binned the spectral elements. In otherwords, we are not in favor of performing spectral binning on the resultant spectrum obtained by
dividing the (target) spectrum by the standard star spectrum. It is true that both methods willcertainly produce very similar results with respect to each other at wavelengths where the Earth’s
atmospheric transmission curve is smooth. However, in the wavelength range where the Earth’s
atmospheric transmission is highly dependent upon the wavelength (e.g., the short wavelength
range of the infrared L-band), the former method must be selected, as the latter procedure isinappropriate. We strongly encourage readers to address this issue.
Consequently, after searching for appropriate binning numbers, we obtained the desired spectrum
as shown in Figure 16. Have you obtained a similar spectrum? If so, you should be ready to write
the paper!

Figure 16: The ﬁnal spectrum of IRAS 00188, adopted from Imanishi et al. (2006 ApJ 637 114).

53