Subaru Data Reduction CookBook: COMICS— Version. 2.1.2e (November 21, 2012)—(e) HD200775

at 18.8 µm(f)18.8 µm PSF(g) HD200775(h).24.5 µm PSF

at 24.5 µm

Based on the textbook in Japanese by Y. Okamoto
Combined eﬀort of the past and current members of the COMICS team and Subaru Telescope
Current Version Editor: Ray S. Furuya

Contents
1 INTRODUCTION OF COMICS
1.1 Mid-infrared Astronomy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 An Overview of COMICS . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
3
4

2 COMICS Data Set
2.1 Chopping and Nodding, and Structure of Raw Data . . . . . . . . . . . . . . .
2.2 Summary of COMICS Data Reduction . . . . . . . . . . . . . . . . . . . . . .

55
93 FITS FILE AND HEADER
3.1 From the 1st to 8th Lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Header Information Common in All Subaru Data . . . . . . . . . . . . . . . .
3.2.1 Header Parameters for the Status of the Telescope . . . . . . . . . . . .
3.2.2 Header Parameters for Setting Optics . . . . . . . . . . . . . . . . . . .
3.3 Header for COMICS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14
14
15
15
16
17

4 Overview of Data Reduction: A Standard Procedure
4.1 Getting Started: the Reduction Software . . . . . . . . . . . . . . . . . . . . .
4.1.1 q series . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.2 IRAF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Which Data Should I Work On? . . . . . . . . . . . . . . . . . . . . . . . . . .

19
20
20
20
20

5 Data Reduction: The Details
5.1 Imaging at N and Q Bands . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.1 Making Dark Frames . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.2 Making Flat Images . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.3 Making a Bad-Pixel Mask . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.4 Reduction of Science Data . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.5 Adding the Images . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.6 Miscellaneous . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Analyzing Spectroscopy (N-band Low Resolution; NL) Data . . . . . . . . . .
5.2.1 Preparing for a Dark Frame . . . . . . . . . . . . . . . . . . . . . . . .
5.2.2 Making a Flat Frame . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.3 Reduction of the Read-out Pattern Noise . . . . . . . . . . . . . . . . .

23
23
23
24
26
27
28
36
38
38
38
40

0

5.2.4 Data Reduction of Scientiﬁc Images . . . . . . . . . . . . . . . . . . . .
5.2.5 Extracting the Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.6 Adding the Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.7 Wavelength Corrections . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.8 Correction of Atmospheric Absorption and Eﬃciency . . . . . . . . . .
5.2.9 Diﬀuse Emission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3 Frequently Asked Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 Contamination of the Ozone Line at 9.6 µm . . . . . . . . . . . . . . .
5.3.2 Should I perform nodding for my observations? . . . . . . . . . . . . .
5.3.3 I have a ”ghost” in my low-resolution image at the N band, what can I
do? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.4 I get errors when I transform my image with geotran . . . . . . . . . .
5.3.5 geomap relevant questions . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.6 It takes a looooooooooooooooong time for task q startrace to complete... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41
48
49
49
49
51
52
52
52
52
53
53
53

Revision History
2001/05/24
2004/11/13
2007/09/24
2007/10/09

Ver. 0.0j
Ver. 1.1β-jVer. 2.0j
Ver. 2.1j

2008/04/01 Ver.2.1.0e
2008/05/15 Ver.2.1.1e
2012/11/21 Ver.2.1.2e

Yoshiko K. Okamoto
Mitsuhiko Honda /Yoko Okada
Takashi Miyata
Takafumi Kamizuka,
Takashi MIyata, and
Ray S. Furuya
Ray S. Furuya, discussed w/
T.M., T.U., & T.K., et al.
Ray S. Furuya
R.S.F. & C.C.P.

The original version
Revision
Revised for Subaru Data Reduction School
Overall revision, particularly for § 5
Deleted description about s series as
it has been merged with q series
English version
Corrected errors
Language corrections

Chapter 1
INTRODUCTION OF COMICS
1.1

Mid-infrared Astronomy

The mid-infrared (MIR) regime spans between ∼ 5 and 30 µm in wavelength (λ); there
exist two windows of the atmosphere surrounding Earth in this MIR regime. These arereferred to as the N and Q bands, which have been long-unexplored because of limitations
imposed by high background noise and technical diﬃculties in the MIR detectors. However,this situation has changed dramatically due to recent progress in the development of nextgeneration 2-dimensional (2D) array detectors, as well as the construction of large-diameter
optical/infrared telescopes. In order to make MIR observations possible at Subaru Telescope,
Cooled Mid-Infrared Camera and Spectrometer (COMICS) has been developed as one of theﬁrst-generation instruments by the team led by Professors H. Kataza and T. Yamashita.
The MIR wavelength regime stands unique due to a number of special characteristics.
First, these longer wavelengths can penetrate interstellar clouds than the shorter ones, e.g.,
optical/near infrared (NIR) do. This allows us to perform deeper searches for deeply embedded
and newly found objects identiﬁed in optical wavelengths, although MIR emissions suﬀer fromabsorptions due to interstellar silicates. Second, the MIR regime is the most eﬃcient band to
detect blackbody (thermal) radiation at temperatures of 100 – 400 K, making COMICS oneof the best instruments to observe warm thermal emission from circumstellar disks aroundyoung stellar objects and warm gas and dust around late-type stars. Third, the MIR regime
is a key band to trace the role and evolution of interstellar dust because there are numerous
emission and absorption features spread across it. The best-known of these features can beseen around 9.7 µm and 18 µm, and are believed to be of silicate origin and detected overrather wide ranges in the two atmospheric windows. In addition, it is known that there exist
unidentiﬁed IR (UIR) bands peaked at λ = 6.2, 7.7, 8.65, 11.25, and 12.8 µm, which can bedetected towards various type of objects.

3

1.2

An Overview of COMICSCOMICS is an instrument at the 8-meter Subaru Telescope which is also optimized to observe
MIR wavelengths. It is expected to explore various astronomical objects such as those in thesolar system, newly formed stars and planets, interstellar dust, as well as nearby and even distant galaxies. COMICS oﬀers two modes of observations: imaging and long-slit spectroscopy.
COMICS is mounted at the Cassegrain focus to minimize thermal emission from the telescope.
In general, one must cool down the instruments for the MIR bands since the black body radiation peaks around 10 µm at room temperature. In fact, all the optics and the associated parts(such as motors) are stored in the dewer and cooled down to less than 30 K. In addition, youneed to cool the detectors (mentioned below) further, down to 4 K. This reduces dark currents.
In order to have diﬀraction-limited image with the 8 m telescope, COMICS is designed witha pixel ﬁeld of ∼ 0.13arcsec/pix for imaging mode. For spectroscopic observations at 10 µm (Nband), COMICS oﬀers low (L)-, medium (M)-, and high (H) resolution modes with R ∼250,
2500, and 10000, which are designated as NL, NM, and NH, respectively. At 20 µm (Q band),COMICS provides medium- and high-resolution modes with R ∼250 (QM) and 5300 (QH),
respectively. COMICS equips four gratings for all the modes except for the high-resolution
mode at Q-band. We adopted Si:As 320×240 IBC array detector which is sensitive between
8 and 28 µm, provided by the Raytheon. The IBC type detectors has a fairly high quantum
eﬃciency which enables us to improve detection limits signiﬁcantly. In the imaging mode,
COMICS uses one of the IBC-type detector. Users can obtain a relationship between the slit
and the observing target(s) with a high accuracy as the detector at the imaging side works
as the slit viewer. For spectroscopic observations, we conﬁgure all ﬁve detectors in parallel toachieve high observing eﬃciency. However, one must slightly rotate the grating to detect photons that fall into the gap between the detectors. Most of the spectroscopic observations haveconﬁgured two sets of such rotation angles to detect all the photons without gaps along the
wavelength, which makes it possible to observe the entire atmospheric window in the medium
resolution mode. In addition, we oﬀer a total of four slits — available in both the 10 and 20µm bands — having diﬀerent widths for possible variation in seeing. Table 1.1 summarizes
the key parameters of the instrument.
Pixel Field of View
Detector
Field of View
Wavelengths

Spectral Resolution

Imaging
Spectroscopy′′
0.130 /pix
0.165′′ /pix
320x240 Si:As ×1
320x240 Si:As ×5′′′′
42 × 3140′′ for long-slit
N:8.8–12.4 µm (5 bands)
N: 7.8–13.3 µm
Narrow bands for UIR 8.6, 11.24, [ArIII], [SIV], and [NeII]
Q:17.7–24.5 µm (4 bands)
Q: 16.7–20.7 µm
R ∼10–20
N band: R ∼250, 2500,10000
R ∼60–70 (narrow band ﬁlters)
Q band: R ∼2500, 5300

Table 1.1: An Overview of COMICS

Chapter 2
COMICS Data Set
This chapter provides an introduction to COMICS data and FITS ﬁle structure, and gives arough idea of how to reduce observing data.

2.1

Chopping and Nodding, and Structure of Raw Data

For most of observations with COMICS, high-speed read-out is used to compensate for thehighly time-variable background level in the MIR bands. For this purpose, we employ ”chopping” of the secondary mirror and ”nodding” of the primary mirror. Most of the atmosphericbackground variations can be removed by the chopping method. However, since chopping isbeing done through tilting of the secondary, a subtracted image obtained (by taking one image
from another) leaves a low-level residual pattern in many cases. This is because optics for thetwo chopped beams are slightly diﬀerent from each other 1 . Such a residual pattern can readilydecrease detection limits, which severely aﬀects the ability to detect faint objects. To avoid
this eﬀect, it is common to ”nod” the telescope so that the source appears alternately in thesignal and reference beams. Nodding scans observe an oﬀset position with a chopping sequencethat is the same as for the target position. This will cause a similar residual pattern to thatseen in the target position. In the data reduction process, we subtract the oﬀset images(s)to remove this low-level residual pattern. Given this purpose, it is suﬃcient to ”nod” witha lower throw-frequency than ”chopping” does. A nodding cycle is typically several minutes
to obtain a pair of images. Lastly, one should keep in mind that the decision to use nodding
depends on the source brightness and/or scientiﬁc goals.
Figures 2.1 and 2.2 illustrate the rough idea of how COMICS obtains a set of images andshows structures of the image FITS ﬁles, respectively. A standard observation will perform
chopping and nodding scans towards slightly diﬀerent pointings so that the target will fall
onto diﬀerent pixel positions of the images (see Figure 2.1). Since chopping will be done ateach nodding beam, the user should obtain a set of data consisting of the four images labeled
1

Although telescope eﬃciency is suggested to be responsible for causing such a residual pattern, amplitudelevel of the residual pattern may be ignored at Subaru Telescope depending on the scientiﬁc cases.

5

A to D (see Figure 2.1). The direction, as well as throws for chopping and nodding, must be
selected by considering not only source structure but also its surrounding environment. For
instance, the user may select a chopping throw of ∼ 10′′ toward an arbitrary direction2 in thecase of a point-like source located in a region without diﬀuse emission. On the other hand, theuser should carefully select these parameters in the case of extended sources and/or deeply
embedded object(s) in diﬀuse emission. The user should take into account several constraints
such as source size, clumpiness of the surrounding clouds, and so on. For many cases, such a
selection gives reasonable atmospheric cancellation and good mechanical performance.
As can be seen in Figure 2.2, each of the four images (labeled A, B, C and D) contains
multiple images. Assuming that the user has repeated n times of chopping scans (i.e., ChopNum= 2n), the image data taken with beam A should have the number of n images with A1 ∼ An .
Each of the A1 , A2 ,...., and An th images contains either a summation image over the mtimes of read-outs (= exposures), i.e., Σi Aij (ADD mode), or sub-images without making thesummation image (RAW mode). In the latter case, Aj th image should include the numberof m raw images of A1j , A2j ,A3j ,...., and Am
j . Another set of the chopping data, i.e., the beam
C images, has the identical data structure. Each nodding image has a substructure withthe two chopped beam images. The nodding images from beams A and C constitute a pair
together with those from the nodding beams of B:B1 ∼ Bn and D:D1 ∼ Dn . In conclusion, astandard set of COMICS data that are ready to reduce contain the other nodded beam dataof B:B1 ∼ Bn and D:D1 ∼ Dn .
We deﬁne names of images in each substructure as follows.
Name of image
nodding beam
chopping beam
beam
ﬁle
frame
exposure

