TriCCS 画像解析(分光版)
http://www.kusastro.kyoto-u.ac.jp/~iwamuro/TRICCS/_index.html

岩室 史英 (京大宇物)


●注意

    この解析手順は i-band での非常に明るい装置内迷光があることを想定した手順です(多分 2023年中
    には解決されています)。装置内迷光が解消された以降のデータに対して適用したい場合は、
    iwamuro@kusastro.kyoto-u.ac.jp までご連絡下さい。連絡頂いてから修正版を準備します。

●ソース(随時更新)

  • headerinfo.c
  • sigclip.c
  • caldx.c
  • caldy.c
  • caldy.plt
  • caldxy.sh
  • cal_g.dat
  • cal_r.dat
  • cal_i.dat
  • fixpix.c
  • fitswarp.c
  • mksflat.c
  • shift_ave.c
  • spcomb.c
  • trcalib.c
  • vecadd.c
  • onespec.c
  • ap2.dat
  • ap4.dat
  • ap6.dat
  • 上記全部

    改変履歴:
    2021.09.01: とりあえず完成
    2021.09.08: バグ出しと改良でほぼ全て更新
    2021.09.12: trcalib.c でスペクトルの傾きに沿って aperture を適用するよう修正
    2021.09.12: sigclip.c で入力が 1 frame でも出力を出すよう修正
    2022.09.13: headerinfo.c で1枚目の平均値と 15000ADU 超え pix 数を表示するよう修正
    2022.09.13: spcomb.c で vflat の縞パターン除去の手順を修正
    2022.11.19: caldy.c で輝線検出の感度を上げ、画像端で輝線判定が出た場合のバグを修正
    2022.11.19: caldxy.sh 後の波長エラー確認用 gnuplot script として caldy.plt 追加
    2024.03.14: sigclip.c で中央 200x200 領域の x 位置を光量重心に変更
    2024.03.14: mksflat.c で配列定義のバグ修正、NaN 領域の扱いを修正
    2024.03.14: caldx.c で、縦縞の検出範囲の横を広げ縦を狭める修正

実際の解析手順は一番下にあります

●ポリシー

  • dark と flat の連続露出の重ね合わせは sigclip で行う
  • flat は、空間・波長特性を除いて用いる
  • Badpix は dark の安定性と flat のカウントで自動判別する
  • 連続露出の天体画像は pixel 単位で位置合わせして重ね合わせ、CR 除去も行う
  • 1 pixel 未満の形状補正は1度だけ行なう

●下準備

  • header 情報確認
    ./headerinfo ディレクトリ [KeyWord1] [KeyWord2] ...

    header 中の OBJECT NAXIS3(連続露出数) EXPTIME1(単一露出時間) GAIN(e/ADU なので conversion factor に相当) をリストする
    EXPTIME1 と GAIN の異なる画像を解析中に混ぜてはダメなので、良く確認しておくこと
    他にリストに入れたいキーワードがあれば幾つでも追加できる
    最後に1枚目の画像の平均値と、15000ADU 以上のカウントを持つ pixel 数が表示されるが、この値が 1000 を超えている場合はサチっている画像の可能性が高いので要注意。

  • dark 作成
    ./sigclip TRCS8桁番号.fits dark.fits [clipping 閾値]

    dark の名前は dark5.fits など、単一露出時間を付けて作る方が安全
    sigclip は連続した露出に対し各ピクセルで1次 fit を行い、fit 値からの3σ clippng 後の平均を計算する
    最も多く clip されたものを以下の様に表示(50回連続露出の場合はこれにより数枚 clipping される)
    Maximum rejected pixel: [XXX,XXX] X/XX rejected
    sigclip はこの時作成したσマップを平均0、標準偏差1に規格化して sigclip.fits として出力
    この dark 安定性マップは mask.fits として置いておき badpix 補正の参照マップとする

    mv sigclip.fits mask.fits

    連続露出枚数が10枚だと3σ ではほぼ clip されないので、30〜50枚あるのが理想

  • flat 作成
    ./sigclip TRCS8桁番号.fits flat.fits dark.fits [clipping 閾値]

    sigclip に2枚目の fits を与えると、これを dark とみなして入力画像から差し引いて解析する
    もちろん露出時間が同じ dark を与えること
    dark が与えられた場合は flat 生成目的だと判断し、出力画像の中央部(現在は x 位置は光量重心) 200x200pix 平均は1に規格化される(このモードでは 入力fits 出力fits Saturated_pix数 平均 σ clippe_pix数 を sigclip.log に出力する)
    上記平均を計算する際はσマップ(sigclip.fits)を内部で参照し、この値が3未満の部分のみを用いる
    カウントが大きくなると線形性がかなり悪化する印象なので、flat は 10000ADU 弱程度が理想

    ここまでは撮像の場合と同じ。


flat 画像(左から g,r,i、0〜2 のスケールで表示)

  • cal 作成
    ./sigclip TRCS8桁番号.fits Hg.fits
    ./sigclip TRCS8桁番号.fits Ne.fits
    ./sigclip TRCS8桁番号.fits Xe.fits
    iraf で
    imarith Hg - dark Hg
    imarith Ne - dark Ne
    imarith Xe - dark Xe
    imcomb Hg,Ne,Xe cal combine=average reject=none

    3種の輝線を全て合わせて1つの画像にする。画像内に写る輝線の本数はまちまちであるため、sigclip は dark なしモードで走らせ、その後で普通に dark を差し引いて重ね合わせる。


cal 画像(左から g,r,i)

●変形マップ作成

  • 変形補正の方針
    TriCCS はスリット分光なので、空間方向の変形トレース(slit edge の凹凸を利用)の際のコントラストがかなり悪く、KOOLS の変形解析のように空間と波長両方の情報を1枚の画像に入れると情報の分離が困難になりそうなので、縦方向と横方向の変形を別々の画像・別々のソフトで解析し、結果を合わせて用いる。画像変形は撮像解析で用いている fitawarp.c を用い、その結果を再度解析して変形マップを修正するということを繰り返して収束させる。

  • Y 方向形状評価
    ./caldy cal.fits ./

    cal.fits から各輝線の曲がりを4次関数で fit し、cal_[g/r/i].dat の情報を参照して各輝線を同定(header の DETID 情報を用いて、第2引数の path の直後に cal_[g/r/i].dat のどれかが自動選択されて文字列連結され、読み込まれる)、波長が

    波長(nm) = dW ✕ (y - Ycen) + Wcen

    で決まる y 座標に移動させるのに必要な y 方向の移動量を計算し、dy.fits として出力する。Ycen は ds9 上で 681、Wcen と dW は cal_[g/r/i].dat の先頭行に書かれた値で決まっており g,r,i それぞれで 460 0.33, 640 -0.32, 800 -0.45 である。cal_[g/r/i].dat の内容は、

    1行目:変形前に期待される y 座標の値
    2行目:空気中での輝線の波長
    3行目:輝線の種類

    4行目/5行目は変形後の位置と理想値とのずれが入っているが、元データとしてはこれらは参照しないため、無くても構わない。pix と波長の関係は3次関数で fit している。

  • X 方向形状評価
    ./caldx flat.fits

    flat.fits には、slit edge の凹凸による薄い縞が見えており、縦方向 21pix でスムージングをかけた後、この縞をトレースして横方向の曲がりを4次関数で fit し、直線にするために必要な x 方向の移動量を dx.fits として出力する。

  • 変形計算と再評価
    ./fitswarp cal.fits cal2.fits dx.fits dy.fits
    ./fitswarp flat.fits flat2.fits dx.fits dy.fits

    得られた (dx,dy) ベクトルで cal,flat それぞれを変形させる。

    mv dx.fits dx0.fits
    mv dy.fits dy0.fits
    ./caldy cal2.fits ./
    ./caldx flat2.fits
    ./vecadd dx0.fits dy0.fits dx.fits dy.fits 0.5

    変形後の画像で再度形状評価を行い、目標形状に必要な追加量を計算する。別名にしておいた1回目の結果に結果を 0.5 倍してから加算して(1倍で加算すると結構振動するので)、新しい変形マップとする。

    caldx, caldy の結果が両方同時に悪化するところまで上記手順を繰り返し、悪化したら1つ戻って最終結果とする。この手順を sh スクリプトにしたものが caldxy.sh で、

    ./caldxy.sh flat.fits cal.fits ./

    で収束するまで計算を繰り返し、最終結果を dx0.fits, dy0.fits として残す。

