●お知らせ
●ソース(随時更新)
実際の解析手順は一番下にあります
●ポリシー
- スペクトルを水平かつ間隔一定、波長は縦で揃えて分散一定となるよう変形(画像変形はこの１回のみ)
- LED flat は空間周波数の高い成分のみを利用(高周波 flat)
- Sky flat で空間周波数の低い成分を補正し、吸収線から波長較正も行う(変形ベクトルマップ作成)
- アンカー画像でスペクトル位置を計測、Sky 取得時とのずれを一次変換で合わせて変形量に加算
- 天体寄与の大きいスペクトルと小さいスペクトルを結合させて天体＋背景光を差し引き CR 検知
- スペクトルは S/N 最大となるようカウントで重みをつけて積算(そのうち重み指定可にする予定)
●Bias 作成
- ３σ clipping で重ね合わせ
./sigclip biasリスト bias.fits
左：CaHK, 右：Hα カメラの bias。950〜1350 で表示。
4x4 で１単位となっている。
sigclip は与えられたファイルリストの画像を 3σ clipping で重ねるコマンド。σ値の画像を sigclip.fits として書き出すが、この解析では特に利用はしない(固定の badpix はほぼ存在しないので)
●高周波 flat 作成
●Sky 解析
- Sky 用 dark 画像作成
./sigclip sky_dark sky_dark.fits
Sky を撮る直前の 100枚(立ち上げ初期は10枚でした)の dark 画像を取得し、３σ clipping で重ねる。Sky はどちらのカメラでも gain x16 で撮るため共通 bias 画像の利用でほぼ問題ないはずだが、このデータは結構長い期間利用されるものとなるのでセットで取得される。
- Sky 画像作成
./sigclip sky sky.fits sky_dark.fits flat.fits
左：CaHK, 右：Hα カメラの Sky 画像。0〜30000 で表示。
sky のリストは、青空を観測した100枚の画像のリスト。第３引数に dark、第４引数に flat を与えた場合は、dark 引きと flat 割りを行う(もちろん規格化はしない)。
- 変形ベクトルマップ作成
./mzwarp sky.fits
左：CaHK, 右：Hα カメラの変形ベクトルマップ画像。上段が dx、下段が dy で、
表示レンジは ±250(左上)、±100(右上)、±100(左下)、±40(右下)。
mzwarp がしていることは
- スペクトル部分を切り出してから 1pixel を 2x2 に細分化し、中央 1/4、上下左右 1/8、斜め４方向 1/16 の重みでカウント再配置
- 変形ベクトルマップに従い 1. の画像を変形(初期値は変形０)
- スペクトル 61本を個々にトレース、４次関数で曲がりを fit、61本での係数変化を２次関数で fit して係数変化が連続的になるようにしてから所定の y 座標(9pix/ファイバー)にする変形ベクトルマップの y 成分を計算
- １周目は波長同定は行わないが、２周目以降はスペクトルの吸収線ピーク位置を検出し、個々の吸収線位置を所定の x 座標にする(CaHK カメラは 0.005nm/pix、Hα カメラは 0.006nm/pix)ための変形ベクトルマップの x 成分を計算、その際 pixel - 波長関係は７次関数で fit し、係数変化を２次関数で fit して連続性をもたせる
- 2.〜4. を繰り返し、収束したと判断したら最後の仕上げに全体の y 方向の微小シフト調整
- 全ての情報を取り込んだ変形ベクトルマップを mzwarp_dxdy.fits として書き出し、それを適用した画像を mzwarp.fits に出力(この画像は正常に完了したことを確認するための画像)
- 分光 flat とファイバー flat 作成
./double sky.fits SKY.fits mzwarp_dxdy.fits
左：CaHK, 右：Hα カメラの Sky 画像。0〜6000 で表示。
double がしていることは
- slit 上下部分を用いて背景レベルの平坦化(短時間露出の場合はほぼ何もしていない)
- スペクトル部分を切り出してから 1pixel を 2x2 に細分化し、変形ベクトルマップで変形
- 各ファイバーに対し、0.2/0.5/0.8/1/1/10.8/0.5/0.2 の重みを付けて一次元化
- 各ファイバーの相対的なスペクトルの違いを smoothing をかけて double_rf.fits として出力
- 分光 flat を適用してファイバー間の違いをなくす
- 波長方向に積分して各ファイバーの相対的なカウントの違いを double_y.dat として出力
- ファイバー flat を適用してから配置並べ替え、バンドルイメージ double_im.bmp として出力
- 宇宙線イベント除去(詳細はターゲット解析のところで)
double_rf.fits を分光フラット (flats.fits)、double_y.dat をファイバー flat (flat.dat) にコピーして再度同じことをすれば、バンドルイメージは完全に平坦になりそうだが、実際は他にも色々行われているため１回では平坦にならない。バンドルイメージが平坦になるまで分光 flat とファイバー flat を更新しながら３回程度繰り返す。
cp double_rf.fits flats.fits
cp double_y.dat flat.dat
./double sky.fits SKY.fits mzwarp_dxdy.fits
x3 回程度
以下、１ヶ月前の分光 flat とファイバー flat を用いた状態から出発して収束させた場合のバンドルイメージ変化の例。１回目で±1% 以内に入っているが、それ以下の量もその後少し変化することがわかる。
 