Description
Nodding images having an identical pointing center
The same as above, but for chopping
Taken at diﬀerent positions due to chopping and nodding
Each FITS ﬁle for the obtained image
An image set taken within a chop beam
Each readout

Notation
All the images in the A and C
All the images in A
A, B, C, D
pairs of A and C, B and D
Subscript of A
Superscript of A

In practice, COMICS creates three types of data sets.
1) RAW mode: keeps all the exposure and all the frames. Images in each ﬁle are Aji andCij .
2) ADD mode: adds all the exposures within a frame (i.e., collapses all the images into
an image), but keeps all the frames. This is the data set for a typical COMICS observation.
Images in each ﬁle are Σj Aji and Σj Cij .

2

We suggest that the user select a chop direction that keeps the secondary mirror at the most stable
position.

D

C
west

chop

north

south

east
B

A
nod

A

B

C

D

Figure 2.1: An illustration of a typical data set taken with COMICS. The user will obtain adata set consisting of four diﬀerent images labeled A – D. These 4 images have been taken atslightly diﬀerent pointing centers in the plane of the sky through chopping and nodding scans.
The user can select an arbitrary direction for secondary mirror chopping and primarily mirror
nodding. Note that each of the images contains multiple images, as shown in Figure 2.2.

S

u

m

m

a

t

i

o

n

o

f

e

a

c

h

e

x

p

o

s

u

r

e

o

r

r

a

e

x

p

o

s

u

r

e

d

a

t

a

w

S

i

t

c

h

i

n

g

c

h

o

p

p

i

n

g

b

e

a

m

s

w

C

h

o

p

(

p

i

n

t

n

g

i

s

m

c

e

a

s

n

s

o

N

d

d

)

i

t

i

h

n

t

g

h

b

e

e

s

a

a

m

m

t

e

a

s

k

t

e

r

s

u

d

c

t

a

u

t

r

a

e

w

E

a

c

h

n

c

o

o

d

n

d

t

i

a

n

i

g

n

c

i

n

r

g

e

a

t

t

h

o

e

s

s

a

e

s

f

o

e

r

t

4

o

f

p

d

o

a

i

n

t

t

a

f

i

l

e

s

Figure 2.2: A schematic image of a typical a COMICS data ﬁle. The chopped and nodded
beam images set illustrated in Figure 2.1 individually contain sub-images. This is because
almost all the data will be taken through multiple chopping scans for each nodding beam.
Each chopping beam image contains either the further multiple sub-images, (i.e., the rawdata) or a summed image of them.

2) ECO mode: Since this mode collapses all the images into an image on-line, and makes afurther summation for all the frames at each beam, the user will obtain a map for the chopping
and nodding at each. Images in each ﬁle are Σi,j Aji and Σi,j Cij .
In practice, some observers may have acquired images without chopping and nodding for
obtaining dome-ﬂat data (see §5.2.3). In this case, Chop=0 appears in their headers. Users
should see multiple frames whose numbers have been speciﬁed with the parameter of ChopNum.
Notice that the number of exposures has been repeated with that speciﬁed in the individual
frames as well. As explained above, all the exposure frames are kept for CoAdd = 0, while thenumber of ChopNum images are saved with CoAdd = 1 because all the exposure images in eachframe are summed. Lastly, the user obtains two images for CoAdd = 2 because all the images
have been summed for each chopped beam.

2.2

Summary of COMICS Data Reduction

This section describes an overview of the COMICS data reduction to obtain fully calibrated
images. The details of each step of data reduction will be described in §5.
1. Subtraction of Chopping and Nodding Data We assume that the users has obtained
a non-calibrated (i.e., “raw”) data set containing four images (§2.1) where the object of interest
can be seen at diﬀerent pixel positions. The user begins by subtracting chopped images taken
sequentially with the diﬀerent chopped beams (note that they must be taken with the same
nodding beam). We refer to it as chopping subtraction. Namely, as shown in Figure 2.2, theuser needs to perform the subtraction as follows,(A1 − C1 ) + (A2 − C2 ) + · · · · · + (An − Cn ) = α(B1 − D1 ) + (B2 − D2 ) + · · · · · + (Bn − Dn ) = β
which gives two images after the chopping subtraction (see Figure 2.3). This operation usually
makes it possible to subtract almost all the background pattern as well as the residual one
originating from the dark current. However, the user may ﬁnd a very low-level residual patternwhich is caused by the use of diﬀerent optics. To eliminate such a pattern, the user would needto subtract an image generated by the chopping subtraction (hereafter, nodding subtraction),
α−β =Ito obtain a nodding-subtracted image I. For many cases, the above operation will give you
images with sky levels close to zero. If you have made multiple nodding observations, youneed to make animage, I, for each nodding.

Chopping Subtraction

Residual Background
Pattern after Chopping
Subtraction

Nodding Subtraction

Figure 2.3: Subtraction of chopped and nodded images: The ﬁrst step in the COMICS data
reduction process is subtracting images between the same nodding beams and between thediﬀerent chopped beams; the resulting four images are stored in a data set. The subtraction
between the images A and C in the upper panel gives the image α in the middle (the samefor the B and D to make the subtracted image β). Subsequently, the user needs to subtract
the two images created by the chopping subtraction to make the desired image shown in thebottom panel (α, β →I); this operation will eliminate residual patterns caused by the useof slightly diﬀerent optics. Notice that the resultant image of I should contain the positive(white) and the negative (black) portion for the source of interest as a result of the subtraction
processes shown above.

2. Flat Fielding The next step is to divide each of the chopping- and nodding-subtracted
images by a ﬂat image to correct any inequity in the gain of the detector pixels. Flat images
are usually made through sky-ﬂat and dome-ﬂat observations. The former is used for imaging
and the latter for spectroscopy. For most COMICS observations, dome-ﬂat data are taken byobserving thermal emission from the dome, windscreen, or mirror covers etc, instead of fromobserving lamps.
3. Image Transformation (Distortion Correction) and Wavelength Calibration forSpectroscopy Observations For spectroscopic observations, the user will need to correct
distortion of the spectrum that shows a curved feature in the images because COMICS guides
the incident wave obliquely with respect to the grating. Figure 2.4 illustrates a standard
procedure to correct such a distortion. First, one has to derive a ”spatial axis” from a pointsource image such as the standard stars. Here, a ”spatial axis” is the axis free from distortionin the fully calibrated image (see Figure 2.4). On the other hand, one should make a spectrum
image for the sky that is created from, for example, summation or/and mean images using all
the chopped beams for the common nodding beam. Keep in mind that all the images used
for calibration must be made with the images found in the data-set ﬁle that contains scientiﬁcimages(s) to be calibrated. The second step is to derive ”dispersion axis” from the sky spectral
image, then transform the image so that the spatial and dispersion axes are perpendicular toeach other. Finally, it is important to calibrate the dispersion axis using a comparison with astandard model of terrestrial atmospheric emission.
4. Adding the Calibrated Images The next step is to add the calibrated images toincrease the signal-to-noise ratio (S/N). This operation must be done only for the images taken
at the same image center. (Check the telescope pointing accuracy during your observations
as well!) In some cases, in order to reduce the number of FITS ﬁles, users may want to
add certain numbers of images taken at the same position with a long observing time. (thisrequirement comes from a practical limitation of image numbers that can be stored in a ﬁle).
Finally, we suggest that users double-check whether or not the target source has shifted by
one or more pixels, even for images taken at the same position.
5. Flux Calibration with Standard Star(s) The ﬁnal step is to convert unit of theimage, namely, from count into brightness. It is not easy to ﬁnd (a) well-test catalog(s) ofstandard stars at MIR wavelengths. To our knowledge, the catalog by Cohen and colleagues(Cohen et al. 1999, AJ 117, 1864) deals with a rather large number of standard stars basedon the IRAS and LRS observations. In addition, the user may use some catalogs available
online, such as the UKIRT list. However, that lists only about 20 bright stars.
Cohen et al. (1999) listed ﬂux densities at N and Q bands measured at low-resolution. To
apply their results to your data, you must properly handle diﬀerences in the slit eﬃciencies
between theirs and that of COMICS. Notice that this eﬀect is particularly important when
you analyze diﬀuse emission because COMICS employs long-slit spectroscopy. For photometry
observations, the user is required to calculate ﬂux of the desired standard sources by integrating

- Making summation or Average of A1 〜An, C1〜Cn- Subtractioning dark current image

Making a Spectral Image of the Star
using Images of A, B, C, and D

Sky Spectrum ImageDistorted
Spatial
Axis
Distorted Dispersion Axis
Intensity

Intensity

Observing various lines from astronomical objects
and atmospheric emission to calibrate wavelength

The desired dispersion axis

Obtaining peak positions by applying a
Gaussian to calibrate the spatial axis

Spatial Axis

Figure 2.4: Determination of wavelength and spatial axes in spectroscopic observations. The
absolute coordinates of the wavelength axis can be identiﬁed from a comparison of the spectral
image obtained from a summation of the raw data with an atmospheric transmission model.
As for the spatial constant axis, the users should determine the axis by applying a Gaussian
proﬁle to the obtained spatial cutting line which passes the point-source along the wavelengthaxis.

the spectra shown in the paper for each N band ﬁlter. This is because almost no systematic
measurements have been done for individual ﬁlters, due to the very high background noise.
Finally, our experience strongly suggests that the user should pay attention not only todiﬀerences in airmass but also the presence or absence of clouds in the sky because they heavily
aﬀect the quality of the images. Given that sky noise level is intrinsically high at MIR bands
regardless of the clouds, some observers may continue to observe with a small amount of cloud
cover. However, it should be noted that even small amount of clouds can easily degrade yourimage quality. In many cases, the user is required to correct the eﬀect of the atmosphericabsorption even for the standard star data obtained with the same airmass. In conclusion, we
strongly suggest making sure that the spectral shape of the standard star(s) has a fairly goodconsistency with that of the atmosphere during the given observing time.

Chapter 3
FITS FILE AND HEADERYour COMICS data are stored in FITS format. COMICS FITS data usually have two types:
one is the series of image ﬁle names starting with ”COMA”, and the other begins with
”COMQ”. For both ﬁles, 8-digit images numbers follow the strings of COMA or COMQ.
The COMA ﬁles are images taken with CoAdd mode, as explained in §2, while COMQ are
those of the online (partially) reduced images using COMA ﬁles. Here, the COMQ data arecreated by collapsing all the images into an image taken with the same chopping-beam image,and subtracting another collapsed chopping-image created in the same fashion. In other word,
COMQ images are the summation images obtained by adding on-beam images and oﬀ-beam
images that are multiplied with −1.
The discussion below gives an introduction on how to read FITS header information ofyour COMICS images. You may see a header with a diﬀerent version that can be foundwith the Q GETVER parameter at the bottom. As for the details of your FITS header andcorrespondence with diﬀerent versions, please ask your support scientist.

3.1

From the 1st to 8th LinesSIMPLEBITPIXNAXISNAXIS1NAXIS2NAXIS3NAXIS4EXTEND

=
=
=
=
=
=
=
=

T / Standard FITS format
32 / # of bits per pixel
4 / of axis in frame
320 / # of pixels/row
240 / # of pixels/column
50 / # of pixels/time
2 / # of pixels/detector
F / ASCII Extension Table