●flat 特性除去画像作成

  • badpix 除去と周辺部 cut
    ./fixpix flat.fits flatr.fits mask.fits flat.fits

    fixpix.c は mask.fits の値が6以上の pixel に対して、周辺部の値を距離の逆数で重みをつけて平均し、内挿する。第4引数に flat.fits を与えた場合は、カウントが負になっている pixel に対応する部分も badpix と判断して同様な処理を行なうとともに、中央幅 1/2 の範囲内での縦方向の平均値が 0.2 以下の左右部分と、この左右範囲内での横方向の平均値が 0.05 以下の上下部分を NaN で置き換える。撮像の flat ではこのような領域は周辺部には現れないが、分光の場合はかなりの面積が NaN で置き換えられる。この範囲は flat 取得の方法により変化するため、NaN 部分は切り落とさずにそのまま残すことにした。

    ./fitswarp flatr.fits flatc.fits dx0.fits dy0.fits


fixpix, fitswarp 後の flat 画像(左から g,r,i、0〜2 のスケールで表示)

  • sflat.fits 作成
    ./mksflat flatc.fits sflat.fits xy [order]

    上記画像を縦横それぞれの方向に1次元化し、それらの積をとって flat 画像を再合成する。天体画像を flat で割って fitswarp で変形後、上記 sflat を掛けることで flat の空間方向と波長方向の特性が除去されるはず。mksflat の第3引数は x/y/xy/xr/yr/xyr/xd/yd/xyd で1次元化の方向と出力を fit値とするか、入力/fit値とするか(r 付き)、入力-fit値とするか(d 付き)を表しており、第4引数を与えた場合はその際に fitting する次数となる(default は0)。画像有効部の中央 200x200 領域の平均値が保存されるように規格化して出力する。


sflat 画像(左から g,r,i、0〜2 のスケールで表示)

    mksflat.c を通す前後の画像の比を取ることで、どの程度の影響が残るのかを確認することができる。

    ./mksflat flatc.fits rflat.fits xyr


flatc/sflat 画像(左から g,r,i、0.8〜1.2 のスケールで表示)

    これらは完全に1となるのが理想的だが、slit に沿っての flat スペクトルの変化は無いと仮定すれば、変形マップの不完全による影響と考えられる。flat の明るさが大きく変化する部分では、少しのずれでも影響が出るので、g-band の下側でうねりが大きいのはそれが原因、i-band では結構迷光が入っているようなので(後述)、i-band のうねりにはその影響が入っている可能性がある。どちらにしても、標準星とターゲットがほぼ同じ位置に入っている場合は気にする必要がないが、スリット全体を使うような観測の場合には、この程度のうねりが flat に起因して入るという事になる。上記は±20% レンジでの表示だが、大部分の領域が±2% の範囲内にあり、周辺部の濃い領域がそれを超える状態となっている。

  • twilight の場合
    flat と同様に fixpix, fitswarp 処理を行った後の twilight は以下。


fixpix, fitswarp 後の twilight 画像(左から g,r,i、0~2 のスケールで表示)

    mksflat.c を通す前後の画像の比は以下の通り。


twilightc/stwilight 画像(左から g,r,i、0.8~1.2 のスケールで表示)

    g,r は lamp よりも twilight の方が slit 内でのスペクトル変化の影響が少なく、ほとんどの空間・波長方向の特性除去がうまくできている。しかし、i-band は twilight の方が圧倒的に迷光が激しく、それが slit 内でのスペクトル変化のように振舞って全く特性除去ができていない。なので、i-band の twilight は迷光が解決するまでは使わない方がいい。g,r に関しては、高空間周波数成分は lamp の flat を、全体の光量分布は twilight を使うのが良さそうなので、以下の処理を行う事にする。

    imarith twilight / flat twilightm
    median twilightm twilightm 5 5
    imarith twilightm * flat flat

    これで 5x5 以内の成分は flat 情報が優先され、それより大きい部分に関しては twilight の情報で置き換えられた flat が作成される。g,r に関して sflat まで作成し、mksflat 前後の比を調べると以下となり、想定通りの状態が確認できた。


flatc/sflat 画像(左から g,r、0.8〜1.2 のスケールで表示)

    う~ん、r-band の方に強い夜光輝線が2本あり、その輝線の跡がかなりうねっている。median smoothing のサイズが輝線幅より少し大きい程度なので、輝線の影響が中途半端に残ってしまった感じだ。median smoothing のサイズを大きくして夜光輝線の影響が無くなるところを探してみると 11pix 程度で以下のようになり、だんだんと flat lamp の波長特性によるうねりが見えてくる状態に。


flatc/sflat 画像(左から g,r、0.8〜1.2 のスケールで表示)

    ここまで来るとやはり大きいうねりのパターンが気になる。x 方向に6次で fit してうねりの成分を取り出す。i-band は flat lamp の flatc 画像を fit する。

    ./mksflat rflat.fits rflat.fits x 6


flatc/sflat 画像を6次で fit したもの(左から g,r,i、0.8〜1.2 のスケールで表示)

    ここで、sflat として flat からキャンセルすべき成分は、波長・空間両方ではなく、波長方向の情報だけにすべきである事に気が付いた。でないと slit edge のギザギザによる光量ムラが最後まで残ってしまう。なので上記 rflat での補正も含めて mksflat 部分の作業は以下となる。

    ./mksflat flatc.fits rflat.fits xyr
    ./mksflat rflat.fits rflat.fits x 6
    ./mksflat flatc.fits sflat.fits x
    imarith sflat * rflat sflat


sflat 画像(左から g,r,i、0〜2 のスケールで表示)

    twilight が無い場合の sflat 画像は以下の通り(i-band は上と同じ)