左：CaHK, 右：Hα カメラのバンドルイメージ画像。
第１ステップは最小-最大レンジ(最大値は青)で表示、以降の表示は ±1% レンジで表示。
一番最後に取得した SKY スペクトルの状況を確認しておく。
sed -e 's/#//g' double.plt > tmp && gnuplot < tmp
左：CaHK, 右：Hα カメラの Sky スペクトル。
スペクトル左上の S/N 値を確認、最終的に必要なターゲットの S/N をかなり上回っていることが必要だが、この画像の情報からターゲットの解析結果に影響する内容は波長較正と分光フラット(48pix smoothing がかかっている)だけなので、多分、S/N 300 近くあれば十分でないかと思う。
●Anchor 解析
●ターゲット解析
- 一次処理
以下、Hα 画像ではc => h、C => H と置き換え
grep OBJECT mzrc.log > OBJECT_C
./reduc OBJECT_C OBJECT_rC bias.fits flat.fits
または、
./reduc 入力fits 出力fits bias.fits flat.fits
左：CaHK, 右：Hα カメラのターゲット画像。
表示レンジは左は-200〜2000、右は -400〜4000。
reduc がしていることは、
- dark (bias) 引きと flat 割り
- 強い宇宙線のピークと周辺部を NaN に置き換え
- スペクトル上下部分の背景光部分を一番上２行に入れる(後のノイズ評価で必要)
- アンプグローを３つのテンプレートパターンを合成して補正
入力は単一 fits ファイル名またはファイルリスト。入力が単一 fits ファイルの場合は出力も単一 fits ファイル、リストで入れた場合は、出力は第２引数で与えた prefix(OBJECT_rC) の後ろに自動的に _３桁番号.fits (_001.fits から)を付けて生成される。
- 細分化と変形、宇宙線処理と一次元化
./double OBJECT_rC OBJECT_dC mzwarp_dxdy.fits anc
または、
./double 入力fits 出力fits mzwarp_dxdy.fits アンカーfits
左：CaHK, 右：Hα カメラのターゲット画像(変形と宇宙線除去後)。
表示レンジは左は-40〜400、右は -80〜800。
左：CaHK, 右：Hα カメラのノイズマップ画像。
表示レンジは左は 0〜2000、右は 0〜4000。
double がしていることは、
- 1pixel を 2x2 に細分化
- reduc 結果の一番上の２行の情報を用いてノイズ２乗画像の生成
- スペクトルの上下部分を内挿してアンプグローの引き残りを除去
- Anchor を１つに指定していない場合は最寄りの前後の Anchor 情報を合成
- 変形ベクトルマップに Anchor 情報から取り出した補正ベクトルを加えて画像変形
- 各ファイバーのスペクトルをファイバーフラットで割り、ファイバーの個性を除去
- 各ファイバーのスペクトルを重み0.2/0.5/0.8/1/1/1/0.8/0.5/0.2で一次元化し、カウントの高いグループと低いグループに分けて平均、この２つのスペクトルの混ぜ合わせで変形後の画像からスペクトルを差し引き、残りの宇宙線イベントを処理
- 7. を３回繰り返して収束させる
- 各ファイバーのスペクトルを一次元化し各スペクトルの積分値で重みを付けて一次元化、prefix_３桁番号_x.dat というファイル名で書き出し
ファイルのフォーマットは、
column #1: 空気波長
column #2: 真空波長
column #3: スペクトル(重み付きカウント平均値)
column #4: １σノイズレベル
column #5: 重みなしで加算したスペクトルから逆算したスペクトル、単位 10-16erg/s/cm2/Å で、システム効率(CaHK〜4%、Hα〜8% 弱)で割り算するとおおまかなターゲットの flux density になる
- 上の重みマップをファイバーフラットで割り、バンドル上のイメージを合成
- ノイズ２乗マップは一次元化までに用いた全ての重みの２乗をかけて同様に一次元化し prefix_x.dat というファイル名で書き出し
- 全ての情報を spec.png に統合して書き出し
入力は単一 fits ファイル名または画像名の prefix。入力が単一 fits ファイルの場合は出力も単一 fits ファイル、入力、出力とも prefix で与えた場合は、ファイル名は prefix_３桁番号.fits (_001.fits から)となる。
左：CaHK, 右：Hα カメラの結果統合画像(spec.png)。
左上に重心(x,y)[arcsec]、2σ[arcsec]、重みなし平均カウント、S/N、
右上に pixel 細分化前の高カウント pixel 数(>30000〜>60000)
- 全データの合成と相対比較
./merge OBJECT_dC
左：CaHK, 右：Hα カメラの相対スペクトル画像(spec.png)。
merge がしていることは、
- 全データの重み付き平均、ノイズは重みをかけて２乗和の平方根を重みの和で割る
結果は prefix_x.dat として保存
- 各スペクトルとノイズを 1. で割る(上図青線)
- 2. を平均１に規格化(上図緑線)
- 各スペクトルの最終ファイル prefix_３桁番号_x.dat に 2. と 3. の情報を追加
ファイルのフォーマットは、
column #5: 全データの重み付き平均に対する相対スペクトル
column #6: 対応する１σノイズレベル
column #7: 平均を１にした相対スペクトル
column #8: 対応する１σノイズレベル
column #9: 元ファイルの column #5
●解析手順まとめ
以下、"OBJECT" => 観測天体名で置換です。
- 解析 directory に全 c ソースとデータファイルを持って来てコンパイル。
gcc ave.c -o ave -lm -lcfitsio -Wall
gcc mkflat.c -o mkflat -lm -lcfitsio -Wall
gcc mzwarp.c -o mzwarp -lm -lcfitsio -Wall
gcc sigclip.c -o sigclip -lm -lcfitsio -Wall
gcc reduc.c -o reduc -lm -lcfitsio -Wall
gcc double.c -o double -lm -lcfitsio -Wall
gcc merge.c -o merge -lm -Wall
- 観測日 directory ２つを作成、データ転送
#全データ持ってくる
#mzrc* (CaHK スペクトル)を全て 2026XXXXC に入れる
mkdir 2026XXXXC
mv mzrc* 2026XXXXC
#mzrh* (Hα スペクトル)を全て 2026XXXXH に入れる
mkdir 2026XXXXH
mv mzrh* 2026XXXXH
#camnum という名前のファイルにカメラ識別番号を入れておく
#CaHK は０、Hα は１
echo 0 > 2026XXXXC/camnum
echo 1 > 2026XXXXH/camnum
#以降、2026XXXXC, 2026XXXXH 両方で行う
cd 2026XXXXC
#データログファイル確認
#ファイル名 日付 時刻 露出時間 ゲイン 検出器設定温度 検出器現在温度 カメラ内回路温度
#装置温度１ 装置温度２ 装置温度３(接触悪くて721とか出ることあり) 装置温度４ 画像平均値
#画像標準偏差 画像最小値 画像最大値 スリット位置 天体名やランプ状態など、の順序で並んでいる
less mzrc.log
- bias 画像作成
grep BIAS mzrc.log | cat -n
grep BIAS mzrc.log | cut -d" " -f1 > bias
../sigclip bias bias.fits
- 高周波 flat 作成
grep INSTFLAT mzrc.log | cat -n
#Fiber Slit を後退させて100枚、前進させて100枚、の200枚で１セット
#何セットか撮れている場合は、100枚ずつ交互に flat1, flat2 にまとめる
grep INSTFLAT mzrc.log | sed -n '001,100p' | cut -d" " -f1 > flat1
grep INSTFLAT mzrc.log | sed -n '101,200p' | cut -d" " -f1 > flat2
grep INSTFLAT mzrc.log | sed -n '201,300p' | cut -d" " -f1 >> flat1
grep INSTFLAT mzrc.log | sed -n '301,400p' | cut -d" " -f1 >> flat2
...
#上記 INSTFLAT 対応の dark は上記各セット直前の10枚
#LED 点灯状態(Cont) + シャッター閉(Shut) で取得しているので以下でリストできる
grep DARK mzrc.log | grep Cont | cat -n
grep DARK mzrc.log | grep Cont > flat_dark
#3σ clipping で平均する
../sigclip flat_dark flat_dark.fits
#２種類の flat はそれぞれ上の dark を引き、明るい中央部平均が１となるように規格化
#sigclip の dark 引きモードでは、自動的に中央部が１に規格化される
../sigclip flat1 flat1.fits flat_dark.fits
../sigclip flat2 flat2.fits flat_dark.fits
#平均してから高周波 flat (flat.fits) 生成プロセスを通す
../ave flat1.fits flat2.fits flat0.fits
../mkflat flat0
mv flat0_xy.fits flat.fits
ds9 flat.fits &
- Sky 解析
grep SKYFLAT mzrc.log | cat -n
#100枚セットで取得されているが、最大カウントが 20000 超えであればどれを選んでも大体同じ
#同一日であれば複数セットをまとめても多分問題ない(以下の sed 部分削除)
#以下は第２セット選択時の例
grep SKYFLAT mzrc.log | sed -n '101,200p' | cut -d" " -f1 > sky
grep DARK mzrc.log | grep -v Cont | grep -v FeNe | cat -n
#上記 SKYFLAT 対応の dark は上記100枚の直前に取得された 10枚
#同一日であれば複数セットをまとめても多分問題ない(以下の sed 部分削除)
grep DARK mzrc.log | grep -v Cont | grep -v FeNe | sed -n '11,20p' | cut -d" " -f1 > sky_dark
../sigclip sky_dark sky_dark.fits
#sky 100枚を合わせて dark を引き、高周波 flat で割る(flat 割りも行う場合は規格化されない)
../sigclip sky sky.fits sky_dark.fits flat.fits
ds9 sky.fits &
#変形ベクトルマップ生成
../mzwarp sky.fits
#分光 flat (flats.fits),ファイバー flat (flat.dat)生成(３回程度繰り返すとファイバー flat が均一に)
../double sky.fits SKY.fits mzwarp_dxdy.fits
sxiv -z 1000 double_im.bmp &
cp double_rf.fits flats.fits
cp double_y.dat flat.dat
../double sky.fits SKY.fits mzwarp_dxdy.fits
cp double_rf.fits flats.fits
cp double_y.dat flat.dat
../double sky.fits SKY.fits mzwarp_dxdy.fits
cp double_rf.fits flats.fits
cp double_y.dat flat.dat
../double sky.fits SKY.fits mzwarp_dxdy.fits
sed -e 's/#//g' double.plt > tmp && gnuplot < tmp
sxiv spec.png &
#SKY.fits 確認(Y=91+9n にスペクトルがあり、CaHK カメラでは X=(WL-394.6)/0.005+2425
#Hα カメラでは X=(WL-654.1)/(-0.006)+2425が波長 WL に対応、x=2425 は画像中心)
ds9 SKY.fits &
#２つ目の画像はノイズ２乗画像だが、説明は天体スペクトルのところで
- Anchor 解析
#Anchor は10枚セットで取得、名前は _３桁番号(_001〜)を付けておく
grep ANCHOR mzrc.log | cat -n
grep ANCHOR mzrc.log | sed -n '01,10p' | cut -d" " -f1 > anc_001
grep ANCHOR mzrc.log | sed -n '11,20p' | cut -d" " -f1 > anc_002
grep ANCHOR mzrc.log | sed -n '21,30p' | cut -d" " -f1 > anc_003
grep ANCHOR mzrc.log | sed -n '31,40p' | cut -d" " -f1 > anc_004
grep ANCHOR mzrc.log | sed -n '41,50p' | cut -d" " -f1 > anc_005
...
#上記に対応する dark は Anchor 直前の10枚
grep DARK mzrc.log | grep FeNe | cat -n
#Anchor と枚数が異なる場合は、どの10枚が対応しているかファイル番号で判断
grep DARK mzrc.log | grep FeNe | sed -n '01,10p' | cut -d" " -f1 > anc_001_dark
grep DARK mzrc.log | grep FeNe | sed -n '11,20p' | cut -d" " -f1 > anc_002_dark
grep DARK mzrc.log | grep FeNe | sed -n '21,30p' | cut -d" " -f1 > anc_003_dark
grep DARK mzrc.log | grep FeNe | sed -n '31,40p' | cut -d" " -f1 > anc_004_dark
grep DARK mzrc.log | grep FeNe | sed -n '41,50p' | cut -d" " -f1 > anc_005_dark
...
../sigclip anc_001_dark anc_001_dark.fits
../sigclip anc_002_dark anc_002_dark.fits
../sigclip anc_003_dark anc_003_dark.fits
../sigclip anc_004_dark anc_004_dark.fits
../sigclip anc_005_dark anc_005_dark.fits
...
../sigclip anc_001 anc_001.fits anc_001_dark.fits flat.fits
../sigclip anc_002 anc_002.fits anc_002_dark.fits flat.fits
../sigclip anc_003 anc_003.fits anc_003_dark.fits flat.fits
../sigclip anc_004 anc_004.fits anc_004_dark.fits flat.fits
../sigclip anc_005 anc_005.fits anc_005_dark.fits flat.fits
...
- ターゲット解析
grep OBJECT mzrc.log
grep OBJECT mzrc.log > OBJECTC
#1〜3 行目は使わない場合の例
#grep OBJECT mzrc.log | sed '1,3d' > OBJECTC
../reduc OBJECTC OBJECT_rC bias.fits flat.fits
rm spec_*.png
../double OBJECT_rC OBJECT_dC mzwarp_dxdy.fits anc
#plot 時に参照されるテンプレートファイルは ../double.plt.template
#yrange など適当に変更しておく
#一番最後の !../cpspec.sh を残しておけば、gif アニメが作れる(最初に rm spec_*.png しておく)
#１天体に対し全ての解析が終了したら以下で gif アニメになる
convert -delay 50 -loop 0 spec_*.png OBJECT_C.gif
rm spec_*.png
../merge OBJECT_dC
#plot 時に参照されるテンプレートファイルは ../merge.plt.template
#こちらの扱いも同様(gif アニメ作る場合は rm spec_*.png しておく)
convert -delay 50 -loop 0 spec_*.png OBJECT_Cr.gif
#OBJECT_C.dat が統合後の最終スペクトル
|