The ﬁrst eight lines in COMICS header follow the FITS format standard deﬁned by the
IAU. The second parameter — BITPIX — shows bit number; BITPIX with minus (−) indicatesthe numbers of ﬂoating decimal points, while one without minus sign indicates the integer,

14which is the case for the current COMA and COMQ images. You will see BITPIX= -64 for
all the FITS ﬁles that are created by the COMICS data reduction package of the q series1 .
The third parameter, NAXIS = 4, indicates the dimension number of the image ﬁles. Notethat COMICS data have a four-dimensional (4D) structure; NAXIS1 and NAXIS2, respectively,
represent the x- and y axes of the detector, NAXIS3 does time-axis (NAXIS3), and NAXIS4
indicates the detectors. This means you obtain one single ﬁle even if you observed with morethan one spectrometer. If you observed using the imaging or low-resolution spectroscopy
mode, you see NAXIS4 = 1 for NAXIS = 4.
Since COMICS has a detector array consisting of 320×240 pixels, you should see NAXIS1=320
and NAXIS2=240. The value found in NAXIS3 would depend on the data format, i.e., clock
settings in COMICS jargon, you gave during the observations. NAXIS4 = 1 should be stored
for imaging mode, while a diﬀerent value other than 1 would be stored for each spectroscopymode.
If you have taken both the imaging and spectroscopy data at the same time, these data
are stored in two paired FITS ﬁles that have adjacent ﬁle numbers; a pair of odd- and evennumbered FITS images correspond to the imaging and spectroscopy data, respectively. In
most of cases, the ﬁle number for the imaging must be smaller by one than that for the
spectroscopy by one. If you observed only with imaging mode, the paired spectroscopy FITS
ﬁle has only the header part.

3.2

Header Information Common in All Subaru Data

The header information after
COMMENT = ’ +++++++++++++++++++ SUBARU COMMON’in the COMICS header parameter is common for all the data taken with the Subaru
Telescope where some basic parameters such as observing time, coordinates, conﬁgurations
for COMICS can be seen. Those who need to have the more detailed information for thecommon part, please see ”FITS no Tebiki” ver. 3.2. (in Japanese) and so on.

3.2.1

Header Parameters for the Status of the Telescope

We advise users to check the following parameters particularly for COMICS data.
INSROT =INST-PA =AUTOGUID= ’ON’
M2-TYPE = ’CS_IR’
M2-TIP = ’CHOPPING’
1-72.158 / Instrument Rotator angle (deg)9.566 / Instrument Rotator P.A. (deg)/ Auto Guider on/off/ Type of the Secondary Mirror (Opt/IR)/ 2nd Mirror tip-tilt on/off

The only exception is the output ﬁles from task q 2val, which quantizes with 2-levels, where the resultant
variable type becomes integer.

M2-ANG1 =M2-ANG2 =M2-POS1 =M2-POS2 =-0.978 / 2nd Mirror Angle-1
0.687 / 2nd Mirror Angle-2-1.104 / 2nd Mirror Position-1-1.571 / 2nd Mirror Position-2

The ﬁrst two lines give information about the ”instrumental rotator ”; the NST-PA parameter shows the rotator angle with respect to the source of interest. AUTOGUID tells you whether
the observations have been done with the auto-guider, which may have used for most of theobservations. The nine parameters starting with M2- give information about the secondary
mirror (however, many users will not need to give these special attention.)

3.2.2

Header Parameters for Setting Optics

You will see the following parameters that describe coordinate information of the images. In
many cases, we believe users cannot execute a command to overlay all the images simply byreferring to the coordinate information described. This is because nominal tracking accuracyof the telescope is not accurate enough to allow such blind image adding. Even if this is not
the case, you cannot add your images solely with the header information. For example, if
you adopted a chopping throw of more than 30”, imaged ﬁelds in the plane of the sky do notagree with those expected from the telescope pointing center because of the large throw of thesecondary mirror.
The header after
COMMENT = ’ ------------------------- COMICS Optics’
are those for COMICS optics.
OBS-MOD = ’spectroscopy’/ Observation ModeFILTER01= ’F01C10.50W6.00’/ Filter name/ID (pre-opt filter-1)
FILTER02= ’H21’
/ Filter name/ID (pre-opt filter-2)
FILTER03= ’F08C11.60W1.10’/ Filter name/ID (img-opt filter)
FILTER04= ’L01L10I’/ Lens name/ID (img-opt)
DISPERSR= ’G01L10L’/ Identifier of the disperser used
SLIT
= ’S02W160’/ Identifier of the slit
SLT-LEN =39.600 / Length of the slit used
SLT-PA =0.0 / Slit Position Angle (degree)
SLT-WID =0.330 / Width of the slit used
SLTCPIX1=
120.0 / Slit center projected on detector(pix)
SLTCPIX2=
160.0 / Slit center projected on detector(pix)OBS-MOD indicates the observing mode of either ”imaging” or ”spectroscopy”. Notice that
OBS-MOD becomes ”imaging” if you observed with the slit viewer.
FILTER01, 02, 03 and 04 are the ﬁlter names for fore-optics 1 and 2, the ﬁlter name for
imaging, and the name of the lens, respectively. The ﬁrst character indicates either ﬁlter (F),lens (L) or grating (G). The next two numbers are IDs for each optical component followedby the strings indicating type of the component. For instance, C.... indicates the center
wavelengths in µm for F, while H?? means ’opened’ has been selected with the ﬁlter wheel(the two numbers with H show the ID for the holes). L10 shows lens ID for 10 µm with the ﬂag
of I indicating Imaging. Similarly, L20 shows lens ID for 20 µm with the ﬂag of P indicating
Photometry. L10 (L20) in the grating mode with G has the same meaning as above, except
for the last string of L (Low), M(Medium) and H (High) resolution (i.e., dispersion) modes,
respectively.
SLIT shows the name of the slit with the two-digit numbers appended for the slit ID and
W.... for slit width in µm. Notice that the slit width of 160 µm corresponds to two pixels of0.33” for spectroscopy observations; these values can be found in SLT-WID. S??W000 indicatessole the mirror without use of any slits. If this is found, the data were taken with the imaging
mode without using any conﬁgurations for spectroscopy. S??D... indicate diameter of the
pin-hole. S?SLT-LEN shows length of the slit, which should be constant if there is no vignettingin the image. SLTCPIX1 and SLTCPIX2 can be ignored because no meaningful parameters havebeen stored for COMICS observations.
The parameters of WAV-MIN and WAV-MAX store some values only for N-band spectroscopy
observations using the current observing system. It should be noted that EXPTIME indicates
not integration time for all the observations, but that for each exposure. DET-TMP shows themeasured temperature of the mount accommodating the detectors. Last, unfortunately, theconversion factor of GAIN does not help you at all with the current system.
COMMENT = ’--- Spectroscopy only’
DISPAXIS=
1 / Dispersion Axis in frame
WAV-MIN =7500.0000 / Shortest wavelen (nm)WAV-MAX =13500.0000 / Longest wavelen (nm)WAVELEN =10500.0000 / Central wavelen (nm)COMMENT = ’ ------------------------- COMICS Detector’
EXPTIME =0.301 / 1 exposure integration time per exp(sec)
DET-TMP =7.65 / Detector temperature (K)GAIN=
350.000 / AD conversion factor (electron/ADU)

3.3

Header for COMICS

The following header lines starting ”Q ” after the line of
COMMENT = ’ +++++++++++++++++++ COMICS ORIGINAL’
are those for COMICS in particular. All the important information needed for your data
reduction such as those relevant to the data type (i.e., clock parameters), observing mode,
window function. In the below, we show the clock parameters that are especially helpful foryour data reduction.

Q_DETST = ’100110’/ Detector Readout Status
Q_CLKFL = ’/home/comics/cbin/clk/clkgen/012/c050.00150.001.00’ / Macro File
Q_PIXTIM=
150 / Clock duration for a pixel (0.1us)
Q_RRSTRT=
1 / Reset Row Start Width (ND)
Q_CHWB =3 / Wipe Exporsure Number in a Chop-beam
Q_CHEB =2 / Exporsure Number in a Chop-beam
Q_CHCN =50 / Chopping Number in this file
Q_CHAM =1 / Add Mode 0:RAW 1:ADD 2:ECO
Q_CHOP =1 / Chopping ON=1 OFF=0
Q_CTYPE =0 / Clock Type 0-9
Q_YSTRT =1 / Readout Region Y start
Q_1EXP =0.301 / Integration time per exp. (sec) = EXPTIME
Q_1FRAME=
0.603 / Integration time per frame(co-added) (sec)
Q DETST shows from which detector the data are being read from among the six detectors(i.e., one for imaging and ﬁve for spectroscopy); the parameter should be 1 for the detector
being read out while 0 for the remaining ones in not use. For instance, the parameter should be
100000 for imaging observations, 100110 for low-resolution spectroscopy at the N-band, and111111 for medium-resolution spectroscopy in the Q-band. The user should use caution in thecase of using a part of spectroscopy detectors such as the case for low-resolution spectroscopy
at N-band. In this particular case of 100110, the 4th dimension of the FITS ﬁles becomes(1,2), although the designation number of the detectors are 3 and 4 (one of the detectors is
for reference).
Q CLKFL shows s name for data type (i.e., clock parameter). The other important parameters to know for the clock status are Q PIXTIM (read-our rate), Q RRSTRT (electrical optimized
ND value), Q CHWB (number of wipes at each frame), Q CHEB (numbers of read-out that arerecognized as usable data at each frame), Q CHCN (numbers of chopping ×2), Q CHAM (CoAddmode), Q CHOP (whether or not chopping is used), Q YSTRT (the starting read-out number of”Y” in the case of a partial read-out), Q 1EXP (integration time per exposure in sec), andQ 1FRAME (integration time per frame).
The ﬁnal part of the header shows the status of chopping as follows,
COMMENT = ’ ------------------------- CHOPPING’
Q_CHTHRW=
60.00 / Chopping Throw
Q_CHDEG =53.34 / Chopping Degree
COMMENT = ’ ------------------------- FITS VERSION’
Q_GETVER= ’4.11’ / FITS header VERSIONEND
Q CHTHRW indicates chopping throw deﬁned by peak-to-peak in arcsec, and Q CHDEG shows
the position angle for the chopping direction in degrees.
Finally Q GETVER shows version of FITS header as described in the beginning of this section.

Chapter 4
Overview of Data Reduction: A
Standard Procedure
In this chapter, we describe a standard procedure to reduce COMICS data. Since we willdeal with the data with observing wavelengths more than one order of magnitude larger thanthe other IR/optical wavelength observations, users are strongly recommended to use theoriginally developed data reduction package by us. These are the q series data reduction
package (hereafter, q series) which are written in C programs. The user can use these
programs ”as is” and/or can modify them.
q series can be categorized into two groups. One contains the basic tools that can beused for a ”quick look” while at the summit. These are mathematical operations on images, extracting a part of the data along column, row, or time, and statistical calculations
for the desired portion of the image. The other group of tasks should be used to prepare
for publication-quality ”ﬁnal” images, e.g., listing header information (q headlst2), determination of spatial axis (q startrace), and wavelength calibration (q sky nlow) and so on.
q series is specialized to reduce COMICS data, but some operations such as transformationand shifting of images are better achieved using other packages such as IRAF. We believe thatboth q series and IRAF can handle almost all the data currently being taken.
To get to the on-line help for q series commands, type the command name from thecommand line of your terminal (currently this is written in Japanese).
In the following, we show the standard procedure of data reduction with the command
names in q series. Each of the commands is described in the subsequent chapter in detail.
We suggest preparing for a summary of your FITS ﬁles that should be reduced with describing
data types, qualities and so on. Such a summary table would be more useful than a list createdby reading FITS headers of the target ﬁles. This is because most users will have a huge numberof FITS ﬁles.

19

4.1

Getting Started: the Reduction Software

We assume that the user will reduce data either on UNIX or LINUX. As described above, youneed to have
• q series• IRAF
Furthermore, the user should have a basic knowledge of script languages such as awk and perl.