twilight 無しの場合の sflat 画像(左から g,r、0〜2 のスケールで表示)

  • cal 画像確認
    ./fitswarp cal.fits calc.fits dx0.fits dy0.fits

    caldy.c は輝線の同定結果を caldy.dat として書き出しており、これによりどの程度目標位置とずれているか確認できる。caldy.dat は cal_[g/r/i].dat と同じフォーマットで、4行目と5行目が評価結果、明るさが足りなかったり、ずれが大きくて同定ミスの可能性が高いと判断された輝線には先頭に # が付いてコメントアウトされる。caldy.dat の4vs2行目と4vs5行目をグラフにしたものが以下。


変形後の cal 画像(左から g,r,i)


pix 波長関係(左から g,r,i)

    注意すべきなのは、caldy.dat や cal_[g/r/i].dat に出てくる波長は空気中での値で、ds9 の WCS や後述の解析の最終結果などは全て真空波長であること。上記ファイルの波長よりも 0.2nm 程度長い値で表示や結果が出るので、ずれていると勘違いしないように。NIST や理科年表の値が全て空気中での波長表示になっているために、そちらとの確認の利便性を考えてこのようにしている。

●天体画像解析

  • dark 引きと flat 割り
    imarith は連続露出の画像全てに適用されるので以下で OK

    imarith TRCS8桁番号 - dark ts8桁番号
    imarith ts8桁番号 / flat ts8桁番号


天体画像例(左から g,r,i、0-100 のスケールで表示)

    ここで問題なのは、i-band の画像が明らかに flat と整合性が悪いことだ。右端のスリット幅が一部広い部分、明らかに flat では割りすぎとなっている事が分かる。この原因は、観測中の方が twilight や lamp での flat 取得時に比べて迷光の寄与が大きく、迷光部分も flat で割られているためその分凹んでしまう。とりあえず、このまま正攻法で進めてみる。

  • badpix 補正
    dark 生成時の安定性マップ(mask.fits)と flat 画像(flat.fits)を参照し、bid pixel を周辺部の値で補間、周辺部を NaN で置き換える。

    ./fixpix ts8桁番号.fits tr8桁番号.fits mask.fits flat.fits

    mask.fits が6以上の部分と、flat.fits が負の部分を補正する。
    縦方向に4pix周期で平均化した値が 2.0 を超える場合も補正(通常はC0の1列のみ適合するはず)
    上下左右の最寄りの正常値を距離の逆数の重みで平均化して代入する


fixpix 後の天体画像例(左から g,r,i、0-100 のスケールで表示)

  • 画像変形
    ./fitswarp tr8桁番号.fits tc8桁番号.fits dx0.fits dy0.fits

    dx0.fits dy0.fits の値を (x,y) とする変形ベクトルで画像変形する。

  • flat 逆補正
    imarith tc8桁番号 * sflat tc8桁番号

    sflat.fits をかけて、flat lamp の特性を除去する。iraf の imarith は NaN 値を0にしてしまうので、その後の処理が1手間増えてしまう...。下右図を見ると、flat ではわからなかった i-band 右側の迷光がかなり影響していることがわかる。撮像時にはこのような背景光の傾斜は見られなかったと思うので、スリットかグリズムが入ったことによる影響の可能性が高い。装置内で移動しているものを良くチェックしてみる必要があると思う(反射する部分は要注意)。


fixpix, sflat 補正後の天体画像例(左から g,r,i、0-100 のスケールで表示)

  • 重ねあわせ
    slit 上の明るい天体の x 座標の位置変化を追跡することで、左右方向に pixel 単位で shift しながら重ね合わせる。初回は連続露出の1枚目を基準フレームとし、直前までの積算フレームを基準フレームとして重ねあわせていく。確認のためのそれぞれの差分画像は shift_ave.fits として出力される。

    ./shift_ave tc8桁番号.fits ti8桁番号.fits

    1枚目を基準とし2枚目との残差が最小となるシフト量を探して重ねあわせる。
    ./shift_ave tc8桁番号.fits ti8桁番号.fits ti8桁番号.fits

    で、重ねあわせられた画像を基準フレームとして位置合わせをやり直す。同時に、各 pixel の σ を計算し画像に追加する。その際の計算はやや複雑で以下の通り。

    天体スペクトルの明るさは天候とは無関係に slit でどれだけケラれるかで決まるため、背景光の明るさ変化とは独立である。なので、明るい天体部にマスクをかけた後、各列横方向に6次関数で3σクリップしながら fit することで背景光と天体スペクトルを分離、天体スペクトルはそれぞれの中央部分の平均値で空間方向の強度変化に合わせて差し引き、可能な限り平坦にした状態で連続露出時の各 pixel の変化を計測、√(画像数-1)✕0.9σ を閾値としてクリッピングする。σ の値は連続露出での各 pixel 値の標準偏差で計算されるが、分光の場合は露出枚数が少ないため統計的な閾値があまり役に立たない。そのため、σ =√(偏差2乗和/(画像数-1)+pixel 値+読み出しノイズ2) でやや大き目に評価し、更に1pixel だけ単独でクリップされる事象が多い場合や、全データがクリップされてしまった pixel が存在する場合などは閾値が低すぎると判断し、閾値を10% 上げて計算をやり直し、除去し過ぎないようにするようにすることで、統計的な評価ができない部分をカバーした。また、この計算を繰り返すことでクリッピングが収束するが、第3引数を "auto" とすることで、これを行なうことにした。auto モードの場合は確認用差分画像 shift_ave.fits は生成されない。

    ./shift_ave tc8桁番号.fits ti8桁番号.fits auto

    宇宙線除去は3回目からで、2回目のみ clipping は sigclip ではなく、最大値1つ削除としている(露出枚数が少ない場合にこれが有効に働く)。3回目以降は上記の説明の通り。auto モードの最終出力の2つ目の画像は、σマップではなく、除去された成分が追加される。また、最終結果の書き出し直前で、-6σ 以下のマイナス点を中心とする -3σ までの範囲を fitting 値で置き換える。どうやら細かいゴミが検出器の直前でちょくちょく移動するようで、ある程度の間は黒い点は動かず、その内移動するという現象であるため、連続露出での違いで除去する方法は役に立たない。
    shift_ave を分光画像に適用した場合は、3つ目の画像として背景光を6次 fit で差し引いた画像が加えられる。

    しかし、i-band はひどい。スリット幅の広い部分は負のカウントとなっているので baxpix 扱いとなって連鎖的に領域が広がって除去されてしまった...


shift_ave 後の背景光差し引きされた天体画像例(左から g,r,i、-20〜20 のスケールで表示)


shift_ave により除去された成分(左から g,r,i、-20〜20 のスケールで表示)

