岩室 史英 (京大宇物)

●ソース(随時更新) headerinfo.c sigclip.c boxcar.c fixpix.c shift_ave.c trpar phot.cl fitswarp.c convcoo.c convphot.c mktrpar.c chkcoo.c fill99.c chkcoo.plt imcomb.c 上記全部 改変履歴： 2021.09.30: とりあえず完成 2021.10.01: ３色位置合わせと測光スクリプト追加 2021.10.02: mktrpar.c の最後、asin → deg の変換忘れを修正 2021.10.03: fixpix.c,shift_ave.c で WCS を引継ぐよう修正、convphot.c で x,y → RA,Dec 変換 2021.10.05: chkcoo.c でカタログとの比較照合ができるようにした(convphot.c, shift_ave.c 修正) 2021.10.06: chkcoo.c の収束条件を修正 2021.10.07: カタログ欠損を 99.99.. で埋める fill99.c を作成、chkcoo.c, fitswarp.c マイナー修正 2021.10.08: convcoo.c で明るい天体中心で複数ピークがある場合に１個に統合するよう修正 2021.10.09: chkcoo.c でカタログとのずれ量(dx,dy)を pix で出すよう修正 2021.10.09: convphot.c で座標算出時に歪曲補正するよう修正、fixpix.c,shift_ave.c マイナー修正 2021.10.10: chkcoo.c,convphot.c で歪曲の式を変更、fitswarp.c で画像の歪曲補正もするよう修正 2021.10.11: chkcoo.c,fitawarp.c で位置合わせ画像のカウントを 1ADU=25mag となるよう修正 2021.10.12: fitswarp.c で CRPIX1,CRPIX2 値を crpix.dat からの読み込みに変更 2021.10.16: chkcoo.c で高緯度域での座標歪みに対応、convphot.c, mktrpar.c, fitswarp.c も修正 2021.10.17: phot.cl 微修正(コメントアウト部分の修正なので多分関係なし) 2021.11.16: chkcoo.c の最後に gnuplot で位置ずれと等級ずれを chkcoo.plt で図示するよう追加 2021.12.07: chkcoo.c で mag-error 図示追加、chkcoo.plt,fixpix.c,shift_ave.c,phot.cl も修正 2021.12.29: sigclip.c で twilight に対応するため１次 fit で clip するよう修正 2021.12.31: phot.cl で daofind のパラメータ変更 2021.12.31: fixpix.c で縦方向に４pix 周期で出る異常値に対する修正を追加 2022.01.02: shift_ave.c で周辺部のデータのない領域を０から nan に変更 2022.01.02: fitswarp.c で歪曲補正オプションを付けた場合に背景光レベルも０にするよう変更 2022.01.02: WCS 位置を初期値として位置合わせ＋平均化を行う imcomb.c を作成 2022.01.02: sigclip.c で１フレームしかない場合に終了するよう微修正 2022.02.24: convphot.c, chkcoo.c で背景ノイズから１σ限界等級も評価するように微修正 2022.03.07: sigclip.c で、dark ありの場合に規格化前の情報をログに残すよう修正 2022.03.07: fixpix.c で、4pix 周期の bad column 検知レベルを 1.5 から 2.0 に修正 2022.03.07: chkcoo.c で、star A マークを使用する場合のバグを修正 2022.04.10: sigclip で一定値しか出さない pix での mask 値を nan から 1000 に変更 2022.04.10: boxcar で smoothing 前に 3x3 median smoothing を追加 2022.04.10: fixpix で bad column の条件を mask 値平均 > 2 から mask >2 の pix 数に変更 2022.08.31: sigclip.c で flat 生成の際の規格化領域を中央 200x200 pix に変更 2022.08.31: fixpix.c で flat のカウントが負の pixel も周辺値で内挿するよう修正 2022.08.31: shift_ave.c で急激な PSF 変化の成分を CR と誤検知しないよう修正 2022.08.31: fitswarp.c でヘッダ情報部分に微修正 2022.09.13: headerinfo.c で１枚目の平均値と 15000ADU 超え pix 数を表示するよう修正 2023.01.30: chkcoo.plt のラベル名(magerror)を修正 2023.01.30: chkcoo.c で、位置合わせ天体自動検出の方法を修正 2023.01.30: fitswarp.c で撮像画像変形の際の出力画像を naxis=2 にするよう修正 2023.11.17: chkcoo.c で、CRVAL1,2 を１枚目での値で統一するよう修正 実際の解析手順は一番下にあります ●ポリシー dark と flat の連続露出の重ね合わせは sigclip で行う 高空間周波数成分は twilight or domeflat で、低周波成分は sky flat で補正する Badpix は dark の安定性と flat のカウントで自動判別する 連続露出の天体画像は pixel 単位で位置合わせして重ね合わせ、CR 除去も行う 回転や 1pixel 未満の shift は行わない(使い勝手重視の最終画像のみ例外) ●下準備 header 情報確認 ./headerinfo ディレクトリ [KeyWord1] [KeyWord2] ... header 中の OBJECT NAXIS3(連続露出数) EXPTIME1(単一露出時間) GAIN(e/ADU なので conversion factor に相当) １枚目の平均値 15000ADU超えpix数 をリストする EXPTIME1 と GAIN の異なる画像を解析中に混ぜてはダメなので、良く確認しておくこと 他にリストに入れたいキーワードがあれば幾つでも追加できる dark 作成 ./sigclip TRCS8桁番号.fits dark.fits [clipping 閾値] dark の名前は dark5.fits など、単一露出時間を付けて作る方が安全 sigclip は連続した露出に対し各ピクセルで１次 fit を行い、fit 値からの３σ clippng 後の平均を計算する 最も多く clip されたものを以下の様に表示(50回連続露出の場合はこれにより数枚 clipping される) Maximum rejected pixel: [XXX,XXX] X/XX rejected sigclip はこの時作成したσマップを平均０、標準偏差１に規格化して sigclip.fits として出力 この dark 安定性マップは mask.fits として置いておき badpix 補正の参照マップとする mv sigclip.fits mask.fits 　　　　　　　dark 画像　　　　　規格化されたσマップ(６以上を表示) 連続露出枚数が10枚だと３σ ではほぼ clip されないので、30〜50枚あるのが理想 flat 作成 ./sigclip TRCS8桁番号.fits flat.fits dark.fits [clipping 閾値] sigclip に２枚目の fits を与えると、これを dark とみなして入力画像から差し引いて解析する もちろん露出時間が同じ dark を与えること dark が与えられた場合は flat 生成目的だと判断し、出力画像の中央部 200x200pix 平均は１に規格化される (このモードでは 入力fits 出力fits Saturated_pix数 平均 σ clippe_pix数 を sigclip.log に出力する) 上記平均を計算する際はσマップ(sigclip.fits)を内部で参照し、この値が３未満の部分のみを用いる カウントが大きくなると線形性がかなり悪化する印象なので、flat は 10000ADU 弱程度が理想 (以下の例は平均が 15000ADU 近くあるので良くないが...) flat 画像 (0.95-1.05表示) Sky flat 作成 観測中の背景光から Sky flat を生成し、低空間周波数成分の flat として用いる flat 生成と同様に全観測天体画像で個々に sigclip を通す ./sigclip TRCS8桁番号.fits sky1.fits dark.fits ./sigclip TRCS8桁番号.fits sky2.fits dark.fits ... imcomb sky1,sly2,... sky combine=median imarith sky / flat sky imedit sky skym 星が写っている部分はσマップが大きい値となるので、星のない部分の平均で規格化される このため、多数の画像を median で重ね合わせることで median sky が生成できる 望遠鏡の pointing 位置には常に天体が入るので、flat で割った後それを imedit で手作業で消す imedit のパラメータは (radius = 50.) Substitution radius (search = 5.) Search radius (buffer = 10.) Background buffer width (width = 20.) Background width という感じで、"b" キーで半径50の範囲が周辺部のノイズモデルで差し替えられる median 画像(0.7〜1.3表示)　　　　imedit で補正後　　　 この画像に boxcar smoothing をかける flat に必要な明るさの 1/100 程度の明るさしかないので、100x100 サイズで平均化 ./boxcar skym.fits skyb.fits 100 100 上下端のケラレと思われる明るさ変化は 100pix よりも小さい変化なので形状が変わってしまう この boxcar では周辺部でカーネル幅を２次関数的に減少させることでこれに対応している 1 pixel に非常に大きなノイズ値が入っていると 100x100 に広がってしまうので 3x3 median smoothing を後に追加した。 smoothing 前の画像を後の画像で割ることで、高周波成分のみを取り出して flat を確認する imarith skym / skyb skyh smoothing 画像　　　　　　　smoothing 前/後 median sky では多くの天体の影響があるが、smoothing で flat への影響は 0.1% 以下になる 上右画像を拡大すると細かい模様が見えるので、flat がうまくできていないこともわかる まあ、仕方がないので両方の情報を合成して最終 flat とする imarith skyb * flat skyflat sky x flat 合成画像(左から g,r,i 0.7〜1.3表示) flat 中央付近に結構な傾斜(上下方向に 10% 程度)があるのが気になる ダイクロイックへの入射角の視野位置による変化(42°〜48°)で g,r の切り分け波長が変化？ ダイクロイックのコーティング材質の屈折率を 1.5 程度だと仮定すると、層内での角度は 26.5°〜29.7°となるので、cos26.5°/cos29.7°=1.03 で3% の波長シフトが考えられる(cos で透過波長が変化する原理に関してはこのページの「ファブリぺロー」の説明を参照)。g,r はそれぞれ波長の 25% 程度のバンド幅なので、3% の切り分け波長の移動で 10% 程度の傾きは説明できる。r バンドは２枚のダイクロイックが逆方向に配置されているので、この影響が２倍で現れて傾きがやや大きいことが考えられる(実際は少し傾きが急な程度なので、これほどの影響は出ていない感じだが)。 → これに関しては、後にバンドパス filter が追加されて、視野内での波長範囲の変化がなくなったため、背景光の傾斜は解消された。 ●天体画像解析 dark 引きと flat 割り imarith は連続露出の画像全てに適用されるので以下で OK imarith TRCS8桁番号 - dark ts8桁番号 imarith ts8桁番号 / skyflat ts8桁番号 badpix 補正 dark 生成時の安定性マップ(mask.fits)を参照し、bid pixel を周辺部の値で補間する ./fixpix ts8桁番号.fits tr8桁番号.fits mask.fits [flat.fits] mask.fits が６以上の部分と、flat.fits が負の部分を補正する。 縦方向に４pix周期で平均化した値が 2.0 を超える場合も補正(通常はC0の１列のみ適合するはず) 上下左右の最寄りの正常値を距離の逆数の重みで平均化して代入する 重ねあわせ 連続露出隣同士のパターンマッチングにより、pixel 単位で shift しながら重ね合わせる ./shift_ave tr8桁番号.fits tc8桁番号.fits １枚目を基準とし２枚目との残差が最小となるシフト量を探す 天候変化に対応するため、基準フレームには最小２乗条件で定数倍＋定数加算をして残差を出す 計算時に表示される値は、サブフレーム番号 xshift yshift 残差σ 定数倍値 定数加算値 の順 これにより、生成された画像を新たな基準画像として再度重ね合わせる ./shift_ave tr8桁番号.fits tc8桁番号.fits tc8桁番号.fits 重ね合わせの際、同じ位置でのカウント値の平均の他にσを以下で計算し出力画像に追加する √(残差２乗和/(データ数－1)＋ピクセル値＋１) 次に走らせた際には、(√(フレーム数-1)x0.90) xσ　を上回る値のデータは異常値として clip する (閾値の下限は 1.8σ、上限は 3.5σ) ./shift_ave tr8桁番号.fits tc8桁番号.fits tc8桁番号.fits ./shift_ave tr8桁番号.fits tc8桁番号.fits tc8桁番号.fits ... 何度も繰り返すことで clip される数が収束し、宇宙線他の突発的イベントをカットできる また、shift_ave.fits には各画像から fit 後の基準画像を差し引いた確認用画像が出力される 第３引数を "auto" とすることで、上記繰り返し計算が自動で行われる この場合は確認用差分画像 shift_ave.fits は生成されない ./shift_ave tr8桁番号.fits tc8桁番号.fits auto 宇宙線除去は３回目からで、２回目のみ clipping は sigclip ではなく、最大値１つ削除 としている(露出枚数が少ない場合にこれが有効に働く)。３回目以降は上記の説明の通り。 auto モードの最終出力では、σマップではなく除去された成分が出力画像に追加される。 また、PSF が大きく変化すると該当部分のカウントが大きく変化し、場合によっては CR の 判定条件を満たしてしまう。その場合は、明るい星に対し同じ側に一律に除去領域が発生 するため、これと位置合わせ基準画像でパターンマッチを行い、100σ以上の明るい部分と 10ピクセル以上一致していれば PSF 変化が起こっているものと判断し、一致部分は除去 せずに重ね合わせる。 CR 除去の様子 flat が良くない様子がまるわかりだが、その内差し替え予定... 天体画像例(左から g,r,i) ●３色位置合わせと測光 g,i を r のサイズと向きに揃え(１回目のみ初期値要)、３色の平均位置に移動(r のみ整数値シフト) 位置合わせに必要な平行移動、回転、拡大は１手順で行う(pixel 情報の混ざりあいを最小にするため) 位置を合わせたら３色を平均して gri 画像を作り、天体位置テーブル作成 天体位置テーブルの座標を逆変換して g,r,i のそれぞれの初期位置の座標に合わせる それを元に g,r,i 各画像で測光、求められた新たな位置を用いて g,r,i の位置関係を再計算 上記 1. から繰り返す 平行移動、回転、拡大 ./fitswarp tc8桁番号.fits tc8桁番号r.fits 回転角 拡大率 中心x 中心y 移動x 移動y 画像変形ソフトをベースに作成しているものなので、与えられた平行移動、回転、拡大量から各 pixel の移動に必要なベクトルマップを生成し、１手順で変形を行うことができる。trpar に入れられている初期値を使ってとりあえず位置合わせする。 天体位置検出 imcomb で平均した gri 画像(tmp1.fits) に daofind を使う daofind tmp1 output=tmp1.coo fwhmpsf=6 sigma=10 verify- fwhmpsf と sigma は seeing や background の明るさにより調整した方がいいかも 表示して検出天体を目で確認 displ tmp1 1 zs+ tvmark 1 tmp1.coo daofind の検出天体確認 天体位置テーブル逆変換 ./convcoo tmp1.coo cooファイルリスト trpar cooファイルリストは g,r,i それぞれのディレクトリを含むファイル名３つで C0/tcXXXXXXX0.coo C1/tcXXXXXXX1.coo C2/tcXXXXXXX2.coo trpar に入っているパラメータを参照して tmp1.coo の座標を逆変換して書き出す 測光 変換された座標リストを用いて phot で測光する phot C0/tcXXXXXXX0.fits coords=C0/tcXXXXXXX0.coo output=C0/tcXXXXXXX0.phot inter- verify- fwhmpsf=6 aperture=6 annulus=12 dannulus=12 zmag=27.5 phot C1/tcXXXXXXX1.fits coords=C1/tcXXXXXXX1.coo output=C1/tcXXXXXXX1.phot inter- verify- fwhmpsf=6 aperture=6 annulus=12 dannulus=12 zmag=27.5 phot C2/tcXXXXXXX2.fits coords=C2/tcXXXXXXX2.coo output=C2/tcXXXXXXX2.phot inter- verify- fwhmpsf=6 aperture=6 annulus=12 dannulus=12 zmag=27.5 測光すると天体位置も修正されるので、g,r,i での結果を用いて再度変形パラメータ trpar を更新 ./convphot C0/tcXXXXXXX0.fits C0/tcXXXXXXX0.phot C0/tcXXXXXXX0.coo ./convphot C1/tcXXXXXXX1.fits C1/tcXXXXXXX1.phot C1/tcXXXXXXX1.coo ./convphot C2/tcXXXXXXX2.fits C2/tcXXXXXXX2.phot C2/tcXXXXXXX2.coo ./mktrpar cooファイルリスト > trpar convphot は phot の結果から x,y,mag,err を抜き出し、coo ファイルを上書き更新する (Sky 部分のノイズから、１σ限界等級も計算し５列目に追加するように機能追加した) 先頭行には # CRVAL1 CRVAL2 CRPIX1 CRPIX2 CD_1_1 CD1_2 CD2_1 CD2_2 DETID の数値が書き出される mktrpar は３つの coo ファイルから相対位置関係を計算し trpar を更新、 (mktrpar.log に３色の回転と拡大にパラメータのログを追加で残す) この部分は、cl のコマンドリスト phot.cl として準備してあり、内部ファイル番号を置き換えて cl < コマンドリスト.cl でまとめて走らせる 4～5回ほど連続で走らせると g,r,i の相対位置関係が大体収束する ３色位置関係の決定精度 以下は同一日に取得された11天体画像それぞれでの g,i を r に合わせるためのパラメータ r バンドは平行移動 (+8,-9) で整数の一定値 番号 g 回転 g 拡大 g 平行移動 i 回転 i 拡大 i 平行移動 INR2-STR 1 +0.202800 1.000640 (-7.792620 +1.458662) +0.253181 0.999717 (-0.268955 +7.534008) 35.366083 2 +0.202504 1.000418 (-7.756078 +1.297923) +0.252048 0.999643 (+0.111387 +7.380376) 23.037639 3 +0.215353 1.000384 (-8.036901 +0.202979) +0.247895 1.000158 (+0.311491 +7.489306) 54.345250 4 +0.206018 1.000354 (-7.712763 +1.228003) +0.253555 0.999579 (+0.373543 +7.978238) 48.116472 5 +0.203450 1.000349 (-7.383922 +1.564718) +0.254244 0.999735 (+0.963462 +7.965909) 68.125667 6 +0.194424 1.000400 (-7.839160 +1.573114) +0.250344 0.999824 (+0.260359 +7.551385) 41.317778 7 +0.198167 1.000109 (-7.956742 +2.500413) +0.248191 0.999668 (+0.275398 +7.329146) 45.488056 8 +0.205979 1.000401 (-7.863840 +2.206353) +0.253551 0.999663 (+0.333951 +7.441117) 54.644944 9 +0.204994 1.000268 (-7.627034 +0.216021) +0.509167 0.999477 (+1.693751 +11.558049) 224.62583 10 +0.204652 1.000209 (-7.915420 -0.256020) +0.439682 0.999682 (+0.325836 +9.669316) 79.479222 11 +0.217733 1.000662 (-7.892544 +0.266374) +0.438091 0.999891 (+0.151074 +10.239834) 73.069333 12 +0.206756 1.000416 (-8.108623 -0.011536) +0.449184 0.999866 (+0.039817 +10.100443) 71.293250 装置が一度逆さまになって以降、位置関係が不安定になっている。特に i バンドは望遠鏡から一番離れているためか影響が大きい。g,r,i の位置関係は毎回ちゃんと計測が必要ということが確認できた。 カラー画像 (歪曲補正後のカラー画像はもう少し下参照) tc8桁番号r.fits は位置合わせが済んでいるので、rgbsun で合体させるだけ 画像のメジアン値(m2,m1,m0 とする)を調べておいてそれを０に設定するのが簡単 imstat C2/tcXXXXXXX2r.fits,C1/tcXXXXXXX1r.fits,C0/tcXXXXXXX0r.fits fields=midpt color rgbsun C2/tcXXXXXXX2r.fits C1/tcXXXXXXX1r.fits C0/tcXXXXXXX0r.fits tcXXXXXXX.ras sw+ rz1=m2 rz2=m2+25 gz1=m1 gz2=m1+25 bz1=m0 bz2=m0+25 !gimp tcXXXXXXX.ras & ３色合成例 ●他のカタログとの比較照合 カタログの取得 VizieR から SDSS カタログなど、該当領域のカタログを持ってくる 手順は「解析手順まとめ」の最後を参照 バンド毎に位置合わせ 同じ領域の SDSS のカタログを tcXXXXXXX.sdss として、 ./chkcoo C0/tcXXXXXXX0.coo tcXXXXXXX.sdss 5 > C0/tcXXXXXXX0.sdss ./chkcoo C1/tcXXXXXXX1.coo tcXXXXXXX.sdss 6 > C1/tcXXXXXXX1.sdss ./chkcoo C2/tcXXXXXXX2.coo tcXXXXXXX.sdss 7 > C2/tcXXXXXXX2.sdss chkcoo は WCS に基づく RA Dec を SDSS などのカタログの RA Dec と比較照合して、位置・拡大率・角度を合わせる(ここの合わせ方はその後、カタログの RA,Dec を検出器上の x,y 座標に変換して合わせる方法に変更、それに伴い以下の dRA,dDec は全て dX dY となった)。その際、等級情報(chkcoo の第３引数の数字は等級を参照するカタログの列番号)も比較し、異常な明るさのずれを示すものはペアリングを間違えているか、明る過ぎるなどの理由で重心検出がうまくできていないものと判断して外し、次に位置ずれが他の天体の位置ずれ分布と比較して大きものを位置合わせの条件から外す。どちらも √(全体数)/2 の値を下限値 1.5 上限値 3.0 で置き換え、この値を σ の閾値としている。 chkcoo で重要になるのは初期パラメータの決定で、これを誤ると正しく位置合わせできない。chkcoo では以下の３種類の位置合わせ方法が可能。 引数３個 上記の使い方で、この場合、ターゲット天体から半径 100pix 内のカタログ中の最も明るい天体に対し、その天体の検出器上での予想位置から半径 100pix 内で最も明るい天体を同一天体であるものとみなし、観測結果の RA Dec (X,Y に変更)を平行移動して位置合わせする。拡大率は１、回転角は０で始める。位置エラーが非常に大きい場合や、同定された天体数が非常に少ない場合は、初めのペアリングが間違っているものと判断し、次に明るい星を使ってやり直す(SDSS では明るい星がカタログに含まれていないため、こういうややこしい方法になった)。初回は20magより明るい天体のみを用い、明るさが異常なものを clip しながら照合対象を徐々に暗い天体に広げて照合天体数に変化がなくなるまで繰り返す。 引数３個、２天体をマニュアル指定 TriCCS の結果(C0/tcXXXXXXX0.coo など) の行の最後に A B のマークを入れ、カタログ中の同一天体の最後にも同様に A B を入れる。A の天体で位置合わせをし、B の天体で拡大率と回転を合わせ、初期値とする。B は無くても大体うまくいくはず。 C0/tcXXXXXXX0.coo などの中身は ... 1775.900 202.500 13.797 0.002 B ... 991.932 764.243 13.849 0.002 A ... tcXXXXXXX.sdss などの中身は ... 350.423901 +27.497135 6 15.370 14.259 13.663 13.584 13.387 B ... 350.510138 +27.551812 6 16.734 15.468 13.590 15.292 13.606 A ... という感じ。 引数７個 ./chkcoo C1/tcXXXXXXX1.coo tcXXXXXXX.sdss 6 -0.00560 -0.00199 1.00430 -0.06808 > C1/tcXXXXXXX1.sdss 上記のように、dRA dDec enl rot の順で初期値を与える。デフォルトの使い方でうまくいかない場合でも、他のバンドでは収束していれば、その最終収束値を引数として与えて走らせると多分収束する。 chkcoo のエラー出力は、 dRA　　　dDec　　enl　　　rot　　　dmag　　sigm　　sigp　rejected ... 0.01589 0.00183 1.00470 -0.01040 -0.04148 0.05805 0.41 35/63/462 のような感じで、dRA と dDec は座標のずれ(deg)、enl は拡大率、rot は回転、dmag はカタログ等級−観測等級の値(天気が良くなるほど増加)で、最も注目すべきはその後の３列の数字。sigm (mag)は等級関係の較正に使われた天体の明るさが、カタログ値に対してどの程度のσでずれているか、sigp (pix)は位置関係の較正に使われた天体の位置ずれがどの程度のσで分布しているか、rejected は sigclip で何個の天体が明るさ異常で除去されたか/更に位置異常除去されたか/全体数、を示している。dRA, dDec, enl の変化量が 1e-5 以下、rot の変化量が 1e-2 以下の条件を全て満たせば収束と判断する。ループ回数の上限は100回。うまくいけば、sigma は大体 0.1 以下、rejected ２数値の合計は全体の半分以下程度のはず。最後の１行の出力は chkcoo.log に追加して記録されるので、多数の領域の比較時に参照できる。また、Limiting mag (1σ) として、Sky 部分から評価した１σ限界等級もエラー出力に表示するようにした。 chkcoo の標準出力(リダイレクトで書き出される部分)は、 X Y RA(修正後) Dec(修正後) 等級(修正後) 等級エラー Xずれ Yずれ 対応する比較カタログの中身 となっている。上記ずれ量の単位はどちらも pix で、検出器の左右の端では光学系の歪曲のためずれが大きくなる(歪曲補正はしていないが、左右200pix より外側は位置合わせ条件から外しており、その数は上記 rejected ２番目の数字に含まれている)。→ 後に、chkcoo 内で歪曲補正も行なうこととなったため、検出器端までずれは小さいはず。但し、両端200pix を位置合わせ条件から外すのは変更なし。 適用例 以下、上記カラー画像の視野に対して SDSS DR12 カタログとの照合結果を示す ●の大きさはカタログ等級、色は観測等級－カタログ等級の値で <-0.15, -0.15～-0.05, -0.05～+0.05, +0.05～+0.15, >+0.15、位置マップでの線分は正しい位置の方向と距離を表している レンズ系だが、星型の歪曲が出ていることがわかる(歪曲補正まだしていません)。同じ日に11天体観測しているので、その内 SDSS カタログが存在する 10天体分のデータをまとめて表示したのが以下。 等級のずれは、サチらない程度の中程度の明るさの天体に関してはカタログよりも明るめ、暗い天体はカタログよりも暗めにずれが出る場合が多いようだ(もちろん多数派は正しい等級を出しているもの)。r バンドの下半分で暗めに観測される天体が多い感じだが、それ以外は特に検出器上の位置での偏りはないので、基本的には全体に共通の特性のようだ(明るく観測される原因はちゃんと調べた方がいいかも...)。 sky flat の傾斜が確認できる可能性があるので、検出器位置 y と等級のずれの関係を plot したのが以下。 r は sky flat に迷光(or 夜光輝線)が混ざっていて、y が大きい領域で割りすぎになっている事が予想されたが、結果は予想外に逆で、y が大きい所の方がカタログ値よりも明るめに観測される傾向が出ている感じだ。う～ん、謎。g,i は意外と平坦になっている感じ。 ちなみに、ダイクロイック切り分け波長付近での夜光輝線の影響を見積もってみた。以下は KOOLS で半月が出ている状態での夜光スペクトル(赤)と 130nm 幅で boxcar smoothing したもの(緑)。白色 LED 光の影響がすごいということが良く分かる図だ。下右図は中央付近を拡大したもの。 546nm の強い輝線(水銀灯から出る Hg 輝線で g,r の切り分け波長がこの辺りのはず)の影響が r バンドの幅 130nm の連続光の 5% 程度の量として寄与していることが、上右図の 611nm での緑線の段差から読み取れるが、夜光輝線の影響で実際のバンド幅の変化よりも背景面輝度の変化が大きくなっている場合でも、明るい側が割りすぎになるはずなので定性的に逆だ。可能性として有り得るのは、実は暗い側に迷光が入っていて背景の傾斜が緩和されているという場合のみだが、このデータだけではよくわからない。 SDSS カタログのない残り１領域は Gaia EDR3 カタログで合わせてみた。VizieR だと G,BP,RP の順でカタログが生成されるが、TriCCS の g,b,r の順で対応させれば大体使えそうだ。等級のばらつきが大きくなるので、収束に時間がかかることに注意(実はこの天域は Dec 58°とかなり赤緯が大きく、後述する高緯度域での座標ゆがみの影響で周辺の歪曲中心がやや下方向にずれたように見えている)。 Gaia だと明るさずれが暗い天体で明るめになる。カタログの問題なのか？他の10領域も全て Gaia で合わせてみる。 r だとばらつきが大きいが、これは Gaia G バンドのバンド幅が非常に大きいせい。全体的には Gaia の方がずれが少ない。phot での測光は単なる aperture photometry なので、Gaia はこれに近く、もっと複雑なことをやっている SDSS とは相性が悪いということか？ SDSS には何らかの model から求めた等級の他に PSF 等級というものもあり、point source だけを考えたらこちらの方が良さそうなので umag, gmag, rmag, imag, zmag の代わりに upmag, gpmag, rpmag, ipmag, zpmag というのを使ってみる。 大体 Gaia と似たような感じになったので、SDSS のカタログと比較する際は PSF 等級の方を参照すべきか。SDSS の方が暗く出ている天体は、gFlags,rFlags,iFlags などを確認してみると、 0000080000000000 = SATUR_CENTER Center is close to saturated pixel(s) のフラグが立っているようだ。ただ、このフラグが立っていても正常な値を出しているものもあるので、このフラグだけで外すのはやりすぎの感じ。まあ、フラグは気にせずに sigclip 任せで除去するので大丈夫そうだ。 SDSS でそれほど明るくない天体が暗めにずれている原因を調べるため、TriCCS 等級>16mag の天体に関して Gaia,SDSS と TriCCS の等級の違い同士を比較してみた。ずれの原因が天体カラーに以依存するのかどうかを確認するため、シンボルの色は SDSS g-i カラーで分け、0.7未満、0.7〜1.0、1.0〜1.3、1.3〜1.6、1.6以上 とした。ずれがカラーに依るのであれば、ずれの原因はバンド端位置の違いだということになる。 水平左右のずれは SDSS に、上下方向は Gaia に、右上左下方向は TriCCS に問題があるものと判断できる。上下方向の Gaia のずれは filter バンド幅の違い(特に G)、右上方向は TriCCS のサチュレーション、左下方向は何らかの原因で TriCCS のカウントが大きすぎ、左にずれているのが SDSS が暗すぎる可能性があるものだが、こういう天体は天体カラーまで考慮すると Gaia とのずれも若干大きい印象だ(大きく左にずれているものは SDSS のサチュレーションで確定だが)。TriCCS と SDSS の等級ずれのカラー依存性は無さそうで、結局、ずれの原因は特定できない。天体の明るさが実際に変化していれば、Gaia との等級ずれも大きい傾向にあることは説明できるが、これだけの頻度で変光天体が存在するのかどうか... とりあえず、ここまでの SDSS カタログとのずれをチェックしてみる。ロテータが反転したのは下から３行目と４行目の間(SDSS のない天域)、最後の１行は空が明るくなってきた影響が出ている可能性が高い。ここからわかることは以下の通り。 pixel scale は 0".35 よりも 0.4% ほど大きい (歪曲補正後では 0.7% になる) 視野回転角は 0.5°程度一時的に変化することがある(天頂付近) 典型的な測光エラーは 5% 程度、位置エラーは 0".5 程度 (歪曲補正後でも同程度) g band 10領域 番号 dRA dDec enl rot dmag sigm sigp rejected 1 -0.00486 -0.00368 1.00479 -0.40952 0.15576 0.03233 0.52 6/10/74 2 -0.00357 -0.00355 1.00492 -0.40378 0.17154 0.04767 0.53 13/19/113 3 -0.00347 -0.00463 1.00547 -0.41729 -0.05008 0.02810 0.28 4/2/17 4 -0.00314 -0.00472 1.00351 -0.38226 0.30730 0.05995 0.53 3/5/33 6 -0.00508 -0.00313 1.00473 -0.35152 -0.15225 0.03644 0.59 7/9/39 7 -0.00606 -0.00199 1.00493 -0.26403 0.00505 0.04566 0.40 11/8/68 8 0.00411 -0.00183 1.00420 0.23938 0.06530 0.03675 0.55 6/8/56 10 0.00180 -0.00445 1.00640 -0.33186 0.15962 0.03220 0.33 9/6/22 11 -0.00043 -0.00461 1.00613 -0.29832 -0.85646 0.05763 0.41 7/2/20 12 -0.00093 -0.00490 1.00446 -0.27340 2.58846 0.13299 0.69 2/2/30 r band 10領域 番号 dRA dDec enl rot dmag sigm sigp rejected 1 -0.00445 -0.00363 1.00409 -0.21296 -0.17627 0.04950 0.48 3/9/74 2 -0.00315 -0.00348 1.00415 -0.19830 -0.17231 0.05844 0.51 12/18/113 3 -0.00346 -0.00453 1.00486 -0.19078 -0.29067 0.03152 0.30 5/1/17 4 -0.00232 -0.00468 1.00320 -0.16155 -0.00482 0.04663 0.50 5/5/33 6 -0.00468 -0.00306 1.00380 -0.15880 -0.48994 0.05726 0.55 6/9/39 7 -0.00560 -0.00200 1.00455 -0.06328 -0.35003 0.05142 0.37 10/9/68 8 0.00459 -0.00184 1.00394 0.44482 -0.27526 0.04516 0.49 10/8/56 10 0.00221 -0.00419 1.00510 -0.14552 -0.10123 0.02331 0.41 9/6/22 11 0.00002 -0.00456 1.00556 -0.07981 -1.17444 0.05903 0.25 4/5/20 12 -0.00098 -0.00478 1.00347 -0.05497 2.19652 0.34184 0.61 1/2/30 i band 10領域 番号 dRA dDec enl rot dmag sigm sigp rejected 1 -0.00441 -0.00378 1.00366 -0.45904 -1.02193 0.03528 0.49 4/8/74 2 -0.00314 -0.00365 1.00397 -0.44884 -0.99037 0.04942 0.48 14/14/113 3 -0.00311 -0.00460 1.00450 -0.42411 -1.10308 0.01463 0.30 4/2/17 4 -0.00235 -0.00479 1.00261 -0.42333 -0.79929 0.08764 0.51 2/5/33 6 -0.00469 -0.00323 1.00358 -0.40781 -1.32851 0.07002 0.54 4/9/39 7 -0.00562 -0.00218 1.00416 -0.31699 -1.18502 0.04261 0.37 13/8/68 8 0.00456 -0.00200 1.00342 0.18928 -1.11441 0.03710 0.46 5/8/56 10 0.00221 -0.00424 1.00494 -0.55400 -0.94500 0.04264 0.34 6/5/22 11 0.00004 -0.00435 1.00437 -0.50116 -2.01233 0.09325 0.52 2/2/20 12 -0.00099 -0.00463 1.00360 -0.52443 1.42304 0.19417 0.55 3/2/30 ●歪曲補正 補正方法 全部重ねた gri 画像で補正するのが最も手間が少なくて楽ちんだが、g,r,i は集光レンズ系が異なるので、歪曲も異なっていることが予想される。また、歪曲の大きい両側幅 200pix をカットした領域で、カタログに対して拡大位置合わせしているため、上で作成した位置残差マップは中心からの距離 r の３乗の成分以外にも、r の１乗の成分も合わせて fit する必要がある。検出器中心を (xc,yc) として、以下の関数で SDSS, Gaia それぞれの g,r,i 結果で fitting を行った。評価値としては、歪曲補正後の位置残差ベクトルの半径方向の成分の残差２乗和を最小とするようにした。 結局、r の５乗の成分も引かないとうまく引けないことがわかったので、以下の関数で fit することにした。 中心距離：r = √((x-xc)2+(y-yc)2)/1000 歪曲ｘ　：dx = (a+br2+cr4)(x-xc)/1000 歪曲ｙ　：dy = (a+br2+cr4)(y-yc)/1000 １次の係数に関しては、歪曲補正後のカタログ位置合わせ時の拡大率調整で吸収されてしまうので、ここで必要なのは３次と５次の係数のみ。SDSS, Gaia それぞれのカタログで計算した位置残差マップから求めた上記３次, ５次の係数 b,c は以下の通り。意外と５次の寄与が大きく、当初間違えて４次だけで大体うまく行ったのは３次と５次の寄与が同程度だったことが原因だった感じだ。 3次 SDSS Gaia 両方 5次 SDSS Gaia 両方 g 3.935 4.599 4.407 g 3.879 3.567 3.651 r 3.548 3.746 3.592 r 3.633 3.630 3.659 i 3.313 3.738 3.515 i 3.416 3.299 3.424 Gaia の g バンドで５次の成分がかなり３次に流れている印象だが、まあ、平均すると２つのカタログで似たような傾向が出たので合っていそうだ。この係数を使って chkcoo とそこから呼び出す fitswarp に歪曲補正を入れることにした。その際、検出器 ID が必要になったので、fixpix, shift_ave, convphot を経由して DETID を引き継いで判断することにした。上記係数を用いて歪曲補正した chkcoo で、再度 SDSS カタログと位置合わせした結果を以下に示す。 これで OK で、歪曲中心は視野中心と固定してしまっても問題なさそうだ。これにより、pixel scale は 0".35 よりも 0.7% 程度大きいということになった。 最後に fitswarp を用いて、歪曲補正と平行移動や拡大・回転を全て補正して WCS の値が画像の端から端までほぼ完璧に合うように修正した画像を生成するようにした。 ./fitswarp tc8桁番号.fits tc8桁番号r.fits 回転角 拡大率 中心x 中心y 移動x 移動y カウント倍率 adj 最後に何でもいいので引数を１つ(上記では "adj")追加すると、header の DETID 情報を参照して各カメラに定められた歪曲補正成分を追加し、(CRPIX1,CRPIX2) を３色共通の固定値(crpix.dat に書かれた値で、このファイルは chkcoo で自動生成される)にしてその分の画像シフトも追加する(CRVAL1,CRVAL2 値が３つの CMOS で毎回異なっていることが後に判明し、この２つの値も固定値で受け渡すようにした)。chkcoo でカタログ座標と観測位置を合わせた際、そのパラメータで fitswarp を自動で走らせて tc8桁番号c.fits という画像として生成するようにした(1ADU=25mag となるようカウント倍率を決定)。r の画像に g,i を合わせて合成した画像との比較は以下の通り。下左が歪曲補正後、下中が r バンドに合わせたもの、下右は下左に対し 3x3 boxcar smoothing をかけたもの。 歪曲補正をすると背景光のうねりが強調された感じになるのが気になるが、これは変形による pixel 混ぜ合わせの場所による変化がノイズ変化として影響し、平均より＋側しか見ていない画像に対しノイズのうねりが背景光のうねりのように見えているだけのようだ。上左図に対し 3x3 boxcar smoothing をかけると、ノイズレベルが大体均一になって、背景光のうねりのように見えていた構造が無くなる(上右図)。上中画像クリックで、上図の左→中→右→SDSS と gif アニメで確認できる。 これで気がついたが、右下に何やら赤い移動天体が写っている。位置合わせの際の残差マップは、明るい天体の場合はピーク位置の検出不良が考えられるが、暗い天体の位置ずれの場合はこのように実際に移動している場合もあるということがわかった。近傍の M-dwarf カタログにある天体で、Gaia EDR3 で 0".1/yr の固有運動があり、SDSS - TriCCS のずれは 5.7pix=2" なので、固有運動の大きい近傍の星ということのようだ。 これで思いつく限りのやるべきことはすべて完了。歪曲パラメータが変わっていかない限り、特に更新は必要ないのではないかと思う。 ...と思ったが、(X,Y) → (RA, Dec) の変換で Dec が大きいときには横縮小変換(RA 座標の短縮のみ考慮)ではまずいことに今更気が付いた。 (X,Y) → (RA,Dec) 変換 fits header 情報での一次変換の後、視野中心での赤経・赤緯方向に (x,y) 座標(球面 → 平面変換で必要な tan は近似で省略、x,y の単位は degree) を設定し、視野中心で (RA,Dec)=(φ,θ)、視野内のある天体で (RA,Dec)=(α,β) とすると、まずは φ=0, θ=0 の 3D 球面上での状況を考えて (cosycosx,cosysinx,siny) を (0,1,0) 方向を軸として -θ 回転させて α,β での式と成分比較 (cosycosxcosθ-sinysinθ,cosysinx,cosycosxsinθ+sinycosθ)=(cosβcosα,cosβsinα,sinβ) より、 z 成分比較 sinβ = cosycosxsinθ+sinycosθ y 成分/x 成分 tanα = cosysinx/(cosycosxcosθ-sinysinθ) z 成分×cosθ　 cosycosxsinθcosθ+sinycos2θ = sinβcosθ x 成分×sinθ　 cosycosxcosθsinθ-sinysin2θ = cosβcosαsinθ ----------------------------------------------------------- 差 siny = sinβcosθ-cosβcosαsinθ y 成分 cosysinx = cosβsinα φ方向に関しては回転対称性により単にφが増えた分αを増やせば OK なので、α-φ を上記αとして用いればよい。tan があるので、θ=90°が特異点となってしまうが、まあ、視野内に極点が来ることはほぼないと思うので場合分けは無しとする。これで (x,y) ⇔ (α,β) 変換ができるはず。 とりあえず、±0.1°範囲内での赤経・赤緯の線がどうなるか引いてみる。下左が Dec=89°、下右が Dec=89.9°の場合。 まあ、いけてそうな感じなので、この関数を間に挟んで (X,Y) ⇔ (RA,Dec) 変換してみる。これまで、RA,Dec ベースでの位置合わせを行っていたが、これにより、x,y での位置合わせにしないとまずいということが判明、chkcoo にかなりの大改修を行った。 chkcoo 内での位置合わせ手順は以下の通り。 カタログ RA,DEC → 上記手順で x,y [deg] に変換 → WCS パラメータで逆変換し X,Y [pix] にする 検出器の X,Y → 検出器中心に対し歪曲補正 → 微小拡大・微小回転 → カタログデータとペアリング ペアリング結果から拡大率と回転角を算出し再度ペアリングを繰り返す(明るさ照合とクリッピングも) 拡大・回転が収束したら、X,Y [pix] → WCS で x,y [deg] 算出 → 上記手順で RA,Dec に変換 新しい chkcoo での収束判定条件は dX,dY,rot の変化が 1e-2 以下かつ enl の変化が 1e-5 以下で、 10領域の位置合わせ結果は以下の通り。回転角の符号が反転したのは、RA と X の符号が逆なので、X,Y での回転角が RA,Dec の回転角とは逆になるため。 g band 10領域 番号 dX dY enl rot dmag sigm sigp rejected 1 48.83 -37.66 1.00811 0.41353 0.15554 0.03251 0.65 7/9/74 2 35.76 -36.48 1.00809 0.41107 0.16928 0.04587 0.49 6/19/102 3 35.54 -47.59 1.00778 0.39331 -0.05345 0.03198 0.47 4/3/15 4 31.46 -48.76 1.00803 0.38414 0.30763 0.06029 0.62 3/5/33 6 50.63 -32.42 1.00831 0.37008 -0.15234 0.03642 0.66 7/10/39 7 55.15 -20.41 1.00798 0.27808 0.00883 0.05044 0.58 10/11/68 8 -35.49 -18.50 1.00797 -0.24546 0.06422 0.03670 0.55 6/9/55 10 -18.93 -46.11 1.00770 0.27756 0.15915 0.03296 0.41 8/7/21 11 3.55 -47.67 1.00791 0.27297 -0.83179 0.06841 0.61 4/3/18 12 9.24 -50.18 1.00781 0.27123 0.14676 0.04696 0.62 5/3/26 r band 10領域 番号 dX dY enl rot dmag sigm sigp rejected 1 44.59 -37.16 1.00700 0.21727 -0.17589 0.04984 0.71 3/9/74 2 31.46 -35.66 1.00712 0.20366 -0.17311 0.05527 0.58 6/17/102 3 34.86 -46.68 1.00660 0.20994 -0.30330 0.02158 0.46 5/2/15 4 23.22 -48.18 1.00724 0.16163 -0.00700 0.04652 0.52 4/5/33 6 46.38 -31.63 1.00709 0.16711 -0.48547 0.05175 0.69 7/9/39 7 50.94 -20.52 1.00719 0.07844 -0.35009 0.05167 0.54 10/11/68 8 -39.68 -18.55 1.00710 -0.44635 -0.27524 0.04511 0.59 9/7/55 10 -22.59 -43.30 1.00703 0.11774 -0.10707 0.02701 0.23 6/10/21 11 -0.74 -46.90 1.00663 0.07061 -1.17880 0.05838 0.64 3/2/18 12 10.45 -49.37 1.00700 0.07332 -0.17325 0.01890 0.58 6/3/26 i band 10領域 番号 dX dY enl rot dmag sigm sigp rejected 1 44.28 -38.71 1.00648 0.46365 -1.02113 0.03601 0.64 4/9/74 2 31.50 -37.56 1.00681 0.45282 -0.99045 0.04208 0.65 11/10/102 3 31.39 -47.39 1.00652 0.44366 -1.10418 0.01071 0.55 4/2/15 4 23.60 -49.46 1.00650 0.42448 -0.79900 0.08706 0.60 2/5/33 6 46.62 -33.38 1.00666 0.41965 -1.32903 0.07014 0.82 4/10/39 7 51.17 -22.41 1.00672 0.33036 -1.18502 0.04265 0.49 13/10/68 8 -39.38 -20.34 1.00651 -0.19381 -1.11318 0.03649 0.53 5/8/55 10 -22.74 -43.99 1.00687 0.53785 -0.96082 0.02298 0.33 10/3/21 11 -0.90 -44.95 1.00701 0.50128 -2.00908 0.04771 0.59 6/2/18 12 10.23 -47.63 1.00674 0.50870 -1.00579 0.02570 0.50 7/3/26 領域10で、かなりの天体が reject されてしまっている感じだ。上表で rejected 欄の数字が赤くなっている２つに対して残差マップを確認してみると以下のような感じ。 そもそも視野内の星が少ないのと、SDSS または TriCCS でサチっているような明るい星が多いため、左右端 200pix を除いた残りの星で緑色でヒゲの出ていない星だけが位置合わせに使われている。他の星の位置は大体合っているので、５天体が残っていれば大丈夫そうな感じだ。さすがに３天体以下はまずいと思うので、位置合わせに使われた天体が３天体以下になったら、Worning の表示を出して一旦停止するようにしておいた。これが出た場合は、phot.cl 内の daofind 付近の部分の条件を変えて自動検出の天体数を増やすことで回避できる。上記２領域に対して、3x3 boxcar smoothing をかけてから閾値を半分にした daofind を適用する方法を試したところ、以下のようになった。 番号 dX dY enl rot dmag sigm sigp rejected 10g -18.50 -46.00 1.00824 0.32032 0.16369 0.06055 1.05 6/7/32 10r -22.85 -42.96 1.00777 0.09801 -0.11400 0.02960 0.85 9/7/32 数的には問題なくなったが、位置エラーがやや大きい値が出ている。それぞれの図で確認してみる。 やはり、暗い天体まで加えると 3x3 boxcar smoothing をかけた後での検出ということもあって、位置エラーが全体的に大きくなるようだ。大きさや方向はランダムな感じなので、特に問題は無さそう。特に星が少ない場合でない限り、現在の daofind の閾値は妥当な感じだ。 あとは、WCS の CD1_2 と CD2_1 の符号だけが未確認で、PA が０でない画像が出てきたときに対応が必要かも。 ●モザイク重ね合わせ 歪曲補正した画像であれば、モザイク観測された画像を繋ぎ合わせて１つの画像に合成できる。 WCS の位置を初期値として、１pixel レベルで最適位置を判断して重ね合わせる。 ./imcomb OBJECT OBJECT.fits [+] 第一引数は重ね合わせる画像名のリスト(tc8桁番号c.fits)。 一番上のファイルが基準になるのでモザイク中央の画像を先頭にする。 第三引数(何でも可)を加えると、最終画像と同じサイズの画像に個々の画像をシフトしてはめ込んだもの(OBJECT_[0-2].fits)を生成する。個々の画像の比較時やカラー画像生成の際に便利。 この視野内に写っている星で y 方向の等級ずれを確認してみた。やはり単一 field (５ dither) の情報だけでははっきりしたことはわからないが、やはり r バンドには縦方向に等級ずれの傾向がありそうな感じだ。 天体のカラー依存性がないかも確認してみる。B−I の値が 0.4未満、0.4〜0.8、0.8〜1.2、1.2〜1.6、1.6以上 で plot してみると、g バンドのみ若干カラー依存性があるようにも見えるが、同一天体が５回観測されており天体数としてはかなり少ないので何とも言えない。i バンドで青い点がなくなっているのは、SDSS の i バンドがサチっていて、TriCCS の i バンド測光値が２mag 程度明るく枠外下方に出ているため。 ダイクロイックでの波長ずれを解消するために、各アームにバンドパスフィルターが追加された。結果をちゃんと確認したわけではないが、r-band の等級ずれ問題は解消されたようだ。 ●解析手順まとめ 解析 directory に全 c ソースを持って来てコンパイル gcc headerinfo.c -o headerinfo -lcfitsio -Wall gcc sigclip.c -o sigclip -lm -lcfitsio -Wall gcc boxcar.c -o boxcar -lcfitsio -Wall gcc fixpix.c -o fixpix -lcfitsio -Wall gcc shift_ave.c -o shift_ave -lm -lcfitsio -Wall gcc fitswarp.c -o fitswarp -lm -lcfitsio -Wall gcc convcoo.c -o convcoo -lm -Wall gcc convphot.c -o convphot -lm -lcfitsio -Wall gcc mktrpar.c -o mktrpar -lm -Wall gcc chkcoo.c -o chkcoo -lm -Wall gcc fill99.c -o fill99 -Wall gcc imcomb.c -o imcomb -lm -lcfitsio -Wall 観測日 directory 作成、データ転送 mkdir 2021MMDD scp ... 2021MMDD/ CMOS 毎に directory を分ける mkdir 2021MMDD/C0 mkdir 2021MMDD/C1 mkdir 2021MMDD/C2 mv 2021MMDD/*0.fits 2021MMDD/C0/ mv 2021MMDD/*1.fits 2021MMDD/C1/ mv 2021MMDD/*2.fits 2021MMDD/C2/ iraf & ds9 起動(以下、cl シェル内での作業) cl !ds9 & cd 2021MMDD/C0 dark と flat 作成 #twilight と対応する dark は 0.01sec 露出であるとして !../../headerinfo ./ !../../sigclip TRCS8桁番号.fits dark0_1.fits !../../sigclip TRCS8桁番号.fits dark0_2.fits !../../sigclip TRCS8桁番号.fits dark0_3.fits !../../sigclip TRCS8桁番号.fits dark0_4.fits !../../sigclip TRCS8桁番号.fits dark0_5.fits imcomb dark0_?.fits dark0.fits combine=average reject=none !../../sigclip TRCS8桁番号.fits tw_1.fits dark0.fits !../../sigclip TRCS8桁番号.fits tw_2.fits dark0.fits !../../sigclip TRCS8桁番号.fits tw_3.fits dark0.fits cat sigclip.log !cut -d" " -f2 sigclip.log > tw !emacs tw & #平均カウントが 10000 以上のものと平均カウントが低いものを外す #Saturated pix (３列目)がほぼ０であることも確認 imcomb @tw flat combine=median reject=none displ flat 1 #以下、典型的な単一フレーム露出時間を 5sec と仮定 !../../headerinfo ./ #dark 5sec !../../sigclip TRCS8桁番号.fits dark5.fits !mv sigclip.fits mask.fits displ dark5.fits 1 zs+ #複数取得している場合 #!../../sigclip TRCS8桁番号.fits dark5_1.fits #!mv sigclip.fits mask_1.fits #!../../sigclip TRCS8桁番号.fits dark5_2.fits #!mv sigclip.fits mask_2.fits #imdel dark5 #imcomb dark5_?.fits dark5 combine=average reject=none #imdel mask #imcomb mask_?.fits mask combine=average reject=none #Bad pixel mask 確認(平均的な pix より６倍以上 dark 安定性の悪いものを表示) #Bad pixel mask は最も多く連続取得している dark で作るべき displ mask.fits 1 zs- zr- z1=6 z2=6.1 Sky flat 作成 !../../headerinfo ./ !../../headerinfo ./ | grep -v ^\# | grep -v DARK | cut -d" " -f 1 > TRCS5 !emacs TRCS5 & #5sec の object 以外は削除 (最後に空行を残さないこと) !sed -e 's/TRCS/tmp/g' TRCS5 > tmp5 !paste -d" " TRCS5 tmp5 | sed -e 's/^/..\/..\/sigclip /g' | sed -e 's/$/ dark5.fits/g' > tmp.sh !sh tmp.sh imexam @tmp5 1 #あまりにも明るい星が写っているものはリストから削除 #!emacs tmp5 & imdel sky0,skym imcomb @tmp5 sky0 combine=median reject=none imarith sky0 / flat sky0 epar imedit #(radius = 50.) Substitution radius #(search = 5.) Search radius #(buffer = 10.) Background buffer width #(width = 20.) Background width #にして ^d で終了 imedit sky0 skym #残っている天体像(pointing 位置には通常残っているはず)の中心で b 入力、q で終了 #残っている天体像が見当たらない場合は q で終了し #imcopy sky0 skym !../../boxcar skym.fits skyb.fits 100 100 displ skyb 1 zs+ #flat の check imdel skyh imarith skym / skyb skyh displ skyh 1 zs+ #flat が正しくできていれば、この画像にはノイズ以外の細かい凹凸は無いはず #筋が入っていたら flat 不良を疑うべき #flat 合成 imdel skyflat imarith skyb * flat skyflat displ skyflat 1 zs+ !rm tmp* dark 引きと flat 割り !../../headerinfo ./ #1sec 露出の天体がある場合 !ls TRCS*.fits > TRCS1 !emacs TRCS1 & #1sec の object 以外は削除 (最後に空行を残さないこと) #以下の 5sec 画像に対する処理を 1sec 画像に対しても同様に行う !sed -e 's/TRCS/ts/g' TRCS5 > ts5 imdel @ts5 imarith @TRCS5 - dark5 @ts5 imarith @ts5 / skyflat @ts5 badpix 補正 !sed -e 's/TRCS/tr/g' TRCS5 > tr5 !paste -d" " ts5 tr5 | sed -e 's/^/..\/..\/fixpix /g' | sed -e 's/$/ mask.fits skyflat.fits/g' > tmp.sh !sh tmp.sh object 重ねあわせ !sed -e 's/TRCS/tc/g' TRCS5 > tc5 !paste -d" " tr5 tc5 | sed -e 's/^/..\/..\/shift_ave /g' | sed -e 's/$/ auto/g' > tmp.sh !sh tmp.sh imexam @tc5 1 #CR 除去確認 (何が除去されたか知りたい場合) #!sed -e 's/$/[*,*,2]/g' tc5 > tc5_2 #imexam @tc5_2 1 C0,C1,C2 それぞれで、上記 4.〜8. の処理を行なう。 ３色位置合わせと測光 noao digiphot apphot cd .. ls #C0,C1,C2 があることを確認 !cp ../trpar . #trpar のコピーは初回のみで OK #以下１天体ずつ、最後の 0,1,2 を除く7桁番号で指定する !sed -e 's/NUM/7桁番号/g' ../phot.cl > tc7桁番号.cl cl < tc7桁番号.cl cl < tc7桁番号.cl cl < tc7桁番号.cl cl < tc7桁番号.cl cl < tc7桁番号.cl cat mktrpar.log #５回程度繰り返す(収束せずに振動する場合も多いが...) #検出された天体にゴミが多い場合は、phot.cl の中の daofind の引数 #(特に background noise に相当する sigma 値)の調整が必要かも #Rot/Enl: +0.003528 +1.000270 +0.000000 +1.000000 +0.004279 +0.999674 #のように最後に表示される数字が g,r,i の回転角と拡大率を表しているので #この数値が一定に収束するのが理想 #C[0-2]/tc7桁番号[0-2]r.fits が r-band に位置合わせした fits (カラー画像作成用) #C[0-2]/tc7桁番号[0-2].phot が位置合わせ前の画像で測光した結果 #C[0-2]/tc7桁番号[0-2].coo は phot の x,y,mag,err に歪曲補正後の RA,Dec を加えたもの #mktrpar.log は３色の回転・拡大履歴で、収束したかの確認がまとめてできる 他のカタログとの比較照合 !../headerinfo C0/ RA DEC http://vizier.u-strasbg.fr/viz-bin/VizieR にアクセス 下側の箱の中に上で表示されている座標を入れる(Dec の + は外す) Go をクリック (Find Catalogs ではありません) 多くヒットするもの、g,r,i の値のあるものを１つ選ぶ(通常は SDSS DR12 か DR16 かな) カタログ名左側のカタログ番号のリンクをクリック 先頭の箱に再度座標を入れる Target dimension: 18 arcmin (必ず 18 以上の数字を入れること、でないと位置合わせに失敗します) RA_ICRS,DE_ICRS,class,upmag,gpmag,rpmag,ipmag,zpmag など選択 Gaia DR2 の場合は RA_ICRS,DE_ICRS,Gmag,BPmag,RPmag (G が r, BP が g, RP が i に大体対応) 左側で max: unlimited ascii text/plain Compute は No sort にチェック Position in: は Decimal にチェック(このタイミングですぐ上の J2000 にチェックが入るので外す) 他のチェックは全て入れない 右側の Submit をクリック 表示されたら右クリックで Save Page As ... を選択、tc7桁番号.sdss などの名前で保存 #Gaia や PanSTARRS 等のカタログは等級値の欠損が多々あるので、それらを 99.99... で置換する #sdss でも低い確率で欠損があるので心配な人は sdss でも fill99 結果と比較しておくことをお勧め !../fill99 tc7桁番号.gaia !diff fill99.dat tc7桁番号.gaia | less #OK であれば、 !mv fill99.dat tc7桁番号.gaia #以下、sdss を gaia に、5,6,7 を 4,3,5 にする(RA_ICRS,DE_ICRS,Gmag,BPmag,RPmag の場合) !rm crpix.dat !../chkcoo C0/tc7桁番号0.coo tc7桁番号.sdss 5 > C0/tc7桁番号0.sdss !../chkcoo C1/tc7桁番号1.coo tc7桁番号.sdss 6 > C1/tc7桁番号1.sdss !../chkcoo C2/tc7桁番号2.coo tc7桁番号.sdss 7 > C2/tc7桁番号2.sdss cat chkcoo.log #crpix.dat は３色位置合わせのため共通の CRPIX1 CRPIX2 CRVAL1 CRVAL2 値を受け渡すファイルで、 #存在しない場合は最適なものが chkcoo で自動生成される(続く２色は生成されたファイルを参照する) #5,6,7 は列番号で、RA,Dec,class,upmag,gpmag,rpmag,ipmag,zpmag という並びであることを想定 #sigma が 1mag 超の大きい値になっていないか、rejected 天体数は全体の半数未満で適切か等に注意 #chkcoo.log でまとめて確認できる #明るすぎる星が視野内にあってうまくいかない場合は、マークA を付けて識別する(こちら参照) #全システム効率は chkcoo.log 第６列の値を dmag として、 #eff*10^(dmag/2.5)/(単一露出時間[sec])*(Inverse_gain[e/ADU])[%] #eff=61.25(g), 71.32(r), 77.30(i) で計算できる #ここでリダイレクト出力される３つのファイルの中身は以下の通り #X Y RA Dec 等級 等級エラー Xずれ Yずれ 対応する比較カタログの中身 #X Y RA Dec 等級はカタログに合わせて修正した値で、以下の *c.fits での検出器,WCS 座標と一致 #この結果は gnuplot を用いて３つの図で表示される #シンボルの大きさ：カタログ等級 #シンボルの色：TriCCS−カタログ値 赤>0.15,紫0.05〜0.15,緑-0.05〜0.05,水-0.15〜0.05,青<-0.15 #矢印：ずれ量の30倍 #WCS を画像全体で完璧に合わせた歪曲補正後の画像 #C[0-2]/tc7桁番号[0-2]c.fits #が自動生成される。これで３色カラー画像をつくっても良い。 #この画像は 1ADU=25mag なので radius=6, buffer=12, width=12 で任意の天体の測光ができる #(epar rimexam で設定) 但し、ノイズは "c" の付いていない画像で計測すること epar rimexam #(radius = 6.) Object radius #(buffer = 12.) Background buffer width #(width = 12.) Background width #(magzero= 25.) Magnitude zero point #にして ^d で終了 !sed -e 's/.fits/c.fits/g' tc7桁番号 > tc7桁番号c imexam @tc7桁番号c 1 #"," ⇔ "n" 繰り返しで上記パラメータ(phot での測光時のパラメータ)での g,r,i 測光ができる ３色カラー画像(飛ばしても可) color rgbsun C2/tc7桁番号2c.fits C1/tc7桁番号1c.fits C0/tc7桁番号0c.fits tc7桁番号.ras sw+ rz1=-1 rz2=10 gz1=-1 gz2=10 bz1=-1 bz2=10 !gimp tc桁番号.ras & #png や jpg 等、好きなフォーマットにエクスポート モザイク重ね合わせ #OBJECT は天体名で置き換え #OBJECTc[0-2] のリストには重ねたい画像のみ残して下さい。 cd C0 ls tc*c.fits > OBJECTc0 !emacs OBJECTc0 & !../../imcomb OBJECTc0 OBJECTc0.fits cd .. cd C1 ls tc*c.fits > OBJECTc1 !emacs OBJECTc1 & !../../imcomb OBJECTc1 OBJECTc1.fits cd .. cd C2 ls tc*c.fits > OBJECTc2 !emacs OBJECTc2 & !../../imcomb OBJECTc2 OBJECTc2.fits cd .. ls C?/OBJECTc?.fits > OBJECT !../imcomb OBJECT OBJECT.fits + !rm OBJECT.fits !rm OBJECT.ras color rgbsun OBJECT_2.fits OBJECT_1.fits OBJECT_0.fits OBJECT.ras sw+ rz1=-0.5 rz2=5 gz1=-0.5 gz2=5 bz1=-0.5 bz2=5 !gimp OBJECT.ras &