4.1.1

q series

Users can obtain the latest version (version 30 as of 2007 October) of q series fromhttp://canadia.ir.isas.ac.jp/comics/open/rbin/rbin.html
After getting the source codes contained in rbin.tgz, extract the ﬁle. You will see a
directory called rbin.030. Please refer to the README ﬁle under 1 README of the directoryfor further information of installation. In this manual, we describe all the commands in theq series without a path. So you should tell your shell the proper path for the commands.

4.1.2

IRAF

Please refer to the other documents or/and your local experts for installing IRAF.

4.2

Which Data Should I Work On?

Before starting data reduction, the user should properly associate data types (e.g., target,
standard star, dark etc) with FITS ﬁle numbers. For this purpose, we suggest listing thefollowing information that will help your reduction. In the list, of course, you should exclude
FITS numbers for unnecessary data and the overhead of the observations (e.g., test exposure,etc.):• Object name• Image number• Data format (i.e., clock type)• Type of observations (Photometry or/and Spectroscopy), and ﬁlter names• Status of the telescope (e.g., Azimuth/Zenis Angles) and COMICS (e.g., InstPA)• Chopping parameters

Table 4.1: Summary of a Standard Procedure: Imaging
0. Data Examination
1. Making Dark

Checking clock and optics
Dark per exposure

2. Making ﬂats for imaging(sky ﬂat)

Subtracting dark from raw sky data
Splitting images for each chopping beam
Statistical operations between
frames and between ﬁlesSmoothing and normalization
3. Making a mask
Finding bad pixels
4. The Target Images
Subtracting chopping data
Analysis of Standard Star(s) Removing pattern noises
Dividing with ﬂat
Shifting & Adding
Measurement of counts for standard starq headlst2
q list statq fcombine
q arithq arithq bsepq list stat
IRAF:gauss, q arith
q badpix
COMQ images
q subchq arithq arith, q fcombine
IRAF:imshiftq photoor IRAF:imexamine

Deriving conversion factor from count to ﬂux density

If you have diﬃculty in making such a list for some reason or/and are analyzing data
retrieved from the archive system, you will probably need to extract such information fromthe FITS header of each ﬁle. Since observers may have many FITS ﬁles, it’s easy to use thetask q headlst2 in q series.
A very useful feature of q headlst2 is that it can read an ASCII text ﬁles containing a listof FITS ﬁle names. For example, assume that your data are located under /data1/20000720A,
type the following commands at the directory to make such an input list ﬁle as follows,
ls /data1/20000720A/COMA* | awk ’{print substr($1,1,16),substr($1,18,12)}’ > 720.lista
ls /data1/20000720A/COMQ* | awk ’{print substr($1,1,16),substr($1,18,12)}’ > 720.listq

Of course, you can make such ﬁle lists by hand. The task can read them as follows:
q_headlst2 @720.lista 720a.clk 720a.opt 720a.temp 720a.tel1 720a.tel2
q_headlst2 @720.listq 720q.clk 720q.opt 720q.temp 720q.tel1 720q.tel2,
This gives all the required header information that will be stored into the ﬁles *clk,opt,temp,tel1,
and tel2.

Table 4.2: Summary of a Standard Procedure: Spectroscopy
0. Data Examination
1. Making a Dark Image

Checking clock and optics
Dark per exposure

2. Read-out Pattern
Removing pattern noise
Reducing noise
3. Making a ﬂat for Spectroscopy Subtracting (Dark) from (Dorm)
Subtraction of detector pattern
Smoothing and normalization, if needed
4. The Target Images
Subtracting chopping data
Analysis of Standard Star(s)Dividing with Flat(a) Determination of spatial axis Making Sky Spectra of Target and Standard Star
1) (Object or Standard Star) − (Dark)
2) If multiple images, needs to average3) Dividing with Flat
Auto-Wavelength Calibration
Determination of spatial axis
(Tracing stellar position within standard star)
Axis determination
Making a conversion table for image transform(b) Image Transformation
Determination of Transformation Parameters
Transformation of target, sky and mask images(c) Calibrating Wavelength
Applying Auto-Wavelength Cal.
and Positions
on the image(s) from the above (b)Determination of disp. axis w. images from (b)(Tracing Stellar Spectrum
w. Standard Star)
Veriﬁcation of the Results(d) Shifting & Adding
Adding images after shifting along
the slit
5. Atmospheric Attenuation
Extraction of spectra
Correction w/Standard Stars

q headlst2
q list statq fcombine
q arithq subchq arithq subch
IRAF:gauss, q arith
COMQ images
q arithCOMA images
q arithq list statq arithq sky nlow (only for NL)
q startrace

q transtable2
IRAF: geomap
IRAF: geotranq sky nlow
q startrace

q arith, q fcombine
IRAF:imshiftq list stat

(Object)/(Standard Star)×(Standard StarTemplate)

6. Correction of Atmospheric
Attenuation using Standard Star

Determination of Conversion Factor to Flux
Determination of atmospheric spectrum
using data and template of standard star
Making 2D Images of Atmospheric Spectrum
Dividing Object Image with the 2D image

q list statq mkimgq arithChapter 5
Data Reduction: The Details
5.1

Imaging at N and Q Bands

5.1.1

Making Dark Frames

For many observations, dark data are taken either before or after scientiﬁc imaging. These
calibration data should be taken with the same exposure parameters (i.e. clock conﬁguration)
for the sources of interest, standard star or/and dome ﬂat. However, it is worth checking
whether or not they are identical because the current COMICS system allows the use of
diﬀerent frame numbers or/and CoAdd mode. Also, many users prefer to specify CoAdd
mode = 0 which gives RAW data or = 1 (ADD) without using chopping. In other words,
making a dark image gives dark counts per exposure for many COMICS observations.
Dark images can be made by (1) averaging images that have been taken with the same
clock conﬁguration, (2) obtaining mean images using the averaged ﬁles, and (3) dividing theimage by the exposure numbers. The following example indicates how to make a dark image
from images of #51219 and 51221.
Image No.
Clock Parameters
Resultant Dark (*)
51219,51221 PP=30, ND=1, Nexp=163, CoAdd=1, Ystart=1, chop=yes imgdark30 1(*):

imgdark30_1: Dark image for an exposure at PP=30, ND=1

The ﬁrst step is to make an averaged image over the observing time, namely, over theframes. This can be done with q list stat. To use q list stat, the user must properly
specify the following options: the ﬁle and detector names, pixel ranges (if you want to use allof them, use ”- -”), and ”:” which selects averaging along time axis.
The second step is to average the resulting ﬁles m51219 and m51221 using the task
q fcombine. Keep in mind that the original COMA00051219 and 51221 images have beencreated with CoAdd = 1 from the 163 exposure ﬁles. We therefore have to divide the image
by 163 to obtain a mean image per exposure (so that all the images have the same integrationtime). This can be done with the task q arith.
23q_list_stat ../DATA/COMA00051219.fits 1 - - : m51219
q_list_stat ../DATA/COMA00051221.fits 1 - - : m51221
q_fcombine m51219 m51221 ave=imgdark30
q_arith imgdark30 / 163 imgdark30_1
The user may have to repeat this procedure for all the types of observations to be reduced.
In most cases, these are for the source of interest(s) taken in the same night, standard star(s),dome ﬂat for spectroscopy (if exists) as well as all the data with identical parameter combinations of img/spc, PP, ND, and Ystart parameters. At this point, we suggest making alog sheet by hand for your records because you have to make many dark images with diﬀerent
clock conﬁgurations.

5.1.2

Making Flat Images

To make ﬂat images, the user may select either dome-ﬂat or sky-ﬂat. For imaging observations,the user must use ”self-sky-ﬂat” because dome-ﬂat cannot be used due to saturation. On theother hand, the user should use dome-ﬂat for spectroscopy observations as sky-ﬂat usually has
severe contamination from atmospheric emission lines whose wavelengths may have shifted.
In the procedure below, we suppose that the user has data available to make ”sky-ﬂat”.
To make ”self-sky-ﬂat”, you ﬁrst have to split the images for each chopping beam. This isbecause you will use ”oﬀ-beam sky” for the ﬂat of the ”on-beam sky” (or the vice versa) bymaking averaged sky images for each chopping beam.
Since COMICS focuses the incident light onto the slit, unevenness of the slit surface isbeing transferred into the ﬂat images in most cases. Therefore, you have to use sky images
taken with photometry mode for ﬂat ﬁelding, and those taken with the slit-width for eachslit-viewer. In addition, if you have used diﬀerent photometry ﬁlters, we suggest making sky
images for each ﬁlter using an oﬀ-beam sky image that have been taken without changing
ﬁlter positions. Again, this is because there are unevenness in gain of the detector pixels.
In principal, it should be possible to make ﬂat images taken with the same ﬁlter conﬁguration and apply it to all the data. However, in practice, this is not recommended because
diﬀerences in telescope pointings cause non-negligible diﬀerences in optics. Hence, we suggest
making ﬂat images using oﬀ-beam sky image for each target image.
It is well known that the result sky ﬂat has a residual gradient through the entire theimage. The amplitude of the gradient is usually ∼10% of the peak. If you use such a ﬂatimage, you will realize that the accuracy of the photometry is most likely worse than thatof the image made with the above ﬂat image. You can see the diﬀerence by comparing two
photometry results toward an object. In one case, the object is at the center of the ﬁeld ofview. In the other case, the object is seen close to to the edge of the ﬁeld of view. Therefore
it is required to remove such a ripple by ﬁtting the low-spatial frequency components for theuse of ﬂat ﬁelding.
0. The following two ﬁles are selected to analyze for imaging observations.

../DATA/COMA00050615.fits../DATA/COMA00050617.fits
1. Subtracting dark from the COMA images, then making a sky image
First, you have to know the numbers of Coadd (Nexp) of the data, which can be obtained
from the header by using Q CHEB. Since Nexp was 16 in this case, you must use a dark image
that has been obtained from an image created by adding 16 exposures. Then, simply multiply
a factor of 16 to the dark image already made for an exposure.
q_arith imgdark30_1 \* 16 imgdark30_16(\* indicates multiplying, not * in UNIX)q_arith ../DATA/COMA00050615.fits - imgdark30_16 sky50615.fits
q_arith ../DATA/COMA00050617.fits - imgdark30_16 sky50617.fits
The above operation has subtracted dark from all the frames. This is because q arith
works for all the Z with the conditions: (1) Z>1 for the images to be subtracted, and (2) Z=1
for the image to subtract. Here we obtained images of ”sky50615.ﬁts” and ”sky50617.ﬁts”
after subtraction of the dark.
2. Splitting images for each beam, then averaging each of them
The sky images at this stage contain the dual chopping beam data, which should be split
beam by beam. This can be done using the task of q bsep which creates images for individual
beam based on the information given in the header.
q_bsep sky50615.fits sky50615_p.fits sky50615_n.fits
q_bsep sky50617.fits sky50617_p.fits sky50617_n.fits
Here ”sky50615 p.ﬁts” corresponds to the data taken with chopping beam A, and ”sky50615 n.ﬁts”
with chopping beam B. You can check the results by referring to the header; the value of theparameter of NAXIS3 for ”sky50615 p.ﬁts” and ”sky50615 n.ﬁts” must be exactly half of thatfor ”sky50615.ﬁts”.
Subsequently we have to average the individual beam ﬁles with q list stat.
q_list_stat sky50615_p.fits 1 - - : sky50615_pa.fits
q_list_stat sky50615_n.fits 1 - - : sky50615_na.fits
q_list_stat sky50617_p.fits 1 - - : sky50617_pa.fits
q_list_stat sky50617_n.fits 1 - - : sky50617_na.fits
2. Imaging smoothing with a Gaussian beam
You may have realized that the ﬂat images that we obtained are most likely to contain alarge-scale residual pattern over the image, which should be eliminated for making a better
ﬂat. We therefore convolve a Gaussian ﬁlter to the above images, and divide them. For thispurpose, we use the task gauss in IRAF.