●最終画像の作成

  • 重ねあわせ
    ./spcomb ti画像リスト 天体名.fits [order]

    ti8桁番号.fits のリストを第1引数とすることで、画像内の最も明るいスペクトルの位置を1番目の画像に合わせて重ね合わせる(なので、1枚目は S/N 最大の画像にすべき)。その際、波長方向のパターンマッチングも行い、上下方向のずれも pixel 単位でずらして重ね合わせる。本来であれば回転も補正すべきだが、sub pixel shift は避けたいので細かい調整はここではしない。装置変形により、2pixel 程度は位置ずれやねじれが発生するので、ターゲットを観測したロテータ角度と同じロテータ角度で flat や cal を個々に取得するのが理想(解析は面倒だが...)。重ね合わせは単なる加算で、選択された画像が何枚の平均画像なのかを元データの fits header を参照して確認し、枚数分の倍率をかけて足し合わせる。その上で、各画像の背景ノイズレベルとシグナルレベルを計測、S/N が最大となるように S/N/N の重みを付けて加算し、重み無しでのシグナル合計値で規格化して最終結果とする。最後に order 次で横方向に fitting を行って差し引く。default 値は 6 で shift_ave も 6 次で fit しているため、6 の場合はほぼ影響ない。最大24次まで指定できるが、多分、倍精度でも桁落ちするため(高次の項から桁落ちしにくい加算方法を取っているが...)、20次を超えたあたりから失敗する行がちらほら出てくるようになる。

    order に負の値を与えた場合は、重み1で単純に加算しこの場合は総積分時間分のカウントに相当する画像となる(即ち、-6 を入れると重み1の単純可算で、最後の fit はほぼ影響ないものとなる)。ターゲットがスリット内で一番明るい天体ではなく勝手に重みを付けられると困る場合は、予め重みを掛け算しておいてから重み無しモードで積算すべき。関与したカウント数と実測ノイズ値を合成して算出したノイズ2乗画像が2枚目に追加されている。

    手動で offset や重みを与えたい場合は、ti画像リストの第2列に x offset、第3列に y offset、第4列に weight を書き足せば OK。spcomb 第3引数に負の次数を与えると重み均一になるがそちらが優先される。

    とにかく i-band がひどいので、ちょっと変則的な方法で解析してみる。


spcomb 後の背景光差し引きされた天体画像例(左から g,r,i、-200〜200 のスケールで表示)

●空間方向の情報を最後に補正する方法

    i-band が全然ダメダメなのは、slit に沿った空間方向の flat 情報を sflat を適用する際にあえて残したために、flat とは関係ない迷光が問題を起こしている事が原因。なので、sflat では空間方向の flat 情報も一旦キャンセルし、背景光を fitting で差し引いてから補正するという方法を取ってみる(実は初めに sflat を作った際に空間情報までキャンセルしてうまく行ったことがこの手法に至ったヒントになった)。

    この場合の sflat の作り方は、

    ./mksflat flatc.fits rflat.fits xyr
    ./mksflat rflat.fits rflat.fits x 6
    ./mksflat flatc.fits sflat.fits xy
    imarith sflat * rflat sflat
    ./mksflat flatc.fits vflat.fits y
    と、3行目の部分を x fit ではなく xy fit にするだけで OK。空間情報の感度を最終的に補正するため y 方向に fit した vflat を作っておく。vflat の中央付近は1ではなくまちまちな値だが、これを用いて最終修正をするソフトの中で中央付近を1に規格化して使う予定。



sflat 画像(上段)と vflat 画像(下段) (左から g,r,i、0〜2 のスケールで表示)

    twilight が無い場合の sflat 画像は以下の通り(i-band は上と同じ)。lamp の vfalt だと端の方がかなり暗くなる(天体が増光される)ので、やはり twilight がある方がいい。



twilight 無しの場合の sflat 画像(上段)と vflat 画像(下段) (左から g,r、0〜2 のスケールで表示)

    天体画像。今度は slit 幅の広い部分がどこなのか分からない状態になっている(slit 幅による夜光の明るさ変化よりも、迷光の寄与の方が上回っている状態)。


fixpix, sflat 補正後の天体画像例(左から g,r,i、0-100 のスケールで表示)

    shift_ave で重ねて背景光を6次の fitting で除いたものが以下。


shift_ave 後の背景光差し引きされた天体画像例(左から g,r,i、-20〜20 のスケールで表示)


shift_ave により除去された成分(左から g,r,i、-20〜20 のスケールで表示)

    背景光差し引きまですると、vflat を適用していない影響で空からの背景光の面輝度ムラが見えてくる。まず、vflat で slit に沿った感度補正をし、残っている背景光のうねりを vflat を6次 fit した差分を用いて列ごとに最小2乗条件で差し引く。これで、明るさ不明の迷光の影響も差し引けるはず。spcomb では vflat の有無を確認し、vflat がある場合には各画像を読み込んだ直後に自動で上記手順が行われるようにした。


spcomb 後の背景光差し引きされた天体画像例(左から g,r,i、-200〜200 のスケールで表示)


i-band で spcomb の第3引数を 20 とした場合