cl> gauss sky50615_na.fits[,,1,1] sky50615_naG.fits 6
cl> gauss sky50615_pa.fits[,,1,1] sky50615_paG.fits 6
cl> gauss sky50617_na.fits[,,1,1] sky50617_naG.fits 6
cl> gauss sky50617_pa.fits[,,1,1] sky50617_paG.fits 6This example shows the case where we have convolved with the Gaussian beam having
a FWHM of six pixels which is the best-searched value for cases, in our experience. The
resultant images are sky50615 nG.ﬁts.
Subsequently, we divide the images before Gaussian convolving by those after convolving,as follows.
q_arith sky50615_na.fits / sky50615_naG.fits sky50615_naF.fits
q_arith sky50615_pa.fits / sky50615_paG.fits sky50615_paF.fits
q_arith sky50617_na.fits / sky50617_naG.fits sky50617_naF.fits
q_arith sky50617_pa.fits / sky50617_paG.fits sky50617_paF.fits
The ﬁles of ”sky50615 naF.ﬁts” etc are the desired ﬂat images.

5.1.3

Making a Bad-Pixel Mask

Bad Pixels on the Flat
You may have realized that there are non-negligible numbers of pixels or/and region(s) where
the count values are obviously lower than those from the surrounding pixels. They are lesssensitive compared with those in the surrounding area, and most likely due to dusts on thedetector and/or unevenness of the gold coating of the slit. In addition, you also need to
examine the resultant ﬂat to see if there are pixels that have signiﬁcantly lower counts thanthose for the surroundings or/and pixels being aﬀected by artiﬁcial line(s). The tasks q badpix
and/or q 2val tell you such bad pixels or/and regions to be eliminated. The former task will
help you to identify bad pixels by making an image where 0 is stored for the pixels exceedinga threshold value you specify and 1 for the other pixels. For instance, if you want to identify
bad pixels with sensitivity less than 80% or more than 120% of the mean of the surrounding
pixels, issue the following commands:
q_badpix sky50615_naF.fits 0.8 1.2 sky50615_naF_bp.fits
q_badpix sky50615_paF.fits 0.8 1.2 sky50615_paF_bp.fits
q_badpix sky50617_naF.fits 0.8 1.2 sky50617_naF_bp.fits
q_badpix sky50617_paF.fits 0.8 1.2 sky50617_paF_bp.fits
Of course, 80% and 120% are for examples, and we advise using trial and error to ﬁnd thebest threshold for your data. Our experience suggests that this method should work for mostcases because the ﬂat images tend to have rather uniform values. However, if it does not work,
you may try to identify bad pixels by specifying threshold values with respect to the standard
deviation for the surrounding pixels rather than using the mean.

Pixels Having No Incident Photons due to Vignetting etc
The task q badpix may fail to identify bad pixels located at the edge of the image area thathave no incident photons. Unfortunately, since there is no general method that can work for
all type of data, you should consult with your local expert or/and the support astronomer,
depending upon the cases.
Making a Mask Image
With the procedure above, we assume that you have successfully obtained an image for bad
pixel masking; the image has values of 1 for the identiﬁed bad pixel while 0 for the others.
You can eliminate all the bad pixels from the target images by multiplying the bad pixel mask
image.

5.1.4

Reduction of Science Data

When you are happy with the ﬂat and bad-pixel masking images, the next step is to apply
these to the images containing the source(s) of interest.
Subtraction of Chopping Images — COMQ Images
You need to subtract the oﬀ-beam images from the on-beam ones for the background subtraction. In practice, however, the user can simply use the COMQ images instead of making
such images by hand. If you want to know image statistics for each or/and between frames
or/and exposures, you should use the task q bsep to extract the desired images before using
q list stat because COMQ images are the summations for each beam images.
Reduction of Read-out Pattern Noise
Since the COMICS system reads the detector current from all 16 channels at the same time,
the readout noise should have the same pattern noise for each of the 16 channels. Recall
that each channel has 20 columns. On the other hand, because the target source(s) is usuallydetected over some portion of the detector, you can reduce the readout noise using the pixels
where the source(s) is not detected.
For this purpose, use the task q subch which can subtract the median of the input image.
If your source(s) is much smaller than the image size, you can simply use all the channels withthe command of
q_subch ../DATA/COMQ00050615.fits obj50615.fits
On the other hand, if your source is extending over a number of images, the above subtraction with median may fail to remove the pattern. You can avoid such cases by specifying
channel number(s). For instance, if you want to use channel 1, namely X = 1–20 in the image,try the following,q_subch ../DATA/COMQ00050615.fits 1 obj50615.fits
For this case, however, you will see that the ”random” readout noise has increased, althoughthe “pattern noise” has been subtracted.
Dividing with the Flat Image
The next step is to divide the chopping-subtracted images, i.e., COMQ images, with the ﬂat.
Your COMQ images should contain both the positive and negative images as the chopping
has been done within the COMICS ﬁeld of view for many cases.
Start with ﬂat-ﬁelding for the positive image. First, divide the target images with the ﬂat
made from the beam-B (the negative) images. After that, remove bad pixels at this stage, as
follows,q_arith obj50615.fits / sky50615_naF.fits obj50615_datP0.fits
q_arith obj50615_datP0.fits \* sky50615_naF_bp.fits obj50615_datP1.fits
You can divide the ﬂat images for the negative ones as same as for the positive ones. But,be sure to multiply by −1 for the positive images. If your negative images do not appear in
the COMICS ﬁeld of view, you can skip this procedure.
q_arith obj50615.fits / sky50615_paF.fits obj50615_datN0.fits
q_arith obj50615_datN0.fits \* sky50615_naF_bp.fits tmp.fits
q_arith tmp.fits \* -1 obj50615_datN1.fits
We assume that you now have completed the calibration processes that are required foreach image and should have the ﬁnal positive images for each image ﬁle (the negative as well
if both of the beams are inside the ﬁeld of view).

5.1.5

Adding the Images

Adding the Images without Shifting
In principal, you should be able to add the calibrated images as long as they have been takenwith the same ﬁlter conﬁguration and their pointing centers agree with each other within themargins of the error. You can probably add the images that have been acquired sequentiallywithin a short time and/or your source(s) cannot be recognized in each images.
Your images can be added using the task q arith. In many cases, however, you must issue
the command several times, which is not always useful. Instead, you can use q fcombine,which is originally designed for averaging images. To use the task q fcombine, you have tosupply an input ASCII text ﬁle where the names of images are written in each line. The
following example shows how to utilize the task by supplying the objlist ﬁle that describes
names of the ﬁles to be added:

q_fcombine @objlist ave=objave.fits med=objmed.fits sig=objsig.fits
which gives you the averaged image (objave.fits), the median one (objmed.fits), and thestandard deviation one (objsig.fits).
If you want to keep the angular resolution (i.e., sharpness) of the original data, we advise
you not to add the image ﬁles that have been taken over a long observing time. This issimply because your source positions when measured with the detector coordinate have shifted
gradually during the observations. If your source(s) has suﬃcient S/N (and ideally is compact)
enough to measure its positions for each ﬁle, we suggest shifting the images before adding.
Shift-and-Add the Images
If your source(s) can be seen in diﬀerent pixel positions or/and if you wish to increase S/Nof the source by concatenating both the positive and negative images, you should consider
to shifting the images before adding them. In practice, you may succeed in obtaining thepositional oﬀset within an accuracy of ∼ 0.1 pixel by comparing the proﬁles with eye. However,there is a more precise method in the following.
1. Getting an oﬀset for image shifting
In order to add the fully calibrated images, you need to measure positional oﬀsets between
images. If you have a bright compact source(s) in your ﬁeld of view, you can use it for apositional reference to calculate oﬀsets using the task q photo. If you don’t have a compact
source but only an extended emission in your ﬁeld of view, you cannot take this approach.
Another method is to compare spatially ”averaged” proﬁles of the emission projected along
the X and Y axes for the images taken with diﬀerent ﬁlters or/and positions. These proﬁlesshould be calculated for a rectangular region which includes the target source in each image(see Figure 5.1.5). Comparing cross-correlations of the spatial proﬁles produced from theimages to be added, one should be able to determine the desired positional shifts for adding.
Figure 5.1.5 illustrates how this method works. Keep in mind that you must exclude columns
and/or rows that contain bad pixels before calculating the cross-correlation.
You have to repeat this process for all the images that you want to add. For the sample
images 50615 and 50617, we have obtained the shifts as follows,
obj50615_datP1.fits
00obj50615_datN1.fits -106.5 +0.8
obj5061_datP1.fits-0.2 +0.8
obj5061_datN1.fits +106.5 +0.8

Reference

2. Shifting images and obtaining a mean
Once obtained positional shifts for all the ﬁles, you can add them by supplying the shifts.
To shift images, we use task imshift.

cl> imshift obj50615_datP1.fits[,,1,1] obj50615_datP1s.fits 0 0
cl> imshift obj50615_datN1.fits[,,1,1] obj50615_datN1s.fits -106.5 0.8
cl> imshift obj50617_datP1.fits[,,1,1] obj50617_datP1s.fits -0.2 0.8
cl> imshift obj50617_datN1.fits[,,1,1] obj50617_datN1s.fits -106.5 0.8
It should be noted that the ﬁrst line above is to makes a template image for the shifted
images, although the given values are zero. The purpose for this is to deﬁne dimensions (i.e.,
pixel numbers) for the resulting image.
Subsequently, you can average the images. Again, if your images still have non-negligible
numbers of bad pixels, you should use the median averaging instead of the simple averaging.
Consider using the task q fcombine for both methods. In practice, it is wise to provide
an ASCII text ﬁle that describes names of the images to be concatenated in each line. For
example, suppose that your list ﬁle has name of ”objlist”, typeq_fcombine @objlist ave=objave.fits med=objmed.fits sig=objsig.fits
This should give you the averaged image (objave.ﬁts), the median one (objmed.ﬁts) andthe standard deviation one (objsig.ﬁts).

Each image(i.e., “exposure” or “frame”)

Dark image forseveral exposures

An averaged image
Desired dark image
between files
per exposure(i.e., “unnormalized” dark)

COMA files of
Dark images

An averaged image
along time axis

Figure 5.1: A typical procedure for making a dark image requires the user to average all thedark images taken with the same clock parameters, and to convert the averaged images intoa dark image per exposure.

Figure 5.2: Standard procedure for making a sky ﬂat image. In this case, the user can get thedesire ﬂat image(s) after subtracting the dark image, beam splitting, averaging, and removing
any large-scale residual pattern. Notice that a ﬂat image created from the one-beam data canbe used sky ﬂat for the data taken with the other beam.

Figure 5.3: Making a mask image with task q badpix: The task produces a mask image
where the bad or/and hot pixels whose sensitivity are high or low compared with those for thesurroundings. The user can specify threshold values to identify unusable pixels; the resultant
mask image has values of 0 for the bad pixels, and 1 for the usable ones.

position 1

position 2

objectI2(X,Y)objectI1(X,Y)

Compression along row

Region 2
X: 181-280
Y: 80-220

Region 1
X: 21-120
Y: 50-190

Count
F1(Y)

F2(Y)

Y

Compression along column

Correlation

Count
F1(X)

F2(X)

Count
F1(Y)

F2(Y-py)

X

Correlation

Count

Y

F2(X-px)
F1(X)X

Shift Values: (px, py)

Figure 5.4: A method to add images for the extended emission

Figure 5.5: A method of photometry data reduction. The ﬁnal images can be obtained bycorrecting for the ﬂat to the COMQ images, masking bad pixels, and adding all the imagesby shifting them (if necessary).

5.1.6

Miscellaneous

Correcting for Sky Levels
We suggest verifying the quality of sky-level subtraction carefully because, in our experience,
gradual level changes over the images tend to remain. This would be due to high backgroundemission, for example, and/or the presence of clouds even for a blank sky. If such a residual
pattern exists, you should consider applying another sky subtraction by making sky images
to remove it.
Conversion of the Image Intensity Unit: Getting the Flux Density Scale
After completing the ﬁnal image construction, you probably want to convert units of the image
from the counts to the brightness using standard star(s).
First, using the data base provided by Dr. Cohen and colleagues, you have to compute thebrightness of the standard star at the ﬁlter wavelength you observed. Select standard stars
which have been observed with an airmass close to that of your target. Since accurate measurements of ﬂux densities of the standard sources have not been generally well-establishedat MIR bands, except for the N-band, you will need to obtain ﬂuxes for the desired ﬁlters
by integrating spectral data provided in the Cohen et al. paper over the desired wavelength
range. For details of the ﬁlter properties, such as the center wavelength, and the passband
characteristics, please refer to the COMICS web page at
http://canadia.ir.isas.jaxa.jp/comics/open/guide/index.html where you can ﬁnd otheruseful links such as the data base by Dr. Cohen et al.
For simplicity in your calculations, you may integrate the spectrum over the ﬁlter bandpass adopting a rectangular shape, as the ﬁrst order of approximation. However, you should
calculate the ﬂux by considering all the parameters, e.g., the bandpass curve of the ﬁlter, that
depend on the wavelength.
Second, the user should know the count of the standard star(s). You can obtain counts
by calculating summation of the counts over a rectangular area containing the standard star,for example. This can be done with the task q photo, or photometry tasks in IRAF. Now
you should be able to obtain a conversion factor from the count unit to ﬂux density throughcomparison between the obtained counts and calculated ﬂux density. To convert source images
into ﬂux density units, you can simply multiply the calculated conversion factor to the fully
calibrated images. If your source is extended, you may want to convert the images intobrightness temperature (Tb ) scale. This conversion can be done by dividing the image with apixel ﬁeld of view. Taking photometry observations as an example, you can obtain the images
with brightness temperature scale by dividing 0.0166 arcsec2 pixel−1 .
Adding Images Taken with the Diﬀerent Filters
If you want to make a color image using multiple images taken with diﬀerent ﬁlters, you
should consider how to properly register these image centers across ﬁlters. One method isto refer to the point source(s) commonly seen in the multi-wavelength images. If this does
not help, you may have to employ the shift-and-add method described in §5.1.5. However,
you should be careful when using the latter method. Make sure of the validity of the image
registrations across the ﬁlters by considering the nature of the source(s) (i.e., emission mechanism). Moreover, we suggest making sure possible position shifts across the images taken
at diﬀerent wavelengths and/or diﬀerent observing times (as telescope tracking accuracy may
have changed during the observations).
Estimate of Image Noise Levels
Given that ﬂuctuation of the incident counts at MIR wavelengths is governed by the sky background emission, which is usually brighter than the target sources, the most straightforwardand most reliable methods to estimate image noise level is measuring the standard deviation
for the blank sky in the target images rather than measuring standard deviation of the target sources across the images. We believe that this approach would be valid for photometry
observations as well, where high accuracy is generally required. Of course, the above method
works only if your ﬁnal images do not have a residual global pattern of the sky backgroundAND the sky background level is considerably brighter than the source(s) of the interest.

5.2

Analyzing Spectroscopy (N-band Low Resolution;
NL) Data

Figure 5.6. shows overall ﬂow of the data reduction process for a point source taken in
low-resolution mode. The basic idea of the data reduction can also be applied to the otherspectroscopy mode observations.

5.2.1

Preparing for a Dark FrameFor reducing spectroscopic data, we have to use not only the science data but also the dome ﬂatimages. For this purpose, the user must prepare for dark images that have been obtained with
exposure parameters (e.g., clock conﬁgurations) identical to that for scientiﬁc observations.
This can be done with the same way described for the imaging mode (see § 5.1.1 for thedetails).
Assume that you have dark image ﬁles of COMA00051862.ﬁts (Nexp=98) for the ﬂat, andthose of COMA00051872.ﬁts (Nexp=32) for the target sources.
q_list_stat ../DATA/COMA00051862.fits 1 - - : m51862
q_arith m51862 / 98 spcdark50_1
dark images for the flat
q_list_stat ../DATA/COMA00051872.fits 1 - - : m51872
q_arith m51872 / 32 spcdark150_1
dark images for the target source(s)It should be again noted that we have used the dark data taken with the same exposure
parameters. For many cases, the users can use dark images that were taken for the ﬂux
measurements. If you don’t have such a data set or/and you have only a dark image obtained
from all the detectors, you probably have to analyze data detector by detector (see §2.1 aswell).
To see the contents — detector information — of each FITS ﬁle, check the Q DETST parameter in the header. This key word is given by a number consisting of 0 and/or 1. The ﬁrst
digit indicates a ﬂag of the imaging or photometry detectors, and the subsequent ones give
status of the spectroscopy detectors. For example, 100110 indicates that you have obtained
a photometry image plus spectroscopy taken with the detectors No. 3 and 4. Notice that thedesignation number of the spectroscopy detectors and the 4th axis, i.e., the w-axis, in the FITS
ﬁles diﬀer from each other. In this example, w=1 corresponds to the spectroscopy detector
No.3 and w=2 does No.4.

5.2.2

Making a Flat Frame

For spectroscopy observations, use dome ﬂat images for making a ﬂat. We suppose thatyou have obtained more than one ﬂat image at each grating position. We strongly suggest
using the dome ﬂat image(s) that have been taken immediately before or/and after the source
integration. If you don’t have such calibration exposure, use those taken as close as possible

Spectral Imageof標準星分光画像
Standard Star

Spectral Imageof 天体分光画像
your target

Flat Image ofDome (Cell) Coverドーム(セルカバー)フラット画像
画像選別
ノイズ低減
Selecting
Imagesflatで割る
Noise Reduction
Flat Fielding

Image
画像変換
(直交化・波長較正)incl.,
Transformation,orthogonization,
wavelength cal.

Darkをひく
Subtructing
Dark
規格化
Normalization

スカイの処理Sky
subtraction

Extractionof spectrum
切り出し(１次元化)

Spectra in 1D

１次元化スペクトル

Atomosperic Corrections

大気吸収補正i.e., dividing by spectrum of(標準星Sptで割り算・
standard star,
templateのかけ算)

multiplyig template

Correction
オゾン補正of Ozone
フラックス・傾き補正
absorption

最終的なスペクトルthefinal spectrum

Figure 5.6: A method to reduce COMICS spectroscopy data. An example for a point-likesource with the low-resolution modeto the scientiﬁc exposure because ﬂat images usually change even with a slight shift of grating
conﬁgurations. In this example, ﬁrst average a COMA image of COMA00051356.ﬁts andsubtract the dark. Since the ﬂat image has Nexp =9, we have to multiply by 9 before thesubtraction as shown in the below.
q_list_stat ../DATA/COMA00051356.fits 1 - - : flat51356.fits
q_arith spcdark50_1 \* 9 spcdark50_9
q_arith flat51356.fits - spcdark50_9 flat51356_d.fits
You will realize that the resulting dome ﬂat contains an un-removed gradient over the ﬁeldof view, which is similar to those seen in the ﬂat for the photometry. Therefore, you need todivide the image with that is convolved with a Gaussian ﬁlter with the same method as we
did in §5.1.2.
cl> gauss flat51356_d.fits[,,1,1] flat51356_dG.fits 20
This will give you an image convolved with a Gaussian with FWHM of 20 pixels. The
width of 20 pixels is recommended, based on our experience.
q_arith flat51356_d.fits / flat51356_dG.fits flat51356_dF.fits
Now we assume that you have successfully obtained a ﬂat image of at51356 dF.ﬁts.

5.2.3

Reduction of the Read-out Pattern Noise

For some observing modes that do not use all the detectors, such as medium-resolution spectroscopy at the N band, users are allowed to reduce the readout pattern noise by subtracting
the readout from the reference detector. This is because one can use read-out from the unused
detector. This will allow you to reduce the readout pattern noise by subtracting the readoutfrom the reference detector.
If you wish to attain high accuracy subtraction, you may try to subtract a median patternimage calculated from the reference detector channels that did not contain signiﬁcant noise.
This is because each detector consists of 16 channels with x=20 pixels and all the channels
are being read at the same time.
To obtain the median for all the channels, you can use the task q subch. For this purpose,
you should ﬁrstly extract a detector image from the reference image of read-out pattern noise.
For instance, if you want to analyze a source spectral image of COMQ00051308.ﬁts (= standardstar), use
q_list_stat ../DATA/COMQ00051308.fits 2 - - 1 chnoize.fits
The next is to calculate median image over the 16 channels from the reference image ofthe read-out pattern noise, usingq_subch chnoize.fits REF chnoize_REF.fits
After that, subtract the noise pattern image of chnoize REF.ﬁts from the spectroscopy
data as follows
q_list_stat ../DATA/COMQ00051308.fits 1 - - 1 obj51308.fits
q_arith obj51308.fits - chnoize_REF.fits obj51308s.fits
which will probably help to reduce the noise levels signiﬁcantly — hopefully down to the
30 – 40% level.

5.2.4

Data Reduction of Scientiﬁc Images

You have completed the noise reduction processes across channels to proceed ﬂat ﬁelding. Now
you are ready to reduce the scientiﬁc images you have completed to make ﬂat image(s). You
can subsequently apply these ﬂat images to the source images.
Adding All the Images
If part of your data were acquired at the same position in the sky, you may consider adding
all the appropriate images, which makes in much simpler for further data reduction. If thisis not the case, we suggest keeping them without adding at this stage. For example, assume
that the two spectral images ﬁles of COMQ00051314.ﬁts and COMQ00051316.ﬁts have beentaken at the same positions in the plane of the sky and that you have obtained obj51314s.ﬁts
are obj51316s.ﬁts after applying the read-out pattern noise reduction described in $5.2.3. The
following command will collapse the two images into an image.
q_fcombine obj51314s.fits obj51316s.fits

ave=objspc.fits

Dividing the Flat
The next step is to divide the added images (i.e., objspc.fits) (or the read-out pattern noisereduced single image, e.g., obj51308s.ﬁts in $5.2.3, depending on your reduction process) by
the dome ﬂat of at51356 dF.ﬁts. Taking the latter case as an example, use
q_arith obj51308s.fits / flat51356_dF.fits obj51308s_f.fits
Orthogonization of the Spectroscopic Images and Wavelength Calibration
After completing reduction of the read-out pattern noise, you have to transform the images toorthogonize dispersion and spatial axes, which will be used to extract spectra of the emission
you desire. Furthermore, you normally calibrate the absolute scale of wavelength at thesame time by ﬁnding pixel-wavelength correspondences. This subsection describes how suchdispersion and spatial axes are determined (see Figure 2.4).