●スペクトル標準星との割り算

  • 標準星の処理
    上記観測天体の処理と全く同一なので割愛

  • カタログスペクトルを用いて大気吸収を補正
    ./trcalib 天体名.fits 最終出力.fits 参照星.fits 参照星.dat > 最終出力.dat

    まず、参照星.fits とカタログスペクトルである 参照星.dat を比較し、参照星.fits 内の恒星吸収パターンから波長方向のシフト量を計算する(参照星.dat 内での #,/ を用いた吸収補正方法は KOOLS 解析の場合と同じ)。シフト後は、参照星.fits から抽出した1次元データを 参照星.dat で割り、大気吸収+装置効率カーブを算出(吸収補正前の結果は trcalib.dat として出力され、吸収補正後の結果は 最終出力.dat の第4~6列に出力される)、天体スペクトルをこれで割ることで最終スペクトルとする。また、背景光 fit していない画像を gain 値で割り、読み出しノイズの2乗を加算した後、上記効率カーブの2乗で割ることで、ノイズ2乗マップを作成、出力画像に追加する。天体スペクトル抽出の際は、このノイズマップを同じ手法で1次元し、√ を取ることでノイズレベルが評価できる(今の所はノイズに比例する値という扱い)。最も明るいスペクトル部分のみに関しては、自動的に1次元化して標準出力に書きだすので、リダイレクトでファイル出力にすれば数値データとなる。このファイルの中身は、

    1列目:波長
    2列目:天体スペクトル(比例値)
    3列目:ノイズ(実測値と比例値の合成値)
    4列目:効率+大気吸収(比例値)
    5列目:標準星スペクトル(比例値)
    6列目:標準星スペクトル(カタログ値)

    「比例値」は観測値のことで、calibration 不可能なため「比例値」としている。1次元化の際に用いた重み(psf の2倍)は、参照星のものは trcalib_psf0.dat、観測天体のものは trcalib_psf.dat と最終画像の一番上の行に出力される。

    slit 分光の場合は全カウントは信用できないと考えているため、天体スペクトルの値は適当(標準星のカタログスペクトルの単位と同じ単位とし、観測天体と標準星の spcomb 時の重み付き総露出時間比をかけたような状態になっている)。また、この結果は画像左端の1~6列目にも書き込まれている。結果は以下の通り。


trcalib 後の天体画像例(左から g,r,i、-50〜50 のスケールで表示)


trcalib 後のノイズ2乗画像例(左から g,r,i、0〜400 のスケールで表示)


左から、効率+大気吸収カーブ(4回の観測結果)、標準星の吸収補正後、観測天体スペクトル、の順

●暗い天体だと...

    slit 分光で暗い天体をちゃんと slit に入れるのは手順が確立していないと大変だが、お試しで試験観測中に quasar を1つ観測してもらったので、KOOLS との比較。

    trcalib まで終了させた画像が以下。i-band の最後の背景 fit は20次でのもの。ターゲットは中央付近に写っている2つの天体の右のもの(暗い方)。自動検出では左の天体の方を取り出してしまうため、

    ./trcalib 天体名.fits 最終出力.fits 参照星.fits 参照星.dat x座標 [ap.dat] > 最終出力.dat

    として強制的に x 位置を指定する。指定された位置±2 pix 内での最大値の場所を探してターゲットの位置とする。Aperture 内の重みを手動で与える場合は ap2.dat (2" Aperture) などを用いる。i-band では、slit が局所的に狭くなっている部分での背景うねりの方が卓越していて、ここに天体を当てないように注意が必要であることがよくわかる。どちらにしても i-band はほとんど意味が無い状態だが...


trcalib 後の天体画像例(左から g,r,i、-500〜500 のスケールで表示)

    スペクトルは以下の通り。KOOLS の結果は下図右。露出時間は TriCCS が 30sx10+100sx10=1300s、KOOLS は,が各600sx4=2400s、が600sx2=1200s。


同一天体(但し変光の可能性高)での KOOLS との比較。露出時間は上述の通り。

    まあ全体の雰囲気は合っているが、やはり slit にちゃんと乗っていない印象。この後で再度 imaging からやり直して try したようだが、そちらは明らかに slit を 2pix 外していたのでこの解析には含めていない。TriCCS 分光モードで暗い天体を観測する場合は slit が長い事を利用して最寄りの天体と同時に入るように PA を決めて同時観測するのがいいように感じた(但し、TriCCS 画像上での北方向の取り付け再現性の確認など必要)。もしくはガイダー上での offset 引き込みでもいいかも。また、i-band は迷光対策ができないと 17mag よりも暗い天体の観測を行うのは無理そうだ。

●標準星同士の比較

    HR7596 を標準星として HR5501 を解析してみた。カタログスペクトルよりもかなり赤いし、明らかに B type ではない連続光形状ででおかしい。また、r と i の接続も良くない。1次元化する際に、空間方向の平均 PSF で重みを付けて積算したので、色収差等でPSF の波長変化が大きい場合に r,i の接続不良などの問題が起こる可能性がある。これを確認するため、積算時の重み PSF の幅を2倍にして積算したものを gif アニメで比較した(下図左をクリック)。アパーチャが大きくなると全体としてまっすぐ繋がる感じになるので、明るい天体の場合は重み PSF は広げた方が良さそうだ。しかし、r,i の接続部は全く改善されず、他の原因を調査する必要がある。

    下図中は HR5501 を HR7596 を標準星として解析した最終画像と、役割を逆転させて解析した HR7596 の最終画像の比較。波長の繋がりを矢印で示している。r,i 間の段差は画像上でも確認でき、積算方法などは関係ないことがわかる。2つの星の観測時の slit 上での位置は結構離れているが、どちらも検出器上の位置としてはほぼ中央に近いため、ケラレは起こっていないはず。可能性として高そうなのは、slit 上の位置が異なるとダイクロイックを通過する角度がほんの少し異なる角度となるため、切り分けられるエネルギー比が変わってしまう可能性だ。実は、この効果のことは解析上は sflat を作る際に flatc/sflat 確認時に空間・波長方向に残った特性を6次で fit して消してしまっている。これがまずかった可能性が高そうだ ... と思ったが、実際にこの6次 fit を行わないで解析しても結果はほとんど変わらず、段差がほんの少しだけ小さくなったかなという程度だった。

    また、大気分散方向と slit の関係を調べてみると、観測時 HR5501 は ALT=46°HA=+2h、HR7596 は ALT=36°HA=-3h で slit 方向はどちらも東西と、HR7596 の方が大気分散の影響が大きいはず。また、この2つの星は、明るさもカラーもほとんど同じだが、観測 flux は HR7596 の方が HR5501 の g で1/2、r で 1/4 となっていて、カウント的には slit にちゃんと入っていないのは HR7596 の方だ。一体、何が起こっているのか...

    下図右は HR7596 による大気吸収と装置効率補正カーブで、上で示したものも同じ天体(観測日は異なる)。同一日の4連続露出でもカラーに若干のばらつきが見られたが、今回の結果も上に示したものと比べて少しだがカラーが異なる。高度がかなり低い事を考えると、多分、これが大気分散の影響でないかと思う。

    あと、下図中ではロテータの角度の違いによる波長分散方向のねじれが確認できる。g,r,i 全てで波長分散方向に対する回転方向が同一なので、多分、グリズムの重い側が常に重力的に沈む状態となっていて、全 band 独立でありながらも共通の特徴が見られるものと考えられる。


左)1次元化時の aperture 大小の影響、中)1次元化直前のスペクトル比較、右) 効率・大気吸収カーブ

    どうやら sflat が怪しいようなので、i-band で使えない twilight は敢えて使わずに、全て lamp での flat だけで再解析してみた。その結果、r,i 境界部の段差は解消されたが(下図左)、どうやら vflat が原因だったようだ。よく考えれば、slit に沿った空間方向の明るさ補正をするためのものだが、その際の slit 方向の光の当たり方の違いまでも反映されてしまうため、それが slit 幅の変化の仕方の違いのように扱われてしまったということだった。という訳で、i-band の迷光対策が完了するまでは twilight は使えないということが確定した。

    ここまで来ると、ロテータの角度変化によりスペクトルが少し傾いていることがやはり気になる。積算時の重みが中央部が1で離れるほど下がっていくので、スペクトルの傾きにより波長帯端で1次元化したときの flux が若干低下してしまう。これを補正するため、1次元化の際はスペクトルの傾きを1次関数で fit して、psf もそれに沿って生成、幅を2倍にして上記傾きに沿ってそれを適用、という方法に修正した。結果の比較が下図中。特に g-band でのスペクトルの反りが減少し、改善されたことが確認できた。暗い天体の場合はスペクトルをトレースできないので、トレース時の rms が1pix 以上ある場合は、傾きは適用せずにそのまま積算する。もともと S/N が悪いので、少し程度のスペクトルの反りは影響ないだろう(rotator の角度とスペクトルのねじれ角の関係が関数化出来た際は、fitswarp の際にその分の補正を入れられるので、この補正はほとんど無くなるはず)。

    スペクトルの傾きに沿って psf 積算するようにしたので、2つの星の psf を比較してみた。下図右はその結果(psf の幅をを2倍に広げたもの)で、HR7596 が実線、HR5501 が破線、はそれぞれ g,r,i に対応する。slit で大半がケラレている場合は、psf も若干小さくなることが予想され、HR5501 の方がやや小さいが、そもそも方向が間逆で高度も違うため、seeing の違いの方が卓越して良くわからない状態だと思う。


左)ハイブリッド vflat の影響、中)1次元化方向最適化の影響、右)psf 比較

    HR5501 のカラーが赤すぎる点について、もう少し調べてみる。この天体は 0.5sec x10, 1sec x10 で連続して観測しており、15000ADU を超える pixel も存在しない。0.5sec の結果を3倍、1sec の結果を1.5倍して個々に結果を plot したものが下図左。0.5sec の方が全体が暗くて赤い事がわかるが、数秒の間の変化なので slit への入り方が変化したとしか考えられない。下図右は HR7596 の 1sec x10, 0.5sec x10 連続観測に対し、互いに参照星として結果を出したもの。こちらは大半が同じカウントで、g だけに変化が見られる。このことから、HR7596 は g だけが slit 端にかかっているのに対し、HR5501 は全波長で slit にかかっているため、かなり赤くなったと考えるのが妥当だろう。似たような明るさの星であるにもかかわらず HR7596 の方が暗かったのは高度が低いことと天候の影響、ということになる。やはり、高度の低い天体を観測する際は、slit は天地方向に向けて観測するのがいいと思う。

●暗い天体に対する aperture 指定

    観測ターゲットが非常に暗い場合、x 位置を引数で与えてその場所から 2pix 以内での最大値の場所を中心として、ターゲットの psf の2倍の幅で重みを付けて1次元化する流れとなっているが、psf すら怪しいくらいの状況である場合は手動で psf を与えて1次元化したい状況が考えられる。その場合は、

    ./trcalib 天体名.fits 最終出力.fits 参照星.fits 参照星.dat x座標 psf.dat > 最終出力.dat

    でファイル読み込みの psf を優先させるようにした。psf.dat のフォーマットは、trcalib で裏出力される trcalib_psf.dat と同じで、

    1列目:ピーク位置を原点とする相対 x 座標
    2列目:1次元化の際の重み
    3列目:関係なし

    で、trcalib_psf.dat には 3列目には psf の実測値が入っているがこれは上記 psf.dat の入力として用いられる場合は無視される。上記の暗い天体の例の場合に適用してみたのが以下で、下左図は psf.dat 無しで trcalib を走らせた場合の trcalib_psf.dat 出力、下右図が g,i で r の trcalib_psf.dat を第6引数に与えた場合のもの。もちろん、手動で数値を修正して与えれば、任意の aperture を適用できる。


左)暗い天体の場合の aperture 例、中) r-band の aperture で統一、右)共通 aperture で1次元化

●ゴースト

    画像中心を挟んで左右対称(正確には中心対称)にゴーストが出る。グリズムの平面側が検出器側を向いているため、検出器で反射した光がグリズムの射出平面に反射して戻っていることが原因。こうならないために平面側は望遠鏡側を向いているのが望ましいのだが仕方がない。中心付近だと反射率が上り、ゴーストが強くなるので(と言っても面輝度 1/1000 程度)明るい天体と一緒に分光する際は中心対称に配置しないように注意しておいた方がいい。


ゴーストの位置

    これで気になったことは全て解析完了したので、分光解析ソフトの開発作業はここで終了。

●迷光調査

    i-band grism 使用の際の迷光を調査するため、dark 取得時に wheel に filter が入っている場合と grism が入っている場合の違いを調べてもらった。通常のコマンドで dark を取得すると、制御ソフトは勝手に分光時であっても i-band にはフィルターを入れてしまうようなので、grism 使用の場合の装置内迷光はどうなっているのかの確認が目的。dark 100sec での grism - filter の結果は以下の通り。


grism を入れた状態での dark と filter を入れた時の dark の差(左から g,r,i、0〜100 のスケールで表示)

    何と、予想とは全く異なる迷光パターンが現れた。1sec 辺りのカウント的には、r-band で <0.4ADU/sec、i-band で <0.8ADU/sec 程度で問題となっている迷光(i-band で 20ADU/sec)と比べてかなり低いので、この迷光は別物だろう。装置内にも迷光を発生させる原因が別にありそうだが、現在問題となっている迷光ではない感じだ。

    次に、夜間のドーム内を 5sec 露出で「観測」してもらった。grism を入れた状態での dark 5sec との差が以下の通り。slit 位置はモーターパルス数で 22900(in 状態), 22000, 21000, 18000, 15000, 10000, -22900(out 状態)。g,r では装置内よりもドーム内が明るく(これは正常)、i では逆に装置内の方が遥かに明るいことがわかる。g,r の表示レンジは 0-10、i の表示レンジは 0-120。


slit を動かしての dome 内分光画像(左から g,r,i、g,r は 0〜10、i は 0〜120 のスケールで表示)

    川端さんからの情報によると、TriCCS 内部の slit 周辺の構造が以下の通りで、迷光が入り放題の感じだ(せめて slit ユニットには周辺部にせり出すバッフル板が欲しい)。シャッタープレートの望遠鏡側がどういう面なのかはわからないが、見えている裏面と同様だと思えば結構反射率もありそうだ(こちらは望遠鏡側の面の塗装を極力反射率の低いものにすべき)。しかし、気になるのはドーム内の迷光が入ったとしても i-band でこれ程明るく光るとは考えにくい。やはり i-band の光源は装置内にあって、それが slit ユニットに反射して直接視野内に入ってくると考えるのが妥当だろう。slit ユニットを入れた状態で filter にして装置内を見れば、何が見えているのがイメージとして確認できる可能性がある。多分、ステージの位置検出用のフォトセンサの LED (画像から予想される波長は 850nm+1100nm の非熱的なもので人間の目には見えない波長なので気づきにくいが、多分スマホのカメラで内部を見れば光っているのがわかるはず)ではないかと思うので、露出中はセンサの通電を切る等の対策が必要になると思う。

    あと、カウントは低いが r-band で slit の位置によらず下の方に迷光が入っていることもわかる。この成分も装置内にあることは確定だと思う。grism を少しだけ動かすことが可能であれば、その位置変化の有無で grism よりも前か後かが判別できるはず。

●暗い天体その2

    再度暗い天体の分光を試してもらったが、ガイダーの調子が悪かったようで、slit viewer を見ながらガイドせずに露出したもの。

    スペクトルは以下の通り。KOOLS の結果は下図右。露出時間は TriCCS が 60sx10x2=1200s、KOOLS は,が各600sx4=2400s、が600sx6=3600s。


同一天体(但し変光の可能性高)での KOOLS との比較。露出時間は上述の通り。

    i バンドは迷光を何とかしない限り使い物にならないが、slit にちゃんと乗っていれば KOOLS よりも良い S/N が得られる可能性は感じられる。ただ、スペクトルの形状は KOOLS の赤のスペクトルに近いものにならないといけないのだが(青はより明るかった時期のスペクトル)、結構青くなってしまっている。大気分散の方向が slit と並行でなく、青い側に slit が当たっているか、標準星の方で逆に赤い側に slit が当たっている可能性が高いと考えられる(object は EL=60°、標準星は EL=50°)。