1. Determination of Dispersion Axis
The dispersion axis is usually determined through comparison to an atmospheric sky emission image. The correction equation for the wavelength can be written as λ[µm] = Ax[pix]+B
for y. Since the coeﬃcients of A and B vary with y as shown in Figure 2.4, ﬁnd their ydependency as A = a0 y + a1 and B = a2 y + a3 . Therefore, a set of the parameters of a0 , a1 ,
a2 , and a3 represent the desired correction equation.
As the ﬁrst order of approximation, a wavelength calibration equation can be expressedby λ[µm] = Ax[pix] + B for y. However, since the coeﬃcient of A and B varies with y formany cases, you have to determine their dependencies on y as A = a0 y + a1 and B = a2 y + a3 ,respectively.
For medium resolution spectroscopy, we suggest splitting data ﬁles for each detector to
make further reduction processes simpler. In practice, ﬁrst make a sky image that can beobtained by averaging images of standard star(s) or/and COMA images of scientiﬁc images(you have to subtract the appropriate dark image and divide it with the ﬂat one). You can
now make a sky image from spectral image of the standard star (COMA00051308.ﬁts).
q_list_stat ../DATA/COMA00051308.fits - - - : skyimage.fits
q_arith spcdark150_1 \* 3 spcdark150_3(Making a dark image
adjusting with Nexp)
q_arith skyimage.fits - spcdark150_3 skyimage_d.fits(Subtracting the dark)
q_list_stat skyimage_d.fits 1 - - - skyimage_d1.fits
q_arith skyimage_d1.fits / flat51356_dF.fits skyimage_d1f.fits (Dividing with flat)
This dispersion axis can be applied to the obtained parameters to the other data sets, unlessthe grating was shifted during observations, because such a calibration curve does not change
for many cases. However, users are advised to double-check wavelength calibration parameters
obtained from each object or/and those form the mean images (if you have sequentially taken
the same conﬁguration images.). On the other hand, if you have grating(s) that may haveshifted during observations, we suggest determining calibration parameters for each grating
position(s) because repeatability of grating parameters is not completely guaranteed.
As for the low-resolution observations at the N band, the dispersion axis can be obtainedusing the task q sky nlow.
q_sky_nlow skyimage_d1f.fits 1 - default 1 2 > skyimage_d1f.res
This command calculates cross-correlation with the database of the atmospheric emission
lines at 20 points in the y axis between y =30 and 220, then guessing a wavelength and pixel
correspondence λ[µm] = Ax[pix] + B. The identical algorithm may be applied for the otherwavelength observations, however, we have not developed such a command for them.
The resulting ﬁtting parameters are stored in skyimage d1f.res and can be checked withstandard graphic packages such as GNUPLOT. In skyimage d1f.res, the third column indicates the extracted y-axis value, and sixth and eighth column the best-ﬁt values of A and B,

respectively. The ﬁtting result at y = 30 may have failed because of the photons at y = 30 are
more likely to be vignetted. If you ﬁnd such result, eliminate them manually from the ﬁtting
data.
Subsequently, apply a linear function to the ﬁtted parameters of A and B by least-square
ﬁtting, namely, use ﬁtting equations of A = a0 y + a1 and B = a2 y + a3 (see Figure 5.7). In
conclusion, use a parameter set of a0 , a1 , a2 and a3 to repeat the best-estimated wavelength
constant axis.
0.02

7.73
"20140skyfit" u 3:6
2.142891e-07*x+1.987821e-02"20140skyfit" u 3:8-4.186469e-04*x+7.737431e+00

0.01998

7.72

0.019967.71

0.019947.7
0.01992
7.69
0.01997.68
0.019887.67
0.01986
7.66

0.019847.65

0.019820.0198

7.64
20

40

60

80

100

120

140

160

180

200

220

20

40

60

80

100

120

140

160

180

200

220

Figure 5.7: A result of wavelength and pixel relationship determination and a plot of linear
function ﬁtting

2. Determination of Spatial Axis
To determine the spatial axis, use images of standard star(s) (or any point-source image
should work as well) that have had detector response subtracted and divided with a ﬂat.
You have to ﬁnd the peak position of the star in the x-axis, then ﬁt it with equation of
y = b0 (x − b1 )2 + b2 to ﬁnd b0 , b1 , and b2 . Since the parameter set of b0 , b1 and b2 is thought
to be stable during an observation, you can use a point source image that has a high S/Nor/and an image that has been obtained by averaging several images without any positional
shifts. In this manual, we use the reduced standard star image of obj51308s f.ﬁts using thetask q startrace.
Assume that the image contains the standard star at the positions of y = 62 –92, use
q_startrace obj51308s_f.fits 1 30-290 62:92 1 > obj51308s_f.str
this searches peak positions and stores the results in obj51308s f.str. If the task fails
to ﬁnd peak positions (i.e., failed to converge) or it takes a quite long time to converge,some tips in the FAQ may help you. The second column of obj51308s f.str stores the xcoordinate and the tenth y of the peak, the user should plot them to verify the results (please
use GNUPLOT etc); see the example shown in Figure 5.8. In some wavelength ranges such as

the O3 bands where very intense emission lines are normally seen, the task may fail to makea ﬁt. If this is the case, we suggest excluding such points before ﬁtting. For the medium
resolution observations, the procedure above should be repeated for all the detectors.
The ﬁtting can be done with a quadratic expression of y = b0 (x − b1 )2 + b2 as shown in
Figure 5.8. The selection of the quadratic expression is due to the structure of the COMICS
optics. In the case of medium resolution spectroscopy, one can approximate distortion as alinear equation because the spectrum spread over the ﬁve detectors. If you try to ﬁt with the
2nd-order equation, you would not be able to get a reasonable result for the detectors at theedge, where intensity of the incident wave usually decreases. We therefore suggest ﬁtting witha linear equation for the medium resolution mode. The resultant coeﬃcients of b0 , b1 , and b2
are the desired parameters for describing the spatial axis.
142"STD.pos2.str" u 2:10
1.4970369e-05*(x-258.51329)**2+141.07376
1.4705852e-05*(x-264.2791)**2+141.06172
141.8

141.6

141.4

141.2

141

140.8
0

50

100

150

200

250

300

Figure 5.8: An example of a space-constant line and the ﬁtting results

3. Image Transformation
Using the above results, you can transform the images with IRAF. The task geomap ﬁnds
the transformation parameters and the task geotran transforms your image. First, produce
an image that will be transformed with geotran as follows,
q_transtable2

a0 a1 a2 a3 b0 b1 trans.dat aa bb

where a0, a1, a2, a3, b0, and b1 are the parameters obtained in the above, trans.dat isa name for the output image, and aa and bb deﬁne a wavelength-pixel relationship (λ[µm] =aa × x[pix] + bb) after the transformation. In principal, aa and bb can be arbitrarily deﬁned.
However, you should use the same parameters obtained from the λ-pixel relation of the target.
This is because one can minimize loss of the information after the transformation by reducing
the numbers of pixels that are transformed outside the ﬁeld of view. In the above example,
we obtained the following parameters for the medium resolution at the N band.