●スペクトルを1つに統合

    gri のスペクトルを1つに統合し、任意の bin のスペクトルとして再合成する。

    ./onespec 最終出力.dat step > 統合最終出力.dat

    onespec は、C0, C1, C2 にある 最終出力.dat を順次読み取って 0.01nm ステップのデータに spline 補間で変換、0.01nm 以上の任意の step のデータに変換してから接続、1つのデータとして表示する。接続部分は、ノイズレベルが 1:3 ~ 3:1 の範囲内でノイズの2乗の逆数の重みで統合している。以下のような感じ。これを行うと1step が1pix ではなくなるので、ノイズレベルが1単位情報あたりでなくなることに注意。step を 0.25 にすると KOOLS の VB/VR との比較ができる。

●解析手順まとめ

    1. 解析 directory に全 c ソースを持って来てコンパイル

      gcc headerinfo.c -o headerinfo -lcfitsio -Wall
      gcc sigclip.c -o sigclip -lm -lcfitsio -Wall
      gcc caldx.c -o caldx -lm -lcfitsio -Wall
      gcc caldy.c -o caldy -lm -lcfitsio -Wall
      gcc fixpix.c -o fixpix -lcfitsio -Wall
      gcc fitswarp.c -o fitswarp -lm -lcfitsio -Wall
      gcc mksflat.c -o mksflat -lm -lcfitsio -Wall
      gcc shift_ave.c -o shift_ave -lm -lcfitsio -Wall
      gcc spcomb.c -o spcomb -lm -lcfitsio -Wall
      gcc trcalib.c -o trcalib -lm -lcfitsio -Wall

      gcc vecadd.c -o vecadd -lm -lcfitsio -Wall

    2. 観測日 directory 作成、データ転送
      mkdir 2024MMDD
      scp ... 2024MMDD/

      CMOS 毎に directory を分ける
      cd 2024MMDD
      mkdir C0
      mkdir C1
      mkdir C2
      mv *0.fits C0/
      mv *1.fits C1/
      mv *2.fits C2/

    3. iraf & ds9 起動(以下、cl シェル内での作業)
      cl
      !ds9 &
      cd 2024MMDD/C0

    4. dark と mask 作成
      #dark は 5sec 露出の場合、
      !../../headerinfo ./
      !../../sigclip TRCS8桁番号.fits dark5_1.fits
      !mv sigclip.fits mask_1.fits
      !../../sigclip TRCS8桁番号.fits dark5_2.fits
      !mv sigclip.fits mask_2.fits
      !../../sigclip TRCS8桁番号.fits dark5_3.fits
      !mv sigclip.fits mask_3.fits
      imdel dark5,mask
      imcomb dark5_?.fits dark5 combine=average reject=none
      imcomb mask_?.fits mask combine=average reject=none
      #他の露出時間の dark も同様に準備しておく。mask は最も枚数の多いもので準備。
      #(10sec 以上の長時間 dark では mask はつくらない方がいいかも...)。
      displ dark5 1
      displ mask 1 zs- zr- z1=6 z2=7
      #dark 0.1sec, 0.5sec など
      !../../sigclip TRCS8桁番号.fits dark01.fits
      !../../sigclip TRCS8桁番号.fits dark05.fits

    5. flat と cal 作成
      #flat は 0.5sec、cal は 0.1sec 左中右の3ヶ所露出であるとして、
      !../../sigclip TRCS8桁番号.fits _flat_1.fits dark5.fits
      !../../sigclip TRCS8桁番号.fits _flat_2.fits dark5.fits
      !../../sigclip TRCS8桁番号.fits _flat_3.fits dark5.fits
      imdel _flat
      imcomb _flat_?.fits _flat combine=average reject=minmax nlow=2 nhigh=0
      displ _flat 1 zs- zr- z1=0 z2=2
      #twilight がある場合 (露出 0.1sec として)
      !../../sigclip TRCS8桁番号.fits _twilight_1.fits dark01.fits
      !../../sigclip TRCS8桁番号.fits _twilight_2.fits dark01.fits
      imdel _twilight
      imcomb _twilight_?.fits _twilight combine=average reject=none
      displ _twilight 1 zs- zr- z1=0 z2=2
      #_twilight の方が _flat よりも S/N が高い場合
      !mv _twilight.fits _flat.fits
      #_flat の方が _twilight よりも S/N が高い場合
      imdel _twilightm
      imarith _twilight / _flat _twilightm
      median _twilightm _twilightm 11 11
      !../../fixpix _twilightm.fits _twilightr.fits mask.fits _flat.fits
      imarith _twilightr * _flat _flat
      displ _flat 2 zs- zr- z1=0 z2=2

      !../../sigclip TRCS8桁番号.fits _Hg_1.fits
      !../../sigclip TRCS8桁番号.fits _Hg_2.fits
      !../../sigclip TRCS8桁番号.fits _Hg_3.fits
      imdel _Hg.fits
      imcomb _Hg_?.fits _Hg.fits combine=average reject=minmax nlow=2 nhigh=0
      imarith _Hg - dark01 _Hg
      !../../sigclip TRCS8桁番号.fits _Ne_1.fits
      !../../sigclip TRCS8桁番号.fits _Ne_2.fits
      !../../sigclip TRCS8桁番号.fits _Ne_3.fits
      imdel _Ne.fits
      imcomb _Ne_?.fits _Ne.fits combine=average reject=minmax nlow=2 nhigh=0
      imarith _Ne - dark01 _Ne
      !../../sigclip TRCS8桁番号.fits _Xe_1.fits
      !../../sigclip TRCS8桁番号.fits _Xe_2.fits
      !../../sigclip TRCS8桁番号.fits _Xe_3.fits
      imcomb _Xe_?.fits _Xe.fits combine=average reject=minmax nlow=2 nhigh=0
      imarith _Xe - dark01 _Xe
      imdel _cal
      imcomb _Hg,_Ne,_Xe _cal combine=average reject=none
      displ _cal 1

      #分光データは先頭に "_" を付けて撮像データと識別する事にする。
      #dark や mask などは撮像と分光で共通なので、これらには "_" は付けない。

    6. 変形マップと sflat 画像準備
      !../../caldxy.sh _flat.fits _cal.fits ../../
      cat caldy.dat
      !gnuplot ../../caldy.plt &
      # ↑fitting 残差をグラフで確認(caldy.dat の 4,5 列目を plot しただけ)
      #最終的な dysg, dxsg 値はどちらも 0.4pix 以下になるはず...
      #cal 輝線の変形後の位置と補正形状を確認したい場合は以下4行
      imdel _calc,_calb
      !../../fitswarp _cal.fits _calc.fits dx0.fits dy0.fits
      blkav _calc _calb 20 1
      displ _calb 1

      !../../fixpix _flat.fits _flatr.fits mask.fits _flat.fits
      !../../fitswarp _flatr.fits _flatc.fits dx0.fits dy0.fits
      !../../mksflat _flatc.fits _rflat.fits xyr
      !../../mksflat _rflat.fits _rflat.fits x 6
      !../../mksflat _flatc.fits _sflat.fits xy
      imarith _sflat * _rflat _sflat
      !../../mksflat _flatc.fits _vflat.fits y
      displ _sflat 1 zs- zr- z1=0 z2=2
      displ _vflat 1 zs- zr- z1=0 z2=2

    7. dark 引きと _flat 割り
      !../../headerinfo ./
      #露出時間毎に異なるリストにする(以下、0.5,1,2,5,10,30sec の露出がある場合)
      !../../headerinfo ./ > tmp
      !cut -c40-47 tmp > tmp1
      !paste tmp1 tmp | sort -n | cut -f2 | cut -c-17 > TRCS0
      !paste tmp1 tmp | sort -n | cut -f2
      !cp TRCS0 TRCS05 && emacs TRCS05 &
      !cp TRCS0 TRCS1 && emacs TRCS1 &
      !cp TRCS0 TRCS2 && emacs TRCS2 &
      !cp TRCS0 TRCS5 && emacs TRCS5 &
      !cp TRCS0 TRCS10 && emacs TRCS10 &
      !cp TRCS0 TRCS30 && emacs TRCS30 &
      #など
      #全ての編集が終了したら
      !sed -e 's/TRCS/_ts/g' TRCS05 > _ts05
      !sed -e 's/TRCS/_ts/g' TRCS1 > _ts1
      !sed -e 's/TRCS/_ts/g' TRCS2 > _ts2
      !sed -e 's/TRCS/_ts/g' TRCS5 > _ts5
      !sed -e 's/TRCS/_ts/g' TRCS10 > _ts10
      !sed -e 's/TRCS/_ts/g' TRCS30 > _ts30
      imdel @_ts05
      imarith @TRCS05 - dark05 @_ts05
      imarith @_ts05 / _flat @_ts05
      imdel @_ts1
      imarith @TRCS1 - dark1 @_ts1
      imarith @_ts1 / _flat @_ts1
      imdel @_ts2
      imarith @TRCS2 - dark2 @_ts2
      imarith @_ts2 / _flat @_ts2
      imdel @_ts5
      imarith @TRCS5 - dark5 @_ts5
      imarith @_ts5 / _flat @_ts5
      imdel @_ts10
      imarith @TRCS10 - dark10 @_ts10
      imarith @_ts10 / _flat @_ts10
      imdel @_ts30
      imarith @TRCS30 - dark30 @_ts30
      imarith @_ts30 / _flat @_ts30
      !ls _ts*.fits > _ts
      #これ以降は異なる露出時間のものも一緒に解析できる

    8. badpix 補正と周辺部 NaN 置換
      !sed -e 's/^_ts/_tr/g' _ts > _tr
      !paste -d" " _ts _tr | sed -e 's/^/..\/..\/fixpix /g' | sed -e 's/$/ mask.fits _flat.fits/g' > tmp.sh
      !sh tmp.sh
      imexam @_tr 1

    9. 画像変形と sflat 補正
      !sed -e 's/^_tr/_tc/g' _tr > _tc
      !paste -d" " _tr _tc | sed -e 's/^/..\/..\/fitswarp /g' | sed -e 's/$/ dx0.fits dy0.fits/g' > tmp.sh
      !sh tmp.sh
      imarith @_tc * _sflat @_tc
      imexam @_tc 1

    10. 重ね合わせと CR 除去
      !sed -e 's/^_tc/_ti/g' _tc > _ti
      !paste -d" " _tc _ti | sed -e 's/^/..\/..\/shift_ave /g' | sed -e 's/$/ auto/g' > tmp.sh
      !sh tmp.sh
      !sed -e 's/$/[*,*,3]/g' _ti > _ti3
      !sed -e 's/$/[*,*,2]/g' _ti > _ti2
      #背景光差し引きされた画像の確認
      imexam @_ti3 1
      #除去された成分の確認
      imexam @_ti2 1

    11. 重ねあわせ
      #天体毎に重ねたい画像(_ti8桁番号.fits)のリストを作る
      !emacs _OBJECT &
      !../../spcomb _OBJECT _OBJECT.fits
      displ _OBJECT.fits 1
      #i-band で背景のうねりが目立つ場合
      #!../../spcomb _OBJECT _OBJECT.fits 20
      #重み無し(6次 fit)モードで走らせる場合
      #!../../spcomb _OBJECT _OBJECT.fits -6
      #天体が暗いなど手動で offset を与えたい場合は _OBJECT 内で第2列に x offset を、
      #必要であれば第3列に y offset を書き足して下さい。更に重みも手動で与えたい場合は
      #第4列に weight を書き足して下さい(重み均一であれば、上記のように第3引数に
      #負の値を与えればそちらが優先されます)。
      #標準星も同様
      !emacs _STAR &
      !../../spcomb _STAR _STAR.fits

    12. 効率と大気吸収補正
      #標準星のカタログスペクトルを STAR.dat として1つ上(観測日)の directory に保存
      #STAR.dat の第1列は波長(Å or nm or μm)、第2列は flux density であること
      #カタログスペクトル中の星の吸収部分は / で、大気の吸収部分は # でコメントアウト
      !../../trcalib _OBJECT.fits _OBJECT_.fits _STAR.fits ../STAR.dat > _OBJECT_.dat
      #_OBJECT.dat のフォーマットは以下の通り
      #column 1: 波長 nm
      #column 2: 観測天体の最終スペクトル(比例値)
      #column 3: 観測天体の最終スペクトル中のノイズレベル
      #column 4: 効率+大気吸収(比例値)
      #column 5: 標準星観測値
      #column 6: 標準星スペクトル(カタログ値)
      #確認用に出力される trcalib_psf.dat のフォーマットは以下の通り。
      #column 1: ピーク位置を原点とする相対 x 座標
      #column 2: 1次元化の際の重み(次列にある実測 psf の2倍の幅)
      #column 3: psf 実測値

      #OBJECT よりも明るい天体が slit 内にあるか、OBJECT が暗くて関係ない場所が選択される場合
      !../../trcalib _OBJECT.fits _OBJECT_.fits _STAR.fits ../STAR.dat x位置 > _OBJECT_.dat
      #更に1次元化の際の重みも指定する場合(ap2.dat は上記 trcalib_pdf.dat と同じで2列目まで有効)
      !../../trcalib _OBJECT.fits _OBJECT_.fits _STAR.fits ../STAR.dat x位置 ap2.dat > _OBJECT_.dat

      #gri スペクトルを1つに統合
      cd ..
      !../onespec _OBJECT_.dat 0.25 > _OBJECT_T_.dat
      #KOOLS VB と合わせる場合
      !../onespec _OBJECT_.dat 0.25 364.25 876 | tac > _OBJECT_K_.dat
      #KOOLS VR と合わせる場合
      !../onespec _OBJECT_.dat 0.25 514.25 1026 | tac > _OBJECT_K_.dat


iwamuro@kusastro.kyoto-u.ac.jp