a0=2.669213e-07
a1=1.990141e-02
a2=-3.758642e-04
a3=7.578833e+00
b0=1.5874755e-05
b1=195.9228
Now we can issue the IRAF geomap command for the directory where login.cl ﬁle exists,clwhich starts IRAF, and you should move to the directory for the analysis. Using epar command, set input parameters for geomap as shown below.
I R A F
Image Reduction and Analysis Facility
PACKAGE = immatchTASK = geomap
input=
database=
xmin
=
xmax
=
ymin
=
ymax
=
(transfo=(results=(fitgeom=(functio=(xxorder=(xyorder=(xxterms=(yxorder=(yyorder=(yxterms=(reject =(calctyp=(verbose=(interac=(graphic=(cursor =(mode=

The input coordinate files
transpar The output database file
1. Minimum x reference coordinate value
320. Maximum x reference coordinate value
1. Minimum y reference coordinate value
240. Maximum y reference coordinate value) The output transform records names) The optional results summary files
general) Fitting geometry
polynomial) Surface type3) Order of x fit in x3) Order of x fit in yhalf) X fit cross terms type3) Order of y fit in x3) Order of y fit in yhalf) Y fit cross terms type
INDEF) Rejection limit in sigma units
double) Computation typeyes) Print messages about progress of task?
no) Fit transformation interactively ?
stdgraph) Default graphics device) Graphics cursor
ql)You should also set the parameters for geotran as follows
I R A F
Image Reduction and Analysis Facility
PACKAGE = immatchTASK = geotran
input=
output =database=transfor=(geometr=
(xin
=
(yin
=
(xshift =(yshift =
(xout
=
(yout
=
(xmag
=
(ymag
=
(xrotati=(yrotati=
(xmin
=
(xmax
=
(ymin
=
(ymax
=
(xscale =(yscale =(ncols =(nlines =(xsample=(ysample=(interpo=(boundar=(constan=(fluxcon=(nxblock=(nyblock=(verbose=(mode=

Input data
Output data
Name of GEOMAP database fileNames of coordinate transforms in database filegeometric) Transformation type (linear,geometric)
INDEF) X origin of input frame in pixels
INDEF) Y origin of input frame in pixels
INDEF) X origin shift in pixels
INDEF) Y origin shift in pixels
INDEF) X origin of output frame in reference units
INDEF) Y origin of output frame in reference units
INDEF) X scale of input picture in pixels per reference
INDEF) Y scale of input picture in pixels per reference
INDEF) X axis rotation in degrees
INDEF) Y axis rotation in degrees
INDEF) Minimum reference x value of output picture
INDEF) Maximum reference x value of output picture
INDEF) Minimum reference y value of output picture
INDEF) Maximum reference y value of output picture1.) X scale of output picture in reference units per1.) Y scale of output picture in reference units per
INDEF) Number of columns in the output picture
INDEF) Number of lines in the output picture1.) Coordinate surface sampling interval in x
1.) Coordinate surface sampling interval in ylinear) Interpolant
nearest) Boundary extension (nearest,constant,reflect,wra0.) Constant boundary extensionyes) Preserve image flux?
512) X dimension of working block size in pixels
512) Y dimension of working block size in pixels
yes) Print messages about the progress of the taskql)

geomap can be executed as follows
im> geomap trans.dat transpar 1 320 1 240

The command above has prepared for the image transformation by setting the parameters.
You need the following: image ﬁles of the targets, standard star(s) that must be subtracted
the detector response and be divided by the ﬂat, sky image (skyimg) that has been used forthe wavelength calibration, and the masking image. The sky image will be used to calibrate
wavelength after image transformation, and the masking image to eliminate the eﬀect of bad
pixels in the transformed image. Since IRAF accepts input ﬁles with an extension of ”.ﬁts”
only, you have to give it explicitly. In the example, we issued the command as follows,
im> geotran obj51308s_f.fits[1:320,1:240,1,1] obj51308s_ft.fits transpar trans.dat
Here, the ﬁrst and second input parameters can be also written with the form of @ﬁlelist.
4. Veriﬁcation of the results and determination of the wavelength calibration
equation
After image transformation, one should verify the results (i.e., the spatial axis should be
horizontal and dispersion one vertical) using the tasks q startrace and q sky nlow. Noticethat the output ﬁle of geotran is a 3D cube, whereas q startrace and q sky nlow deal withthe 4D image only. We therefore convert the geotran output ﬁles into a 4D image.
Next, check the results on the spatial axis. For this purpose, we should multiply q startrace
to the transformed image (obj51308s ft.ﬁts) and plot the results (using e.g., GNUPLOT) as
shown in Figure 5.8 where the peak position of the signal should be constant regardless of thex-value. In practice, one can do such a procedure using the following commands:
q_chgaxis 4 obj51308s_ft.fits obj51308s_ft4.fits
q_startrace obj51308s_ft4.fits 1 30-290 62:92 1 > obj51308s_ft4.res(Verify the plot of obj51308s_ft4.res)
The ﬁrst line converts obj51308s ft.ﬁts to a 4D image, and the second one is the command
to multiply on q startrace.
The next step is to check the results of the dispersion axis, which can be done by transforming the sky image (skyimage d1f.ﬁts) and applying q sky nlow to the results. Make surethe values of A and B are constant with respect to the y-axis.
im> geotran skyimage_d1f.fits[1:320,1:240,1,1] skyimage_d1ft.fits transpar trans.dat
q_chgaxis 4 skyimage_d1ft.fits skyimage_d1ft4.fits
q_sky_nlow skyimage_d1ft4.fits 1 - default 1 2 > skyimage_d1ft4.res(Verify the plots of skyimage_d1ft4.res)
Now we can determine the ﬁnal version of the coeﬃcients for the wavelength-calibration
equation. (i.e., the input parameters of aa and bb to the q transtable2 command) using theresults from skyimage d1ft4.res.

Further subtraction of the residual pattern
It is best if all the sky level variations could be removed using chopping. However, in practice,
chopping frequency tends to be slower than sky level variations. If this is the case, you should
see some residual pattern in the chopping-subtracted images. Since you have completed theprocess of orthogonizing the images (i.e., x-axis: wavelength, and y-axis: space), you shouldbe able to estimate such a low-level residual sky pattern for the emission-free region(s) toapply a further correction. For this purpose, we should remove the median value with respect
to the y-axis (i.e., dispersion axis), which can be done with q submedrow, as follows:
q_submedrow obj51308s_ft.fits 1 320 150 240 obj51308s_ftc.fits
which gives the ﬁnal chopping image corrected for the residual sky pattern noise.

5.2.5

Extracting the Spectra

Now you can extract the desired spectra from the obtained 2D data. You may want to use tryand-error to ﬁnd an appropriate aperture width (e.g., 3, 5, 7, and 9 pixels) along the spatial axis
to maximize the S/N of the spectrum. Here, you can measure the noise level through standard
deviation of the emission-free region (i.e., the blank sky region). For instance, suppose thatyou have an image which accommodates the source at y = 72. In this case, you can estimate
the noise level for the blank region of e.g., y=100 – 130.
q_list_stat obj51308s_ftc.fits 1 15-295 110:130 1 >! obj51308s_ftc.noise1pix
In this example, the sixth column of obj51308s ftc.noise1pix indicates the noise level foreach x-value which corresponds to the wavelength. Next, extract spectra of the source at any
y position. Assuming that you want to have a width of three pixels (i.e., y = 71 – 73), thefollowing gives the desired spectrum.
q_list_stat obj51308s_ftc.fits 1 15-295 71:73 1 >! obj51308s_ftc.signal1pix
The 5th column of obj51308s ftc.signal1pix is count value at the x-position (corresponding
to the wavelength). Notice that this procedure calculates averaged value over the speciﬁed
pixel
range. With this reasoning, multiply a factor of 3 to obtain the added count value, and√3 to estimate the corresponding noise level.
awk ’{print $2,$5*3}’ obj51308s_ftc.signal1pix >! obj51308s_ftc.signalAdd
awk ’{print $2,$6*sqrt(3)}’ obj51308s_ftc.noise1pix >! obj51308s_ftc.noiseAdd
It would be good idea to store these results into a ﬁle using the paste command on UNIXlike,
paste obj51308s_ftc.signalAdd obj51308s_ftc.noiseAdd | awk ’{print $1,$2,$4}’ >!
obj51308s_ftc.SN

You now have a ﬁle of objspc Ftc.SN that stores the signal [ADU] and noise [ADU] at thecorresponding x[pix].

5.2.6

Adding the Spectra

As explained above, one can simply add any 1D spectra which have been created with thesame extracting parameters in the case of point-like sources. For diﬀuse emission, one has
to perform image registration prior to concatenating imges. The spatial registration between
images is normally done with the other object(s) seen in the slit viewer or measuring peak
position on the spectroscopy side. Recall that a pixel corresponds to 0.13′′ for the slit viewer,i.e., the imaging side, and 0.165′′ for the spectroscopy side. Because of this diﬀerence, thismethod is valid only for objects that have shifted along the slit. If your object has shiftedalong the perpendicular direction, you have to deal with this as another slit position data.
In summary, the shifting and adding procedure for the spectroscopy data is essentially thesame as for that in the photometry observations; the only diﬀerence is not to shift along the
wavelength direction, i.e., the x axis.

5.2.7

Wavelength Corrections

The spectral and noise data you have obtained are represented by count with respect to pixel
numbers, which should be converted into wavelength units. Since the relationship between
wavelength and pixel are already known from the transformed sky image (recall that the bestestimated relationship is written in the skyimage d1ft.res ﬁle), use those parameters. This
gives a = 1.9900 × 10−2 and b = 7.5774, then the wavelength correction equation can be
written as,
λ[µm] = 1.9900 × 10−2 × x + 7.5774
If you are familiar with awk, typecat obj51308s_ftc.SN | awk ’{print 1.9900e-2*$1+7.5774,$2,$3}’ > obj51308s_ftc.SN.um

5.2.8

Correction of Atmospheric Absorption and Eﬃciency

Dividing with the Spectrum of a Standard StarAs shown in Figure 5.9, both the spectra of standard star(s) and the object(s) have beenequally attenuated by the atmospheric transmission (Tλ ) and aﬀected by the instrumental
eﬃciency. Assuming that Tλ is identical both for the standard star and the object,
Fobsobj = Ftrueobj × Tλ ,Fobsstd = Ftruestd × Tλ

Assuming that Ftruestd is known, we can estimate Ftrueobj from Fobsstd
Ftrueobj = Fobsobj /Fobsstd × Ftruestd
In conclusion, you can obtain an ”atmosphere-free” spectrum by dividing the source spectrumwith that of the standard star and multiplying the ”atmosphere-free true” spectrum (alsocalled ”true” spectrum) of the standard star (See the spectrum of (C)/(A) × (B) in Figure
5.9).

STD

STD
(B)

λ

(A)

atmosphereslit efficiency
other efficiencyµµ

objectobject(D)

λ(C)

atmosphereslit efficiency
other efficiency

??
µ

µ

Figure 5.9: A method to obtain spectrum of the object using the standard star. The spectrum obtained by (A)/(B) is a spectrum of the atmosphere including the instrumental response. If you observed a standard star and the object(s) with the identical slit eﬃciency,(D)=(C)/[(A)/(B)] should give the desired spectrum. However, because the slit eﬃciencies
for a standard star and an object diﬀer from each other in many cases, the user must multiply
some scaling factors, see the text for the details.

”True” spectra of the standard stars (Ftruestd ) can be found in Cohen et al. (1999) as
templates. They can be obtained fromhttp://canadia.ir.isas.jaxa.jp/comics/open/guide/index.html.
Correction of Slit Eﬃciency
Since COMICS adopts the slit spectroscopy method, all the photons do not fall onto the slit.
Therefore, you have to consider slit eﬃciency. In general, slit eﬃciency depends not only on
the incident angle but also wavelength dependency because the size of stellar images vary withwavelength. Recall that slit width is comparable to diﬀraction limit. Moreover, one cannot
estimate slit eﬃciency if the stellar image is ”dancing”. With these reasonings, we have tocorrect slit eﬃciency not only for diﬀuse emission but also for point-like sources. Keep in mind
that slit eﬃciencies of the source and standard star are not identical even for the former cases.

In order to correct slit eﬃciency properly, you need a multi-color data set. Assume that(1) we have two color images, (2) the spectrum obtained so far can be written as Fobsobj (λ),
and (3) the eﬃciency equation is deﬁned by ε(λ) = a × λ + b when a and b are constant. In
this case, you should ﬁnd the coeﬃciencies of a and b so that Fobsobj (λ), the spectrum obtained
from the analysis so far keeps a relation of′
Fobsobj(λ) = 1/ε(λ) × Fobsobj (λ)

with those from the photometry observations. Finally, it is worth pointing out that slit
eﬃciency correction works well with a linear correction in many of the cases.

5.2.9

Diﬀuse Emission

In principal, you can analyze spectral line images of the diﬀuse emission with the same manneras done for the point-like sources. Namely, you can extract the desired spectrum at arbitrary
positions. However, if the emission is widely extended across the images, it would be ineﬃcientand painful to extract spectra at every position for correcting all of them one by one. In thiscase, it would be wise to apply the correction directly to the 2D image before extracting
the spectra. For this purpose, we ﬁrst have to calculate a correction function (which can be
described as
F (λ) ≡

Ftruestd (λ)
Fobsstd (λ) × ε

where Fobsstd denotes the resulting spectrum of the standard star, Ftruestd its ”true” spectrum, and ε slit eﬃciency.
Since the correction factors are functions of the wavelength, F (λ) can be applied to thedesired row by inverting the equation as is done for the wavelength. This allows you to obtaina relationship between pixels and wavelength to make 2D images that will be multiplied withthe reduced scientiﬁc image (objspc Ftc.ﬁts). The procedure to make a 2D correction function
can be done with task q mking.
For instance, if your correction function is stored in a ﬁle called calib.dat where the ﬁrst
and second column are X=1 and 2, respectively, such a correction should be done as follows,echo 1 | awk ’{for(i=0;i<240;i++){print"cat calib.dat"}}’ | sh > tmp
q_mkimg tmp calib.fits 320 240
q_arith objspc_Ftc.fits \* calib.fits objspc_Ftc_calib.fits

5.3

Frequently Asked Questions

5.3.1

Contamination of the Ozone Line at 9.6 µm

You will realize that it is not easy to divide source spectra by that of standard stars around
9.3 – 10 µm caused by the absorption lines of due to the ozone in the upper atmosphere. On
the other hand, you may not encounter such a problem if both the spectra have been takenwith very close airmasses. Our experience suggests it is likely that airmass mismatching would
be responsible for a wavy artifact, or/and bump/dip around 9.3 – 10 µm.
Such an eﬀect may be recovered partially (but not completely) with the following fashion. Using ATRAN (see Lord, S. D. 1992, NASA Technical Memorandom No. 103957), we cancalculate the expected atmospheric transparency at the wavelength with the observed zenith
angle, then apply it to the spectra. Of course, we have to wisely select the input parametersof ATRAN, such as the representative parameters of the atmospheric conditions such as PWV
and number of layers. It may be not trivial to ﬁnd the appropriate numbers of these parameters, but you can guess them from the weather logs (including from other observatories suchas nearby sub-mm telescopes). If you cannot, the standard values would be a good starting
point, and you will become known as an expert 10 µm observer from the ground!

5.3.2

Should I perform nodding for my observations?

For a standard MIR imaging, we usually employ not only secondary mirror chopping but
also nodding. This observing method will produce a set of data containing four images. The
reason to employ chopping and nodding is to reduce the residual pattern in the image after
subtraction of the two images taken with the diﬀerent optics. Such a residual pattern can beremoved by subtracting the nodding images. However, since the COMICS system on Subaru
telescope has a rather low background noise level originated from system the system, theresidual pattern is usually not so obvious. Therefore, chopping should work pretty well
ONLY if the expected intensity of the target is strong enough compared with theresidual pattern. In practice, our experience also suggests that a typical residual patternat the N 11.7 band is 35 mJy arcsec−2 . If your source is brighter than this value, we suggest
employing only chopping mode. Of course, if your target is expected to be fainter than thelevel, you will deﬁnitely want to employ nodding.

5.3.3

I have a ”ghost” in my low-resolution image at the N band,
what can I do?

Unfortunately, the current system often causes such artifacts. As a result, it’s almost impossible to use such images for scientiﬁc analysis. The only way to reduce such ”ghosting” is to
conﬁgure your imaging ﬁlter at through the ”hole” position (but please be aware that thisdoes not guarantee perfect removal of the ghost) for future spectroscopic observations. Whenyou observe with through hole, please keep in mind that you have to decrease the gain so as
not to saturate the slit viewer.

5.3.4

I get errors when I transform my image with geotran

For many cases, the error is caused by mismatching the image dimensions or/and the image
dimension is not being recognized properly. If this is the case, we suggest explicitly giving thedimensions of the input ﬁle, e.g., [1:320,1:240,1].

5.3.5

geomap relevant questions

If you are working on geomap several times on the same ﬁle, it may cause problems because
the task is trying to write the same key words to the created database ﬁle. We suggest deletingthe output ﬁle (or rename it) before repeating the geomap procedure.

5.3.6

It takes a looooooooooooooooong time for task q startrace
to complete ...

This is likely because the initial guess given in the program is far from the ﬂux of your target.
Does your source have a rather large signal value? If this is the case, try to divide the image
by certain number, e.g., 1000. On the other hand, if the sky background level is high and S/Nof the source is low, the program may fail to converge around the x-value. The key is to avoid
such x value(s) and try to run q startrace again after removal of such points.