| KOOLS 画像解析 |
|---|
岩室 史英 (京大宇物)
●重要なお知らせ (2024.10.19 にお知らせの追加があります)
2024.10.09 |
 | スペクトルが明るいと、分光器内部の散乱による隣のファイバーからの光の漏れこみで合成像にスパイクが入ったように見える場合があります。元画像のスペクトルのピークがサチっていないか確認して下さい。 |
つづく
●観測天体の処理
- １次元化までは標準星の解析と全く同じ
観測天体スペクトル (OBJECT.fits) CR 処理後 (OBJECTc.fits)
- １次元化と Sky 差し引き
./extract OBJECTc.fits OBJECTe.fits 704 724観測天体は通常は Sky よりも暗いので、スペクトル全体で明るさを評価すると個々のファイバーに写っている天体の明るさを正しく評価できない。輝線などがある場合は、その範囲を x 座標で指定すれば、その間のカウントで明るさが評価される。とりあえずは範囲指定なしのデフォルト設定で解析し、16nm 毎に切り分けた IFU イメージを見て最も S/N が大きい波長範囲を判断(上の例では x=[704:724])して再度走らせる。
x=[704:724] で明るさ評価して切り出し後 (OBJECTe.fits)
- 縞パターン軽減 (最終的にはより完全な手法ができたのでここは飛ばしても OK です)
./mkstripe 3.5925
./klwarp sin.fits sin.fits dx_VB.fits dy_VB.fits
./klwarp cos.fits cos.fits dx_VB.fits dy_VB.fits
sin.fits, cos.fits 両方とも flat2.fits で割る
./extract OBJECTc.fits OBJECTe.fits 704 724
(sin.fits, cos.fits が存在すれば自動で読み込んで適用する)ここまでくると、非常に弱い縞パターンが気になる。人工的に特定周期の縞画像(４つの quadrant で位相共通を仮定、今回は周期3.5925pix が最適だった)を生成し、同じように画像変形と flat2 割りをして差し引いてみる。縞画像は sin と cos を生成し、extract の中で横方向 5pix メジアン差し引き後の画像(extract.fits: 上半分が縞補正前、下半分が縞補正後)との最小二乗法で強度と位相を計算する。extract を走らせた際の表示値
zsdav: 5.32013 => 4.70358
は 5pix メジアン差し引き後の画僧のσが縞処理により 5.32013 から 4.70358 に減少したことを示している。縞パターンの周期は、この値が最小にするものを選ぶ必要があるため、決定は結構面倒...(自動化できないこともないが)。ノイズ源が 120Hz に起因するものであれば、縞の周期が毎回上記の値と一致するので上記の値を定数として使える。次回の観測で縞パターンの周期に変化があるかわかると思う。
縞画像で処理後の観測天体スペクトル (OBJECTe.fits)
- 分光標準星との割り算
./calib OBJECT_list STAR_list STAR.dat 10 > OBJECT.dat観測天体の最終画像のリストと標準星の最終画像のリストを与えれば、それぞれの line 321-323 のスペクトルを用いて (line 320 までの旧バージョン画像では S/N の最も高い sky 差し引き後のスペクトルを自動で選んで) ノイズレベルを合わせた後、平均カウント(S/N x noise 値)で重みをかけて足し合わせる。その際、対応するスカイ差し引き前のスペクトルとスカイのみのスペクトルもそれぞれのノイズレベルの２乗を掛けて足し合わせから平方根を取ることで、スペクトル内の相対ノイズレベル(最小値を１として規格化)として出力する。calib を走らせた際の表示値
OBJECT.fits: Line=274 S/N=160.1 Noise=0.37 Phofac=23.4 Phot=1389 Weight=432.7
は左からファイル名、抽出 line、columm 1 値(SN 相当値で真の SN ではないことに注意)、column 2 値(ノイズファクター)、column 3 値(測光ファクター)、測光値(column 1x2x3 値)、積算時の重み(column 1÷2 値) を表している。標準星スペクトルは ESO Spectrophotometric Standards 等からスペクトルが入手できればそれを与える(STAR.dat, ESO から ftp ダウンロードできない場合はこちらを参照)。第１列は波長(Å or nm or μm で自動判定)、第２列は flux density(最終結果の単位はこの単位と同じになります) であること。与える標準星スペクトルのデータ間隔は 0.25nm 以下であることが望ましい。これより間隔が荒い場合は、解析者の責任で補完したデータを生成して与えること(間隔が荒い場合は calib 内部で直線補完されます)。スペクトルデータが無ければスペクトルタイプを与える(B0-M3)ことで恒星大気モデル(arXiv:1807.06049)を元に、主系列の Teff, logG の値で合成した B0〜M3 までの主系列テンプレート大気スペクトルから同じタイプのものを選ぶ(但し、星間吸収の影響は未補正となります)。この標準星スペクトルは、観測されたスペクトルと共に 10mn 幅でメジアンスムージングをかけ、この比を観測天体と相対ノイズそれぞれに掛け合わせることで最終スペクトルを得る(標準星スペクトルに関しては古いデータが多くて波長方向のキャリブレーション精度が悪かったりするのと、吸収線の形状が変化している可能性もあるので、細かい吸収パターンは信用できないらしい)。その際、観測天体・標準星それぞれの測光ファクターと、単一フレームでの観測天体の露出時間/標準星の露出時間の値(上の例では 10)を用いて、観測天体の Flux の単位を標準星のものに合わせる。
観測中の同時解析など標準星のスペクトルそのものが無い場合は、calib_VB.dat または calib_VR.dat (良い条件で取得できた VR のスペクトルがあれば準備できますが、まだありません)を自動選択し、とりあえずの結果を出す。OBJECT.dat の内訳
column 1: 波長 nm
column 2: 観測天体の観測スペクトル
column 3: 観測天体の観測スペクトル中のノイズレベル
column 4: 標準星の観測スペクトル
column 5: 標準星のカタログスペクトル
column 6: 標準星の観測スペクトル(10nm 範囲内のメジアンスムージング後)
column 7: 標準星のカタログスペクトル(10nm 範囲内のメジアンスムージング後)
column 8: 観測天体の最終スペクトル (Flux の単位は与えた標準星スペクトルのものと同じ)
column 9: 観測天体の最終スペクトル中の相対ノイズレベル
B0〜M3 の恒星大気モデル 観測スペクトル 標準星のスペクトル
観測スペクトル中の緑線はノイズレベル。
標準星のスペクトルは、黄/緑がカタログ値、水/青が観測値、桃/赤がその比。
-
傾斜途中の凹凸にメジアンスムージングをかけると、傾斜の影響で凹凸の両端形状が若干影響を受けるが、10nm 幅のスムージングでは 1% レベルでの事なので、通常の観測では特に問題とはならないと思われる。それよりも、補正できるはずの大気吸収がメジアンスムージングによって補正できなくなる方が気になる。試しに、異なる２つの標準星のデータを用いて、片方の結果をもう１方でキャリブレーションしたのが以下のグラフ。
標準星同士でのキャリブレーション 大気吸収部分修正後
-
傾きが急なので縦軸を対数軸にし、上述の通りメジアンスムージングをかけてキャリブレーションしたものが赤、メジアンスムージングをかけずにキャリブレーションしたものが青となっている。大気吸収のある５つの波長帯の内、590, 628 は赤のままで問題無さそうだが、687,720,761 の３つの波長帯は明らかに青の方が正しそうだ(761 は吸収が深いため airmass の違いで完全には補正できていないようだが)。この３つの波長付近では、距離 5nm 以内でメジアンスムージングなし、その外側距離 10nm まではメジアンスムージングなし→ありに線形変化するようにした(上右図緑線)。
→ その後、VPH-red で 823nm の大気吸収も修正が必要であることが判明し、これも加えて補正する波長帯は４つとなった。
最後に、ノイズレベルを正しく出力するため、最終スペクトルを幅20nm で部分的に2次関数で fit し、3σクリップで外れ値を除去、再度残りで 2次関数 fit を３回繰り返してその位置でのノイズ値として最終スペクトルからのノイズ実測値を計算。このノイズ実測値を上記の画像から評価したノイズレベルで割り、その値を再度主要波長範囲で2次関数 fit、3σクリップで外れ値を除去、再度 fit を３回繰り返して、ノイズレベルと実測値の比を 2次関数として評価、この関数でノイズレベルを規格化することで、これまで最小値で適当に規格化していたノイズレベルが正しい値となった(即ち、大局的なノイズレベルは天体スペクトルから実測、局所的な変化はカウントの２乗積算値から出した相対ノイズレベルで決まっている状態)。下左グラフの緑線はこの方法で計算された 3σノイズレベルを表している。
観測スペクトル(大気吸収修正後) 標準星のスペクトル(大気吸収修正後)
●その他気づいたこと
- ３つのファイバーバンドルの透過率のカラー特性が、標準星観測時とフラット取得時で結構異なる変化をしているようだ。Az の回転によるファイバーバンドルのねじれの変化が関係している可能性があるので、様々な Az 位置でフラットを取得することを２～３セット繰り返し、安定性をチェックしてみるべきだと思う。(Az 軸中心でのファイバーバンドル固定方法を色々試しているようですが、改善はしていません)
- bias 画像に関して
3.5925pix 周期の縞は 10枚の bias 画像平均を個別の bias 画像から差し引いたものでも確認できるが、縞の強弱が単なる唸りでは説明できそうにないので、10枚の bias の標準偏差フレームを作ることで縞が出やすい部分と出にくい部分を調べてみた。すると下図右のような画像になり、各 quadrant の左と右が繋がらないのでこれでは唸りモデルで再現できない。なぜこれほどまでに左端と右端が不連続になるのかの原因は全く不明だが、とにかく、唸りは使わずに強度変化を決めてやる必要がある。
bias 画像 bias の標準偏差画像
-
縞の成分には垂直な 3.6pix 周期のもの(下左図)と、傾斜した 7.8pix 周期のもの(下右図)の２種類があるので、この２成分(sin,cos の４枚)を10枚の bias 画像から自動で抽出するソフトを作った。
./mkbias biaslist bias.fits
出力される画像は以下の通り。
bias.fits: 入力 bias リストの２σクリッピング平均画像
mkbias_[番号].fits: ２σクリップ後の各入力画像から上記平均を引いた画像
mkbias.fits: ４つに畳んだ標準偏差画像に M パターンに沿った boxcar スムージングをかけたもの
sin1.fits, cos1.fits: 3.6pix 成分縞テンプレート
sin2.fits, cos2.fits: 7.8pix 成分縞テンプレート
mpat.fits: M パターンマスク
unity.fits: 変形後の M パターンマスクのカウント修正に使う全面１の画像
上右図の bias 標準偏差に出ている "M パターン" を数値化した "M パターンマスク" を作り、この値が同じものに沿って cosmic ray 除去のタイミングで２種類の縞成分を除くことにした(実際は M パターンに沿った幅21pix で２成分最小２乗計算をし、結果を中央１列にのみ適用)。とりあえず、バイアス画像で動作確認 (以下は変形・フラット処理後の画像)。
3.6pix 周期成分縞画像 7.8pix 周期成分縞画像
平均引き・変形・フラット割り後のバイアス 縞軽減後
3.6pix 周期成分振幅分布 7.8pix 周期成分振幅分布
-
目で見ると 7.8pix 周期のほうが目立って見えるが、実際は 3.6pix 周期のものの方が倍以上振幅が大きい。細かい縞が見づらいこともあって、処理前後の画像の違いはまあ、こんなもんかなといった感じ。天体が写っていても大丈夫か確認してみる。
上記縞処理方法で処理後の観測天体スペクトル(１次元化後)
以前の処理法とほぼ遜色ない感じ。念のため、以前の手法も重ねて適用してどれだけ追加で除けるか調べてみた画像がこれ(上半分が適用前、下半分が適用後)。数値的にも zsdav: 4.87853 => 4.87477 でσ値がほとんど変化しないので、ほとんど縞が除けない状況になっている事が確認できた。最後に、できるだけ多くの bias 画像を用いて、縞の状態が過去どのように変化しているのか調べてみる。
以下、5/5-5/16 の bias 画像の縞変化(傾斜は縦軸に対する傾き)
| 日付 | M 傾斜 | 3.6成分周期 | 3.6成分傾斜 | 7.8成分周期 | 7.8成分傾斜 |
|---|---|---|---|---|---|
| 5/5 | 2.28571 | 3.6046965 | 0.0000000 | 7.7796849 | 0.4799947 |
| 5/6 | 1.89640 | 3.5149401 | -0.0006024 | 7.4597499 | 0.4853059 |
| 5/7 | 1.82258 | 3.6040059 | 0.0041910 | 7.7578923 | 0.4756021 |
| 5/7 | 1.92941 | 3.6042770 | 0.0033505 | 7.7790700 | 0.4692092 |
| 5/8 | 2.10048 | 3.5791595 | 0.0000000 | 7.7871561 | 0.4905279 |
| 5/9 | 14.4000 | 3.2905546 | 0.0000000 | 7.7121872 | 0.5625702 |
| 5/10 | 1.54753 | 3.6028866 | -0.0008505 | 7.7823822 | 0.4559666 |
| 5/10 | 2.09852 | 3.6034753 | -0.0008497 | 7.7769599 | 0.4723620 |
| 5/11 | 2.24684 | 3.5559916 | -0.0026949 | 7.6423381 | 0.5436668 |
| 5/14 | 3.04494 | 3.5781561 | -0.0015738 | 7.6221545 | 0.5893466 |
| 5/15 | 14.3125 | 3.6045019 | -0.0029958 | 7.7773389 | 0.4740607 |
| 5/16 | 1.46642 | 3.6048799 | -0.0029878 | 7.7715149 | 0.4722291 |
最も重要なのは 3.6成分周期の変化で、同じ日にこれが変化すると始めの縞除去のやり方の try & error での縞除去しかなくなり、自動化が難しくなる。今の所、同じ日なら大丈夫そうだが、観測の始めと終わりに10枚ずつ取得しておくほうが安全そうだ。上表で黄色の行は bias 取得枚数が５枚しかなく、２σクリッピングで異常値が除けない場合があるもの。また、5/15 の M 傾斜の異常値は、bias 標準偏差画像(４つにたたんでスムージングをかけた画像)に暗い横縞成分が入っているため。
- 標準星のデータがない状態での calib
観測中の同時解析などの際、通常は標準星を撮るのは後回しになるので、そういう場合のために標準的な "効率曲線" を calib から参照してやる必要がある。今の所、VPH_blue で取得された２つ分のデータしかないので、それらを平均したものを calib_VB.dat として書き出しておき、標準星が無い場合にはそれを参照して結果を出すように calib.c を修正した (VPH_red の場合はヘッダー内の波長情報から calib_VR.dat を自動選択するが、VPH_red で取得された標準星データがまだないので準備できていない)。
２つの標準星のカウント => Flux density 変換係数
- rmcr での２成分 fit による差し引き
非常に明るい天体を観測した場合、天体の写っているファイバーと sky のみが写っているファイバーのスペクトルが大きく異なるため、単に中間的な明るさのものを平均して スケーリングして差し引きするだけでは奇麗にスペクトルを消すことができない。127本のファイバーの内、明るいものから 40本とそれ以外の 87本でそれぞれ平均を取り、この２つのスペクトルの一次結合で各ファイバーのスペクトルを合わせてやることで、明るい天体の場合でも奇麗に消すことができた(←上下端の縞除去を行わない領域と比較することで、同時に処理した２つの縞成分も消えている事が確認できる)。
明るい天体のスペクトル (STARX.fits) CR 処理後 (STARXc.fits)
- rmcr でのスペクトル外迷光の処理
上図左上端にも薄く写っているが、ファイバースリット端から分光器室内の赤色光(780-820nm)が漏れ込んでいるようだ。分散素子を VPH-blue から VPH-red にすると中央付近に移動して更に明るくなり、CR 除去や縞パターン除去に影響を及ぼすので、単純な幅15pix (3.6 周期の約４倍で 7.8 周期の約２倍)のメジアンフィルター自己画像差し引きで、できるだけ縞情報を残して強制的に平坦に矯正する事にした。
780-820nm の迷光 メジアンフィルター導入後
●標準星カタログスペクトルの修正
- 標準星のカタログスペクトルと観測スペクトルが強い吸収線付近などで一致しない場合、両者の比に不要なうねりが現れる(下図緑○部分)。その場合、calib に与える標準星のカタログスペクトルのデータを以下のように "/" でコメントアウトすることで、対応するカタログスペクトルと観測スペクトルの両方が両側のデータで内挿されて接続される。
4708.0000 0.25209E+07 0.18638E+06 16.0
4724.0000 0.24923E+07 0.18552E+06 16.0
4740.0000 0.24619E+07 0.18450E+06 16.0
/ 4756.0000 0.24318E+07 0.18348E+06 16.0
/ 4772.0000 0.24045E+07 0.18264E+06 16.0
/ 4788.0000 0.23796E+07 0.18197E+06 16.0
/ 4804.0000 0.23573E+07 0.18147E+06 16.0
/ 4820.0000 0.23074E+07 0.17881E+06 16.0
/ 4836.0000 0.21930E+07 0.17108E+06 16.0
/ 4852.0000 0.17627E+07 0.13842E+06 16.0
/ 4868.0000 0.16831E+07 0.13305E+06 16.0
/ 4884.0000 0.21207E+07 0.16873E+06 16.0
4900.0000 0.22163E+07 0.17750E+06 16.0
4916.0000 0.22080E+07 0.17799E+06 16.0
4932.0000 0.21977E+07 0.17832E+06 16.0 - 標準星のカタログスペクトルの中に、明らかに大気の吸収と思われるものが残っている場合、観測スペクトルに大気吸収が残ってしまう(下図青○部分)。この場合は、カタログスペクトルを以下のように "#" でコメントアウトすることで、対応するカタログスペクトルのみが両側のデータで内挿されて接続される。
7540.0000 0.56974E+06 0.10804E+06 16.0
7556.0000 0.56577E+06 0.10775E+06 16.0
7572.0000 0.56027E+06 0.10715E+06 16.0
# 7588.0000 0.47661E+06 0.91538E+05 16.0
# 7604.0000 0.23353E+06 0.45040E+05 16.0
# 7620.0000 0.30179E+06 0.58452E+05 16.0
# 7636.0000 0.31210E+06 0.60701E+05 16.0
# 7652.0000 0.39854E+06 0.77839E+05 16.0
# 7668.0000 0.47891E+06 0.93929E+05 16.0
# 7684.0000 0.51718E+06 0.10186E+06 16.0
7700.0000 0.52752E+06 0.10433E+06 16.0
7716.0000 0.52825E+06 0.10491E+06 16.0
7732.0000 0.52703E+06 0.10510E+06 16.0
標準星のスペクトル(補正前) 観測スペクトル(補正前)
標準星のスペクトル(補正後) 観測スペクトル(補正後)
●クロストーク
- 生データの特定の bit にのみ quadrant １と４の間にクロストークがある事が判明。
以下は、明るい標準星の生データと bit 6,5,4,3,2 の画像。bit 5 にのみクロストークが出ている。
-
以下は、bit 5 の画像を左右反転させて反転前の画像とブリンクしたもの。もちろんどちらも 0 の場所、どちらも 1 の場所はあるが、逆関係のものが多い。
-
この bit だけ単純に反転して加算してみたのが下図左。一番左の quadrant だけ見ればいいが、明らかに足し過ぎでやはりそう単純な話ではない。and を取って、bit 5 が 0 の部分にのみ左右反転した bit 5 を加算したものが下図右。足し過ぎ感は無くなったが、マシになった程度。もう少しいいアイデアは無いか...
-
結局、完全な条件は探せなかったが、一番左の quadrant のみ以下の条件で修正することにした。
1. 左右反転した bit 5 画像で、周辺を含む 3x3 領域のどこかに 1 がある
2. 左側３つの平均値 - 16 よりも小さい
3. 左側２つの外挿値 - 8 よりも小さい
4. 右側３つの平均値 - 16 よりも小さい
5. 右側２つの外挿値 - 8 よりも小さい
6. 左右２つの平均値 - 8 よりも小さい
まず、左から上記 1,2,3,6 の条件全てにあてはまるものに 32 を足し、次に右から 1,4,5,6 の条件全てにあてはまるものでまだ 32 を足していないものに 32 を足す。これを２回繰り返す(最大で64加算される)。このやり方で修正したものが下図右、適用前のものが下図左。
- 2019/10/17 に CCD 読み出しクロックの修正でこの問題が解決したため、ヘッダ情報を見てこれよりも古いデータにのみ上記を適用するように修正。但し、１分以上の露出の場合、宇宙線イベントが増えて(なぜか宇宙線イベントはクロストークが出ない)この補正が逆効果になるので、長時間露出には適用しないことにした。
●rmcr に１成分追加＋閾値２次元化
- CR 除去ルーチン(rmcr)では、明るいスペクトルのグループ(object+sky)と暗いスペクトルのグループ(sky)の２つの平均スペクトルを最小二乗条件で重ねて観測スペクトルから差し引く、ということをしているが、天体が非常に明るい場合やフラット取得時と観測時の装置状態の微小変化によるずれなどにより、残差が大きくうねっている事がある。これがあると CR 除去の際に誤って天体スペクトルを引っ掛けてしまう場合があるため、残差のうねりを除くため最後に以下の手順でもう１成分差し引くことにした。
- 横方向のσが最も大きい行を探す
- 他の行を上記σ最大行と順次掛けあわせそれぞれの行の平均を計算
- 上記の値が負のものは符号を反転させて最大値を１つだけ除き重ね合わせて平均する
- 上記残差テンプレートを用いて最小二乗条件で各行に合わせて差し引く
これよる残差画像の変化が下図。暗い天体ではほぼ変化はないが、特に明るい標準星などの場合は以下のように改善される。
- これまで、残差画像に対して定数の閾値で CR 判断を行なっていたが、夜光輝線のある場所や天体の輝線が強い場合には輝線の先端を CR と誤判定する確率が高くなる。それを避けるため、全体の閾値を上げるのもつまらないので、以下の手順で場所により閾値を変えることにした。
- CR 除去後の残差画像に差し引いたスペクトル＋残差うねりモデルを加算して元に戻す
- スペクトルが写っている部分の各行と各列のカウントの平方根の平均を別々に出す
- 各列の平均値(横方向の変化)で全体の平均が１となるよう規格化後、１以下の部分を１にする
- 各行の平均値(縦方向の変化)で全体の平均が１となるよう規格化(こちらは１以下でもそのまま)
- 上記２つの値の積に√(スペクトルカウント平均)x30を掛けたものをその場所の閾値とする
- これにより CR 除去マスクを作り直し、1. から繰り返す(３回まわす)
これにより、輝線などが縦に並んで写っている部分は上記３の情報で閾値が引き上げられ、天体が入っているスペクトルが明るい部分は上記４の情報で閾値が引き上げられることになる。但し、スペクトル両端部ではかなりカウントが下がり、標準星スペクトルでの誤検出を増やしてしまうため、１以下の部分を１にすることで波長方向への変化を抑えた。以下は非常に明るい天体(下左)と中程度の明るさ(下右)の天体の場合の閾値分布の例。
- CR 除去効果の確認が毎回必要であることがわかったので、CR 除去後の画像の確認の他、CR 除去前後の画像を引き算して目で画像チェックする(以下「解析手順まとめ」で #CR 誤検出チェック として２箇所に追加済)。以下の例では、CR 除去後の画像(下左上)で default 閾値では黄色○部分に CR 引き残しがあることがわかる。前後の差し引き画像では夜光輝線の先端などが入っていない事がわかるので、閾値を default 値(rmcr を走らせた際に Threshold = で表示される値)の半分にするため rmcr の第３引数に 0.5 を与えて走らせたものがその下の図。夜光輝線を引っ掛けること無く、CR 引き残しを消せたことがわかる。
●Aperture 制御
- Object と Sky のファイバーを自動選択して決定してしまうため、空間的にどこのファイバーが Object と Sky として利用されたかがわからなかったので、これをわかるようにしたことと同時に手動で使用するファイバーを決定できるように extract.c を修正した。使い方は、まずこれまで通り
./extract OBJECTc.fits OBJECTe.fits 761 820
とすると、OBJECTe.apt というファイルが出力される。このファイルの中身は、
# Weight (0:min 99:max) + Aperture (+: Object, -: Sky, x: Object&Sky) // Sky factor 0- 1- 1- 2- 2- 2- 1- 1- 1- 3- 6 6 5 4- 3- 0- 2- 3- 8 16 18 11 6 4- 0- 1- 4- 10 28 47+ 37+ 15 6 1- 0- 1- 4- 11 35+ 62+ 80+ 46+ 13 4- 2- 0- 1- 2- 8 30+ 76+ 91+ 83+ 33+ 10 3- 2- 0- 0- 1- 5 15 49+ 83+ 99+ 64+ 23 6 2- 1- 0- 0- 3- 7 20 43+ 77+ 62+ 34+ 9 3- 2- 0- 1- 3- 7 16 28 30+ 25 12 5 2- 0- 1- 2- 6 5 10 9 8 4- 2- 0- 1- 1- 2- 3- 3- 3- 2- 1- 0- 0- 1- 1- 1- 1- 1- 1- 0- 0- 1- 0- 1- 1- 0-のようになっており、数値は各ファイバーの重み値 (0: 一番暗いファイバー、1: 一番明るいファイバー)を100倍して 0.0001 を引いた値を整数化したもので(0-99 の整数に変換される)、右側に + が付いているものは最大S/Nスペクトルで Object+Sky(line 321) として使われたもの、- が付いているものは Sky(line 322) として使われたもの、x が付いているものは両方に使われたもの(まずあり得ないと思うが、自動判定なのでもしもの場合に対応できるようにしておいた)を表している。
一番最後の実数は夜光輝線強度を合わせるための factor。 - 手動でアパーチャを変更したい場合は、上記ファイルをエディタで開いて +/- を追加/削除、必要であれば + 部分の数値の重みを変更する。以下の例は、自動選択されたアパーチャ内を均等な重みにしたもの(Sky 部分の重みは常に１なので変更しても意味なし、天体部分の数値は 99 とする必要はないが２桁だと上書きモードで修正できるので楽)。
一番最後の夜光輝線強度を合わせるための factor が削除してあるのは、改めてその値を自動計算させるため。# Weight (0:min 99:max) + Aperture (+: Object, -: Sky, x: Object&Sky) // Sky factor 0- 1- 1- 2- 2- 2- 1- 1- 1- 3- 6 6 5 4- 3- 0- 2- 3- 8 16 18 11 6 4- 0- 1- 4- 10 28 99+ 99+ 15 6 1- 0- 1- 4- 11 99+ 99+ 99+ 99+ 13 4- 2- 0- 1- 2- 8 99+ 99+ 99+ 99+ 99+ 10 3- 2- 0- 0- 1- 5 15 99+ 99+ 99+ 99+ 23 6 2- 1- 0- 0- 3- 7 20 99+ 99+ 99+ 99+ 9 3- 2- 0- 1- 3- 7 16 28 99+ 25 12 5 2- 0- 1- 2- 6 5 10 9 8 4- 2- 0- 1- 1- 2- 3- 3- 3- 2- 1- 0- 0- 1- 1- 1- 1- 1- 1- 0- 0- 1- 0- 1- 1- 0-このファイルを extract の最後の引数として与える。
./extract OBJECTc.fits OBJECTe.fits 761 820 OBJECTe.apt
apt ファイルが与えられた場合は指定されたアパーチャでの計算も追加で行い、OBJECTe.fits の一番上３行をこのアパーチャで計算したもので置き換える(下図)。この３行の意味は下から Object+Sky, Sky, Object に対応し、257-264,265-272,273-280行に相当するもの(先頭３列の意味も同じ)となっている(但し Sky は夜光輝線強度を合わせる factor をかけた後のもの)。calib.c では、この３行が追加されている場合はこの部分だけを参照するように修正している。
- 夜光輝線が object の輝線のすぐ隣にあるような場合、夜光輝線を合わせる factor の自動計算が上手く行かない場合がある。その場合は、apt ファイル最後にある夜光輝線強度を合わせる factor の値を手動で修正し、apt ファイルを extract の引数に追加して読み込ませることで factor を指定可能。また、apt ファイルで指定した aperture での夜光輝線を合わせる factor を自動計算させたい場合は、上記の例のようにこの値を削除して読み込ませる事で自動計算が行われる。
- apt ファイルを引数で与えた場合は apt ファイルの書き出しは行われず、入力 apt ファイルがそのままの形で残る。引数で与えないと新しい apt ファイルを上書きしてしまうので、上書きしたくない場合は修正した apt ファイルの名称を変更しておく方が無難。
つづく↓
●大気分散補正
- これまでは大気分散補正を真面目に考えていなかったが、そろそろ本ロテータへの搭載が見えてきたのと、暇な GW 中でちょっと時間ができたので検討を始めることにした。観測天文学「大気と波長帯」にある大気の屈折率(Owens 1967, Appl. Optics, 6, 51)の式に、岡山の標高から算出した気圧(968.2hPa)と湿度50% での水蒸気分圧(8.53hPa)を入れ、15℃における高度45°での大気分散(arcsec)を計算すると以下のようになる。
この値(d(w) とする)と、波長別の合成イメージの重心位置変化から最適な airmass と分散方向を評価し、apt ファイルにパラメータを書き出すことにする。手順は、
- 明るさ評価をする波長範囲内(中心波長 wpsf)での合成イメージにファイバー１本分の幅でスムージングをかけたものを psf とし、apt ファイルで "+" (または "x") で選択されている部分で切り取ったものを重み付きアパーチャとする。
- 適当な大気分散ベクトルを仮定し、上記 d(w)-d(wpsf) を掛け合わせて波長 w での大気分散によるアパーチャずれ量とする
- 上記の量で波長ごとにずらしたアパーチャを波長別の合成イメージにかけ合わせて空間的に積分し、全波長分足し合わせた際に最大値を与える大気分散ベクトルを探す
この手法で高度40°程度の天体の大気分散を補正したものが以下。青が補正なし、赤が補正あり。高度が低い場合、500nm 以下では結構違いが出る感じだ...
extract では、psf が出たところで一旦 apt ファイルを書き出し、それを再度読み込んでから大気分散ベクトルを探し、最後に夜光輝線強度 factor を合わせてその結果を apt ファイルに以下のように追加している。
238.381 3.83167
1.02816数値は順に大気分散ベクトルの方向(ImR を使わない場合は望遠鏡の 280-ALT 値になる感じだ)、大気分散ベクトルの絶対値、
夜光輝線強度 factorの順。apt ファイルを extract の最後の引数として与えた場合、大気分散ベクトルの方向が与えてある場合は方向を固定して探し、絶対値も与えてある場合は探さずに指定された値を用いる。夜光輝線強度 factor まで全て書かれている場合は、夜光輝線強度 factor の計算も行わずに指定された値を用いる(大気分散ベクトルを指定せずに夜光輝線強度 factor のみ与えることはできない)。最後の引数に apt ファイル名の代わりに "NOR" (No rotator の意味)と入れると、同一 directory で kif(同一ファイル番号).fits を開き、280-(ALT-STR + ALT-END)/2 値を大気分散ベクトルの方向として用いる。ImR を用いていない場合はこれで解析するのがいいと思う(ImR 使用時でも header から計算可能なはずだが、必要な状況になった場合に準備することにする)。→ header 情報で IN/OUT を自動判断し、OUT の際には NOR 相当の動作をするよう改訂(2020/11/16)。
./extract OBJECTc.fits OBJECTe.fits 761 820 NOR
●新ファイバー用の klwarp
-
基本的には同じだが以下の点を変更・追加
- flat で光が来ないファイバーにも対応(但し、３つのバンドルグループの端が消える場合には未対応)
- 波長 fitting を３次から４次に変更
３次 ４次
- 少し収束が早くなるようにアルゴリズムを修正
新ファイバーでの dome, twilight, dome/twilight の flat_VB2 画像(表示レンジ 0.7-1.3, 比は log で ±0.3 表示)。twilight でも Sky 用 fiber は透過波長特性が異なるようだ(光の当たり方が異なるのであれば問題だが...)。dome flat は twilight よりも更に差が激しく、やはり使えない感じだ。
●新ファイバー用の mkflat
- sky 用ファイバーの flat に関して
sky ファイバーは天体用ファイバーの位置からかなり離れたところにあり、拡散板を用いてもランプによる flat ではファイバーに対する光の入射角度に違いが出て、sky 用ファイバー部分は明るさだけでなく波長特性も異なる flat となりこのままでは使うことができない。twilight flat は、光の当たり方は同一になるが大気の吸収パターンや短波長側での光量不足など、波長により S/N が不足する状態となる。これを解決するため、mkflat で通常の flat の他に twilight の画像も読み込み、sky ファイバーの光量の違いと波長特性の違いを補正するモードを設けた。./mkflat flat.fits flat2.fits twilight.fits
mkflat を３引数で用いた場合には、flat2.fits の他 mkflat.fits という確認用画像が出力される。
-
この画像は、下から
A: 一次元化した flat.fits
B: 一次元化した twilight.fits
C: A を A 全体の平均スペクトルで割ったもの
D: B を B 全体の平均スペクトルで割ったもの
E: C/D で、通常 flat で問題となる成分と大気吸収残差成分の重ねあわせ
F: E に横方向 51pix のメジアンスムージングをかけたもの
上記 F のスペクトルで各ファイバーの flat スペクトルを割り、その後はこれまで通り波長特性と明るさ特性を除いた flat2.fits を生成する。これで sky 用ファイバー部分の情報も使えるはず。
●新ファイバー用の extract
- 各ファイバーの重み付けに関して
新ファイバーでは、直前のレンズによるケラレ等でファイバー毎にかなり実効効率が異なるため、以下の問題に対処する必要が出てきた。- ある程度より暗い光しか入らないファイバーは解析に加えない方がいい
- 各ファイバーの平均カウントだけで天体がどの程度入っているかを判断することができない
各ファイバーの相対効率(平均１)：t_n (n=1～117)
各ファイバーの平均カウント：c_n
各ファイバーの補正平均カウント：f_n = c_n/t_nとし、t_n>0.4 のファイバーに対して (t_n<=0.4 のファイバーは Bad fiber として使用しない)
各ファイバーの天体光含有率：s_n = (f_n - min(f_n))/(max(f_n) - min(f_n))
(f_n 値の最小値を０、最大値を１にスケーリングしたもの)
各ファイバーの実効天体光含有率：e_n = s_n*t_n
実効天体光含有率(再スケーリング後)：(e_n - min(e_n))/(max(e_n) - min(e_n))として、これまで t_n=1 の場合の天体光含有率を各スペクトルの重み(column 1 値)として用いていたものを、上記再スケーリング後の実効天体光含有率で置き換えた(新しい extract 出力画像の column 1 値)。これにより、効率が悪いが天体光が多く入っているファイバーと、天体光の入っていない効率のいいファイバーの順位付けが明確となり、天体光含有率ではなく再スケーリング後の実効天体光含有率を各スペクトルの重みとすることで効率の悪いファイバーの寄与を小さくした。
解析に含めなかった t_n<=0.4 のファイバー位置は apt ファイル内で "/" で示した。
- Sky 差し引き時の修正 factor に関して
これまで、Sky 差し引きの際は夜光輝線のパターンが最小になるようにその都度計算していたが、Sky 部分に Object 光が混ざっていると正しい値が算出できなかったため、twilight flat での相対効率を信用して差し引くことにした。これにより、各 fiber の特性代表値を１つ増やし column 1～4 に以下の値を置くことにした。column 1: 各スペクトルの重み、すなわち実効天体光含有率
(0: object 光が無いと思われるもの、1: 最も object 光が入っていると思われるもの)
column 2: ノイズファクター (fiber 幅 7pix の重みの２乗和の√で、1.612)
column 3: 規格化ファクター (ここでは column 1 と同じ値)
column 4: ファイバーの相対効率値line 1-117: ファイバー 1-117 を一次元化したもの
line 118 : ファイバー 1-117 の均等な重みでの平均値
line 119-235: 各スペクトルの重み(一番左の column 値)の大きい順に下から並べたもの
line 236 : 空行
line 237-244: 下から重み >0.2 >0.3 >0.4 >0.5 >0.6 >0.7 >0.8 >0.9
のスペクトルを各スペクトルの重みを付けて平均化したもの。
column 1: 重みの合計、column 2: ノイズファクター、column 3: 規格化ファクター
column 4: 相対効率の重み平均値(ノイズファクターは各スペクトルのノイズファ
クター×重みの２乗和の平方根、規格化ファクターは重みの重み平均値で、この値で
各 line を割る事で重み１のスペクトルのカウントに合わせることができる値)。
line 245-252: 下から重み <0.4 <0.35 <0.3 <0.25 <0.2 <0.15 <0.1 <0.05
のスペクトルを均等な重みで平均化したもの。column 1: 積算スペクトル数、
column 2: ノイズファクター、column 3: 1 と同じ、column 4: 相対効率の平均値
line 253-260: line 237-244 に対しノイズファクター/相対効率平均 の値が小さい sky スペクトルを
line 265-272 のできるだけ上のものから選んで、相対効率の平均値の比で修正して
引き算したもの。
column 1: 引き算したスペクトルの平均値を最終的なノイズファクターで割った値
column 2: ノイズファクター、column 3: 測光ファクター(line 237-244 の
重み÷規格化ファクターで、この値をスペクトルにかけることで対応する aperture
内での全カウント合計値に変換できる)、column 4: 相対効率の重み平均値
column 1 値は S/N に大体比例する値で、最大となる結果で apt ファイルを生成する
line 261 : 空行
line 262-295: スペクトルを 16nm 毎に切り分けて IFU 上のイメージにしたもの。
一番左は明るさ評価をする波長範囲内(extract の第 3,4 引数)での積分値で、ピーク
値がスペクトルの最大値と一致するように規格化したもの(単に見やすくするため)。
t_n>0.4 のファイバーのみ、t_n で割って効率補正したものを表示している。
line 296 : 空行
line 297-299: 下から、S/N 最大となる Object+Sky, Sky, Object のスペクトル。sky は大気分散
補正後の Object+Sky スペクトルから、波長ごとのファイバー相対効率平均値を
参照して補正 factor を算出し、修正したものとなっている。line 297 の値は、apt ファイルに出力されている S/N 最大となる重みとアパーチャを元に計算しているが、大気分散補正の計算の際、アパーチャ内での重み分布がスパイク状になっているとうまく大気分散がトレースできないため、apt ファイルで与えられた重み分布に対しファイバー１本分の幅でスムージングをかけてから指定されたアパーチャを適用して計算している。そのため、事前計算で算出された最大 S/N を若干下回る S/N となる。
- 大気分散補正に関して
大気分散補正は天体光の大半が視野内に入っていることが前提で、はみ出ている場合は大抵の場合大気分散０と判定される。その場合、大気分散による傾斜がスペクトルに残ることになるが、新バンドルでは視野が狭いためこの状況に陥りやすい。そのため、データごとに大気分散量を fit するのではなく、ALT, INR の角度から予想される値を用いて、視野外に出ている割合を波長毎に補正するようにした。
-
以下は、従来の方法と上記測定値を用いて視野外に出ている割合を補正する方法の比較。同じ晩に南天にある同じ標準星を ALT〜(45°,50°,35°) で観測した結果を２種類の方法で解析し、それぞれの平均値で規格化した自己相対スペクトル(下図左&中：時間の経過とともに天候が悪化している)。ある程度は補正が働いている事が確認できるが、大気分散方向で大きく視野外に出ていると補正が困難になる。下図右は別の日に同じ天体を視野中央に入れて新方式で解析した結果。少なくとも高度角45°までは視野中央に入れればほぼ問題ないことが確認できた。
これらの確認をしやすくするため、calib.c で自己相対スペクトルを calib.dat として書き出すことにした。column の順序は以下の通り。
column 1: 波長(μm)
column 2〜N+1: object の自己相対スペクトル(N 回分)
column N+2: object スペクトルの単純平均(自己相対スペクトルの分母部分)
column N+3〜N+M+2: 標準星の自己相対スペクトル(M 回分)
column N+M+3: 標準星スペクトルの単純平均(自己相対スペクトルの分母部分)
column N+M+4: ↑ とカタログスペクトルから算出した効率(%)で、大体カウントに
比例する重みで積算されているため、真の効率の６割程度の値であることに注意
-
旧バンドルでの大気分散の大きさと方向を確認するため、2020/4〜9 のいくつかの標準星データを用いて確認してみた。大気分散の強度や方向の fit 結果にはかなりのばらつきがあるが(気圧や湿度の違いによるばらつきがどの程度効くのかは不明)、旧バンドルと新バンドルとの見かけの分散強度比が 1/1.16:1/0.81 であることも考慮すると以下のようになる(天頂角<30° の点は大半が分散０に集中)。ImR を使わない場合の NOR を付けて extract を走らせるモードでは、これまでは方向のみを決めて強度は個々のデータの fitting で決めていたが、方向・強度どちらも下図緑線で決めたものを使うことにする。即ち、旧バンドルでの NOR モードと新バンドルの通常モードはどちらも大気分散の方向と強度として緑線で示された測定値を用いることにする。また、ImR の IN/OUT を header で判断して、OUT の場合には extract が NOR モード相当の動作をするよう改訂した。逆に fitting したい場合には FIT と入れることで大気分散の方向と強度を fitting で決定する(新旧バンドル共通)。
●calib で観測天体と標準星の露出時間比を自動計算
-
これまで、観測天体と標準星の露出時間比を手動で与えていたが、露出時間の異なるデータが混在する場合にも簡単に対応させるため、EXPTIME 情報を header から読み取り、スペクトル重ね合わせの際の重みを用いて平均し、露出時間の比を自動評価するようにした。以後、第３引数を省略した状態で走らせる使い方を default とする。
./calib OBJECT_list STAR_list STAR.dat > OBJECT.dat
...
OBJ: S/N=526.723 Noise=0.16837 Photfac=29.5731 Phot=2622.67 Exptime=552.708 Alt=30.9593
...
REF: S/N=381.308 Noise=0.138322 Photfac=32.8325 Phot=1731.69 Exptime=60.0041 Alt=29.0575
上記のように、これまでの情報に加えて加重平均で求めた露出時間と高度角情報が表示される。この露出時間比と標準星のカタログ flux から観測データの flux が較正される。高度角情報は参考のために表示しているもので、観測天体と標準星が大体同じ高度角で取得されているか確認するものである。この値が大きく異なる場合は、大気吸収の補正が不完全になっている可能性を考慮する必要がある(自動補正などは一切していない)。
●VPH683 に対応
- VPH-blue と比較
-
青が VPH-blue、赤が VPH683、波長分解能の違いは２倍強程度だが、V683 の方は 650nm 以下はほとんど使えない感じ。ほぼ同じ露出時間でもかなり S/N に差があり、V683 はかなり損な印象だ。
●露出時間が短い場合のカラー安定性に注意
- 2つの標準星で比較
Feige110 (高度角50°, 露出時間60秒)と HR153 (高度角70°, 露出時間1秒) の比較。どちらも７枚ずつ取得し、それぞれの平均スペクトルで規格化して plot したものが下左図(Feige110, HR153)。高度角が高い HR153 の方がカラーが安定していないことがわかる(３つ上の項目中のグラフよりも縦軸レンジが狭いことに注意)。
-
それぞれの31枚の IFU 上の波長別合成イメージをそれぞれの波長帯での平均カウントで規格化し、一番長波長側の１枚を除いた残りの30枚を10枚ずつに分け、RGB 画像にしてみた。下左が HR153, 下右が Feige110。やはり、露出時間が短いと PSF が安定せず、それがスペクトルの color 変化を引き起こしている可能性が高い(PSF がいい瞬間が最もスペクトルが青い)。こういう場合は S/N 最大条件でアパーチャを決めるのではなく、できるだけ重み一様で大きく広げたアパーチャにすべきとは思うが、自動判定はなかなか難しそうだ...(試しに手動で重み一様かつできるだけ広げたアパーチャにしたものが上右図)。
-
ところで、高度角70°でも結構大気分散が大きいことが気になった。大気分散は天頂角の tan に比例するので(高度が低くなると tan 3乗の成分も出てくるが)、KOOLS での観測結果に合わせて天頂角の２次関数(下図緑線)で決めていた大気分散は、やはり天頂角の tan (下図青線)で決める方がいい気がしてきた。(下図左は新ファイバー、右は旧ファイバーのデータ)
-
この大気分散強度を用いて extract からやり直してみたのが以下。
-
まあ、以前のものと似たり寄ったりで良くなった部分もあれば悪くなった部分もある。一応、こっちの方が理論的にも正しいはずなので、こちらのソースで差し替え(再解析が必要なほどではないと思う)。そもそも、HR153 の方は、ALT+INR の値では、青い光が４時の方向に伸びていないといけないはずが、合成イメージでは６時位の方向になっていてちょっと大気分散の方向が変(Feige110 の方は11時の方向で、こちらは ALT+INR の値と一致している)。
●Xe ランプが合わない
- これまで、Xe ランプと He,Ne が良くないことを薄々感じていたが、ちゃんと調べてみた。以下、Hg と Ne で波長較正した klwarp で Hg_VB,Ne_VB,Xe_VB 画像を変形、extract で一次元化して line 118 (117 本の単純平均)を plot したものと、NIST&理科年表のデータの比較。Xe は全体的に NIST 波長データが短波長側にずれており、波長帯端でずれが大きくなっている感じ。Xe II では合わないので、Xe I であることは間違いなさそうだが、この文献によると、「放電気体の圧力が増加すると,原子相互の衝突などが増すために波長幅は増大し,中心波長も長波長側へずれていく」とあり、発光効率を高めるために高圧で封入されている場合に波長シフトが発生する可能性がある(輝線幅は特に広がっている印象は無いので、高圧ではないと思われるが...)。青が NIST (Xe の緑は XeII)、ピンクの短い線が理科年表のデータ。
| 385-445nm | 445-505nm | 505-565nm | 565-625nm | |
|---|---|---|---|---|
| Hg |  |  |  |  |
| Ne |  |  |  |  |
| Xe |  |  |  |  |
| 625-685nm | 685-745nm | 745-805nm | 805-865nm | |
| Hg |  |  |  |  |
| Ne |  |  |  |  |
| Xe |  |  |  |  |
-
一番長波長側の Xe の強い２本の輝線のピークが見えていないが、この２本は各段に強度が強く完全にサチっているので無視。それにしても Xe は Hg,Ne と相性が悪い。どうしてこうなるのか...
- fitting 次数の変更
これまで、通常の分散素子の分散は３次程度で歪曲の次数も高くても５次だろうという予想から最大で６次までしか確認していなかったが、大塚さんより７次で fit できるという情報をもらったので klwarp.c の大改修を行った。fitting 次数を上げると少しの揺らぎでも大きく影響が出るので、元となる輝線データの精査・追加と検出数が足りない場合に次数を５次に落とす等、様々な例外処理を加えて何とか安定して動くようになった。また、NIST も理科年表も空気中での波長で表示されるため、
同定する輝線リストは空気中での波長を参照し、読み込んだ直後に klwarp.c 内部で真空波長に変換している。最終結果も空気中の波長での表示となる。まず、NIST で使われている空気の屈折率を確認するため、default で表示される空気中での波長と、設定変更して真空中の波長で表示した場合の比(下図赤線)から NIST で使われている空気の屈折率を確認した。Owens 1967 の式で、気圧 1012.5hPa の乾燥空気で計算(下図緑線)すると大体合うことが確認できた。(現在は空気中の波長のままで計算してため、この部分は関係ない)
ここで、NIST のデータは default では空気中の波長が表示されている事に今頃気が付いた... う～ん、ここはじっくり考えないといけないところ...
-
とにかく何とか klwarp.c を改修できて、結果は以下の通り。左から VPH_blue, red, 683 の順。ΔW.L. の単位は pix。
-
手持ちのデータでもやり直してみた。直近のものから並べてみる(取得日はファイル名に記載)。日によって Xe のデータはかなり状態が変わる印象だ(露出時間が短い場合もあるので、180sec 露出で撮っていれば大丈夫そうだが)。
VPH_blue
VPH_red
-
これまでの波長較正結果と比べてどの程度変化したのかを調べてみた。ΔW.L. の単位は nm。
-
上図左より VPH_blue,red,V683 で、ピンク線は 1pix の幅、取得日の違いにより色を変えている。少し上にある Hg,Ne,Xe のスペクトル比較に用いたデータは 20210802 のもので、上左図の青線。この日の４次 fit の結果はちょっと他とは異なる振る舞いをしているが、旧ソフトでは自動では同定できない輝線があったようで、その影響で輝線がない 480nm,800nm 付近での振る舞いが変わるようだ。
全体が上っているのは、NIST のデータを真空中の数値に変換した影響(これまでのものは波長が短くなってしまっていた)で、これまでのデータの波長を 0.2nm 増やせば、波長帯中央付近に関しては大体正しくなる(→ 現在は空気中の波長のままで計算しているため、2nm の増加分は無くなっている)。また、実際のデータが短波長側で折れ曲がっているのは、fitting 次数を７次に上げたことにより参照した較正用輝線の外側では急速に値が変化するため、波長補正を最後の輝線の少し外側で打ち切っていることによるもの。
新しい klwarp を走らせた際に出力される波長較正確認用の klwarp.dat の内容は以下の通り。
１行目：中心波長、Δλ/pix、有効波長範囲の画像上での右端の位置
２行目以降：
１列目：輝線を画像上で探す際の初期値 [pix]
２列目：空気中での波長 [nm] (内部では全て真空波長で計算している)
３列目：輝線の種類
４列目：中央のファイバーで、この輝線の変形後の画像(flatcal_VB.fits など)上での位置 [pix]
５列目：中央のファイバーで、この輝線の目標位置(1025-(波長-620)/0.25)とのずれ [pix]
６列目：この輝線の輝線位置の横方向のばらつき [pix] (0.3pix 以上のものは較正には用いない)
７列目：この輝線を検出できたファイバー数 (ファイバー数の半分以下のものは較正には用いない)
８列目：この輝線の目標位置からのずれの標準偏差 [pix] (1を超える場合は同定ミスの可能性大)
1,4,5 列目の値が nan となっているものは、中央のファイバーではこの輝線のずれが大きく、クリッピングされてしまったことを表している。較正に用いなかった輝線は、先頭に # が付いてコメントアウトされた状態になる。この klwarp.dat の４列目と２列目、５列目×0.25 を plot したものが、上の pix - 波長関係図。
また、この作業の際、CAL 画像に強い宇宙線がスペクトルに近い角度で入っていると、宇宙線につられて初期のトレースで失敗する事例を発見したので、CAL 画像の重ね合わせを median ではなく minmax rejection (最大値１つを削除) での平均の計算とした。
●縞パターンの変化
- いつ頃からか、検出器に現れる縞のパターンがかなり変化してきた。以前は検出器の読み出しクロック起因の高周波ノイズと思われるほぼ垂直な細かい縞が出ていたが、それがなくなり、左右45°位の角度の縞が４つの quadrant 間の強度相関が無いような感じで現れる。今後のことも考え、縞パターンが変化しても対応できるように mkbias と rmcr の大改修を行った。
まずは、bias の連続画像平均からの差分画像を４つに畳んで、bias 取得中の縞を見やすくする(下図一番左)。水平な細かい縞と左上から右下方向への縞があることがスムージングをかけると確認できる(下図左から２番目)。スムージング前の画像を２次元フーリエ変換してパワースペクトル(の平方根なので振幅スペクトル、下図左から３番目で表示は log スケール)を調べてみる。様々な縞成分があるが、中央上下端にある水平縞成分と、右下がり縞に対応する右上方向に幾つかの周波数で強い成分が存在することがわかる。bias は30枚あるので、それぞれで計算したスペクトルを重ね合わせてみる(下図左から４番目)。すると、上記の成分だけが共通して存在し、他の成分は取得中に変化した感じだ。この情報から２つの縞テンプレートを作成するが、水平方向の縞はスペクトルに影響を与える可能性があるので、画像上で水平から±30°の角度内の縞は無視して、強度の強い順に２つの成分を取り出す。quadrant 間の強度相関が無くなったようなので、これまで mkbias で作成していた強度相関の参照ファイル (mpat.fits と unity.fits) は不要となった。
-
取り出した２つの成分が以下の通り。これを画像変形して rmcr での処理で用いる。
-
rmcr では、まず天体スペクトル成分を差し引き、宇宙線ノイズを除去(下図右)する。
-
この画像に画像変形・フラット割り後の縞テンプレート２種(sin,cos パターンあるので４枚)を 30x30 pix を単位として最小２乗条件で４成分の強度を算出、これをずらしながら繰り返して４成分の強度マップを作る。これから縞成分を合成したもの下図左だが、表示スケールを合わせると非常に薄い縞となり、大して除けていない(というか、そもそも bias 取得時の縞パターンと観測中の縞パターンが異なる...)。まあ、悪化はせず若干の改善という程度だ。縞成分を差し引き後の画像が下図右。上図右と比べても大して変化はしていない...(前後の画像を拡大してブリンクすると、テンプレートパターンの縞が軽減されていることは確認できるが)
-
まあ、やっていることに間違いはないと思うので、観測中の縞とバイアス取得時の縞パターンが同じであれば、ちゃんと除けるのではないかと思う。
●標準 bias
-
大塚さんが 1000枚の bias を撮ってくれたので、早速それで標準 bias (bias1000.fits) を作ってみた。一応、当日の bias と比較して、形状に変化無ければこちらを用いることにする。以下、これまでの 20枚 bias との比較。もっと早く取得しておくべきだったと思う...
-
３週間後に磯貝くんが 800枚の bias を撮ってくれたので、上記と同様に bias0800.fits を作り、比較してみた。まあまあ一致しているような雰囲気だが、差分を取ると違いが見られる。数値としては、
# IMAGE NPIX MEAN STDDEV MIN MAX
bias0800.fits 3358720 1.076 2.425 -0.6719 163.1
bias1000.fits 3358720 1.04 2.414 -0.4482 163.4
biasdiff.fits 3358720 0.03595 0.2045 -1.275 1.323
という感じで、0.2ADU 程度の違いはある。下左図は bias0800 (-2〜5 ADU 表示)、下右図は biasdiff = bias0800 - bias1000 (-1〜1 ADU 表示)。まあ、この程度の違いは仕方がないという感じか... 両方を枚数の重みで平均したものを bias1800 として使うことにする。
-
観測当日に取得された bias は縞テンプレートの生成に必要だが、よく比較してみると bias パターンの縦縞の位置は観測日の検出器の状態により変化するようだ(２つ上の図の左右比較)。低空間周波数成分は観測日のものを反映させる必要がある事がわかったので modbias.c を追加した。
- 観測当日bias (bias0.fits) と標準bias (bias1800.fits) の融合
./modbias bias0.fits bias1800.fits bias.fitsmodbias が行なっていること
- 観測当日bias - 標準bias を計算
- 1.の結果に 3x3 median smoothing をかける
- 2.の結果の４つの quadrant を、平均値を計算しながら左右反転重ねあわせ
- 3.の結果を 11x51 範囲で boxcar smoothing (領域端では z 歪みを無くすため幅を減らす)
- 各 quadrant の平均を再現させながら4.の結果を４つの quadrant に展開
- 5.の結果に標準bias を足す
これで、11x51よりも大きいサイズの構造は観測当日bias の情報が反映され、縦縞の位置を動かさないままスムーズにすることができた。2. の median smoothing は、1pix だけ異常値が出た場合でも、その影響を受けないようにするためのもの。これがないと 11x51 に異常値が広がることが予想される。
下図左と右はそれぞれ上記 6. と 5. の結果。
●ロテータ原点移動
- 大気分散方向の確認
TriCCS の分光モード追加に伴い、較正光源の移動方向と TriCCS の Slit 方向を合わせるため、TriCCS の取り付け角を 90°回転させそれをキャンセルする方向にロテータの原点を 90°回転させた。その結果、KOOLS のファイバー上での方位と天体の大気分散の方向も 90°回転したため、その確認を行った。
-
分散方向は 90°回転角が増える方向、分散強度は以前出した関係式の通りで、上左図の青線を用いることで問題無さそうだ。2022/06/27 以降の取得データに関して、この方向を大気分散方向として自動適用するように extract.c を修正した。extract 出力画像上での方位は、左が北で下が西(ロテータ回転前は上が北で左が西だった)となっている。
●効率計算
-
この解析ソフトは、S/N 最大ポリシーで作られているため、カウントの低い部分の重みはどんどん減らされて積算されてしまい、そのままでは正しい効率を評価することが難しい。extract でファイバーの空間方向の重みを全て１とし、apt ファイルで最も光量の多いファイバーに対し 5% 以上の光が入ったファイバーを全て同じ重みで足しあわせることでできるだけカウントを減らさないようにして評価してみた(conversion factor 1.389)。以下は大気の透過率も含む VPH-blue(４回), VPH-red(３回) での測定値。VPH-blue 結果の短波長側のうねりは、flat lamp のスペクトル特性が完全に取り切れていない影響によるもの。短波長側ではスペクトル変化が激しく、かつフォーカスも変化していくため、平均的なスペクトルで割り算するだけでは不足のようだ。天体スペクトル算出の際には、ターゲットと標準星に共通する部分なのでキャンセルされるため、気にしていなかったが、効率計算の際は結構気になる...
●Sky 有観測の場合の Object - Sky 差し引き
-
Object と Sky が同数の場合は単にリスト同士で差し引きすればよいが、Qbject でノイズが決まるような明るい天体の観測の場合は、Sky の方が Object よりも少ない場合もある。その場合は、個々の Object 画像に対して、最寄りの Sky 画像を探して Object と同数の Sky リストを作る必要がある。
./obj_sky OBJECT_list SKY_list
OBJECT_list SKY_list どちらも 天体名_ファイル番号.fits のリストで、Object のファイル番号に対し、最も近いファイル番号の Sky がペアとなるように SKY_list の中身を修正する。最も近いファイル番号の Sky が複数ある場合は、最も小さいファイル番号が Sky か Object かで番号の小さい方か大きい方かを決定する。
その後、Sky の露出時間を Object の露出時間よりも短くして S/N のバランスを取る場合に遭遇したため、これに対応するために、Object と Sky の露出時間比を元画像の EXPTIME を参照して調べ、その比を補正した Sky 画像を Object 画像１枚に対して１枚ずつ生成することにした(Object 画像の "_番号" 前に "s" を付けて生成)。
また、一番最後に noise レベルを正しく与えるために、上記 Sky 画像を Object 画像に２度足して同様に一次元化と flux 較正までした後、この結果の相対ノイズレベル情報のみを抽出し、Sky 差し引き画像を用いていられた最終結果からは絶対ノイズレベル情報のみを抽出してこの２つを融合させ、正しいノイズとして最終ファイルの第９列を差し替える replnoise.c を準備した。
./replnoise Sky差し引きでの最終結果 Sky２倍上乗せでの最終結果
これにより、Sky差し引きでの最終結果のノイズ列(第９列)が差し替えられる。
●APT ファイル自動修正
-
これまでは、point source の計測を前提にしていたため、アパーチャの形状にはあまり気を遣わずに S/N 最大となる事だけを考えてきたが、apt ファイルを修正して決められた半径での固定アパーチャでの結果が出せるように apt ファイルを自動修正するソフトを準備した。
./modapt aptfile 半径 タイプ [-]
半径の単位は arcsec で、タイプは B(Box)/G(Gauss)/C(Cone) のどれか１つ、"-" は Sky 領域無し(予め Sky が引かれている場合に用いる)となる。Box は半径内が１で外が０、Gauss は半径内が Gaussian の１σまでと同じ形状で外が０、Cone は半径内が１～0.5で直線的に変化し外が０、という重みのアパーチャ。中心は与えられた apt ファイルでの重心で決定する(各ファイバー形状は六角形としている)。指定されたアパーチャに対し、各ファイバー内で三角形グリッド上の約３万点での重み評価を行い、その平均でファイバー１本の重みを決定する。apt ファイルは、アパーチャ内に関係する部分のみ書き換える。以下は一例。
修正前：
5 9 15 20 18 9 6 3- 1- 0-
5 11 20 36+ 39+ 27 14 6 3- 1-
9 25 55+ 79+ 65+ 30+ 14 5 3- 1-
7 21 56+ 99+ 97+ 59+ 24 10 4- 2-
13 31+ 71+ 99+ 54+ 37+ 13 6 3- 1-
6 15 35+ 71+ 32+ 34+ 16 7 3- 1-
5 14 29 33+ 29 15 7 3- 1- 0-
3- 6 12 16 16 11 5 3- 1- 0-
3- 4- 6 7 7 4- 2- 1- 0- 0-
1- 2- 2- 3- 3- 3- 2- 1- 0- 0-
0- 1- 1- 1- 1- 1- 0- 0- 0- 0-
半径2" タイプB スカイ無し 条件で修正後：
5 9 15 20 18 9 6 3 1 0
5 11 20 36 0+ 27 14 6 3 1
9 6+ 80+ 99+ 93+ 27+ 14 5 3 1
7 11+ 99+ 99+ 99+ 99+ 42+ 10 4 2
13 66+ 99+ 99+ 99+ 96+ 0+ 6 3 1
6 15 48+ 99+ 99+ 73+ 1+ 7 3 1
5 14 29 2+ 29 15 7 3 1 0
3 6 12 16 16 11 5 3 1 0
3 4 6 7 7 4 2 1 0 0
1 2 2 3 3 3 2 1 0 0
0 1 1 1 1 1 0 0 0 0
赤色部分が修正されたアパーチャで、その外側の数値は修正前と同じだが何のフラグも付いていないのでただの参考数値。スカイ無し条件なので外側の "-" は全て外れている。modapt は結果を標準出力に書き出すので、入力とは別名のファイルにリダイレクトして用いる。
●大気吸収補正
-
非常に S/N の高い観測を行う場合、airmass の違いによる大気吸収の出方の違いが気になってくる。ターゲットと標準星の高度はできるだけ同じになるように観測をするのが基本だが、完全に同じにはできないので S/N の高い観測を行った場合は、ある程度でも吸収の違いを補正したくなってくる。結構明るい標準星を何度も観測する機会があったので、EL の違いによる大気吸収の変化を調べてみた。
以下、EL=73°,76.5°,85°,76°(但し翌日) のある標準星のスペクトルを平均値で規格化したもの。アパーチャは 3" と十分に広げたが、露出が短いためかフレーム毎のカラー変化が大きい場合もあるようだ。下中図は、グループ毎に平均したもの。
-
700nm 以下の部分で、airmass の違いが原因と思われるカラーの違いが出ている(flat lamp の状態の違いの可能性もあるが...)。airmass との関係は空気中のエアロゾルなどの微粒子の量によって変化すると思われるため一般的な関係を出すことはできないが、この時期(12月)の空気は微粒子は少ないことが予想されるので今回の結果による補正で平均的に過補正になることはないと思う(今後しばらくは確認を続けるが)。また、760nm の吸収も airmass に関係している感じなので、これも適当にテンプレートを作って airmass に比例した修正を加える(上右図)。補正後の状態は以下。翌日はやや過補正気味になっているかも...
-
緑のスペクトルで 700nm 付近に見られるうねりは水蒸気による吸収の可能性が高く、同一の日でも変化が速いと思われるのでこれは補正しない。また、400nm 付近がかなりうねっているのは生画像でスペクトルが傾いている事によるモアレ効果の出方の違いの可能性が高いと思うが、とりあえず原因追及は保留。上右図は高度と airmass の関係。高度20°までなら単に sec(z) で計算してもほぼ問題ない感じだ。
これでうまく行ったかと思ったが、更に数日後の同じ標準星の結果を見ると(下左図)、ほぼ天頂なのに数日前とはかなり異なった結果が得られた(分光器内で焦点面が傾いた状態だったものが改善された可能性が高いが)。ここには出していないが、ターゲットにも同様な傾向が現れており、スペクトル反りは最終的には相殺されるが、ターゲットの大気吸収補正が過補正気味になっており、大気吸収補正をしない方がいい状態となっている。吸収補正をした方がいい場合もあるので、もう少し様子を見ることにする。
ところで、ついでに VPH495 に関しても標準星が相対的にどう変化するかを見てみたら結構深刻な問題を発見した。下右図は、赤桃緑青それぞれ10枚程度の連続露出の平均で、赤桃は連続して観測が行われている(但し、赤は観測開始直後)。連続観測の最中ですら、波長帯両端は全然安定しないことがわかる。
-
この原因を調べるために、１次元化した赤桃データの平均に対する比を見ていくと、観測開始付近でかなりの長時間検出器が暴れていて、下左図のような状態になっていることがわかった。一番右の quadrant の非常に暗いカウントが入っている辺りだけが選択的にダークが増えたような状態になっている(再生像がリング状に！)。これは、辿っていくと生画像から入っているため(下下図)、検出器側の問題(ビット不良の一種？)のような感じだ...
-
詳しく見ていくと、quadrant 4 の bit 7 だけが、1088ADU (２進数で0000010001000000) に近づくと 0 であるべきところが 1 になってしまう、という不安定性がある感じであることがわかった。問題はこれがいつから発生しているかだ。クロストークの出方とも関連がある可能性もあるので、いつもサチっている Xe の画像で経年変化を見てみる。
 2022/10/20 |
 2023/04/27 |
 2023/05/23 |
|---|---|---|
 2023/09/13 |
 2023/11/05 |
 2023/12/03 |
 2023/12/21 |
 2023/12/22 |
 2023/12/25 |
-
上図は、表示レンジ 1024 - 1152 で固定した Xe 画像の quadran 4 の中央のファイバー部分の画像。2020にも bit 7 異常が若干発生していた可能性もあるが、2023/9 - 2023/11 の間頃から起こり始めているようだ。クロストークの出方とは関係ないようだ。coros で、2023/10/01 以降の画像の quadrant 4 部分に対し、1120 - 1151 ADU のピクセルに対し、該当ピクセル含む左側 9pix のメジアン値+64 に近い値(±16以内)を持つ場合は、-64 しておくという処理を入れて quadrant 4 左側から縦方向を先、横方向を後にして処理してみる。結果は以下の通り。これでうまく行きそうだ。
-
その後、左側の quadrant 1,2 でも 1088ADU (２進数で0000010001000000) に近づくと、bit 6 が０になるという bit 異常が発生していることが判明した。基本的には上記と同様に対応するが、一度誤った修正が始まると芋づる式にどんどん繋がって伸びてしまうため、完了後に逆側から確認して間違った補正は元に戻すようにした。非常に運悪く cosmic ray がこの地帯の両端に入ると、その間の誤補正が修正できないが、可能性は低いのでこれには特に対処しないこととする。
気温の低下に伴い bias のカウントも下がってきた感じで、とうとう quadrant 2 の bias のカウントが危険領域に入ってしまった。この状態では、オーバースキャン領域でも bit 異常が発生するため、オーバースキャンの値の利用の仕方によっては注意が必要だ。オーバースキャン領域(quadrant 1,2 と 3,4 の間の部分のみ)では、bit 異常が発生するカウントでかつ左右の pixel のカウントに比べて異常値のジャンプが見られた場合は即修正、ということにすることでほぼ解決できた。
-
２つの宇宙線が同一 quadrant の同じ列に入った場合に、その間が誤判定されてしまうことに対処するためにかなり複雑なアルゴリズムで対処した。あまりに複雑でわかりやすく説明できないのでここでは割愛。この問題は、bias レベルを 1200 まで引き上げることで解決したようなので、観測日が 20231001～20241118 の間のデータにのみこのアルゴリズムを適用することにした。
...と思ったが、どうやら 1088+128=1216ADU (２進数で0000010011000000)に近づくと bit 異常が発生する場合があるような感じなので、20241118 以降のデータに関してはこのレベルに近づいたら注意ということで bit error を修正することにした。coros を走らせた際、bit error が100個以上検出された場合に８個の整数が表示される。左側の４個が順に quadrant 1〜4 のオーバースキャン領域で補正した bit 異常の数で、右側の４個が画像領域での bit 異常の数を表している。安定している場合は全て０になるので、これら８個の整数は表示されない。
更に、64 と 32 の bit 異常が混在して現れる現象が確認されたので(結構な頻度で混在することが判明)、これにも対処できるように修正する。幸い、bit 異常はなぜか以下のパターンでしか現れないので、「bit 異常境界値」さえ決定しておけば(2024/1/18 より前は 1088、以降は 1216)そこから±(32～64) の範囲さえ探せば良く、周辺の状況から 32 なのか 64 なのかを判別できる。ただ、32 や 64 の段差は、状況によっては正常な画像であってもたまたま起こりうるため、bit 異常による段差なのか、輝線や連続光スペクトルの裾野や弱い輝線のピークで起こっている事なのかを判別する必要がある。
-
また、ほぼ正常な画像には bit 異常修正はしない方が無難なので、第０段階としてどのタイプの bit 異常がどの程度出ているのかを各 quadrant で解析し、優先度の高いものから２タイプだけを修正することとした。優先度は -64,+64,-32,+32 の順。また、ほぼ全ての pixel で bit 異常が出ている場合は、１タイプに統一される傾向があるため、それも考慮して各 quadrant でどのタイプの補正を行うのかを決定してから補正を始める。
- 各 quadrant にどのタイプの異常がどの程度出ているかを大まかに調査、２タイプまで選択
- bit 異常候補を全て一旦修正
- 通常のオーバースキャン領域による bias 合わせと quadrant 間の gain 修正
- 輝線検出と輝線周辺部で bit 異常修正結果を再修正（輝線周囲は判断アルゴリズムが異なる）
- その他の領域で bit 異常修正結果を再修正
- 全体で再修正結果を再々修正
1.5ヶ月にわたり、2023年11月から2024年2月までの手持ちの1000枚以上の画像全てに適用できるようにほぼ毎晩改良していったので、最終版はソースが複雑になりすぎて時間が経ったら自分でも理解不能になりそうなレベルにまで成長した。修正手順の概要は以下の通り。
という感じで、初めに可能性のあるものを全部修正し、後から様々な条件でその結果を見直す、という流れとなっている。
./coros 生画像.fits 新画像.fits +
で、coros.fits も生成される。この画像は bit 異常をどのように処理したかを表すファイルで、1 は高優先度で bit 異常処理されたもの、2 は低優先度で bit 異常処理されたもの。10 を超えているものは、一旦処理されたが、後に元の状態に戻されたもの、20 を超えたものは輝線周辺での異なるアルゴリズムで判断したもの、を示している。多分、私以外は見ても良くわからないと思う。
coros を走らせると入出力ファイル名の後に以下のような表示が出る(実際は改行なしに１行表示)。
この例では、Q1/Q3 で -32 が、Q3/Q4 で +64 が出ており、その補正数が表示されている。末尾の () 内の数字はオーバースキャン領域内での補正数で、これがあまりに多いと補正無しでは横縞が入るはずだが、これに関しては詳細には調査はしていない。
|Q1| -32: 3459/+32: 162( 2) |Q2| -32: 23/+32: 172( 0) |Q3| +64: 3870/-32: 19367( 0) |Q4| +64: 53076/---:------( 0)
この例の実際の画像が以下。-32 の補正を行うと、上下の輝線の間が繋がってしまう場所が発生するが、スペクトル抽出の際はその部分はほぼ重み０となるため、他の部分の補正が行われることを優先している。
-
2024/3/5 に読み出し回路のボード交換が行われ、bit 異常問題は解決した。この日以降に取得されたデータに対しては coros.c では bit 異常処理は行わない。
●Cal lamp 取得時の背景光
-
VPH495 を使うと、背景に原因不明の連続光が発生する。slit 以外の部分も光っているが検出器全面ではないので、slit unit 付近からの光の漏れ込み(光源はファイバーバンドル自身？)があるようだ。VPH-blue や VPH-red ではこの成分は弱いのでこれまではあまり気にしなかったが、VPH495 ではこの光が flatcal でのスペクトルと輝線位置のトレースに影響を与えるので、ソフト的に大まかに差し引くことにした。
./bkgsub cal画像
で、輝線が写っている行と写っていない行を自動識別して、上下の輝線が写っていない行の１次 fit を内挿して差し引く。これで、VPH495 でのトレースが安定した。
-
Cal で cross talk がやたら大きい事案が発生(VPH-blue と VPH495 の Hg と Ne で確認)、上記背景連続光補正の前に 下250pix, 上150pix 部分の３σ clipping 平均を別々に求め、その間は線形に補間して差し引くことで対応したが(こんな場合もあるので、線形補間では完全には対応できないが...)、輝線のピーク位置の判定に影響が出る可能性もあるので、しばらく様子を見てみる。その後、この症状はシャッターのトラブルであることが判明、露出前も露出後もシャッターが開いている場合に読み出し時の電荷転送で上下に像が流れ、シャッターのトラブルが始まったか終了した時の画像が片側だけに流れる像になっている、ということで説明できる。下250pix 部分で５σ clipping で除かれたピクセル数を調べ、＋側の clip 数 - ー側の clip 数 > 1000 の場合にこのエラーが発生しているものと判断することにした。
●地球運動によるみかけの赤方偏移変化の修正
-
同一天体をかなりの頻度で分光モニタする機会があり、毎回[OIII]輝線が赤い側にずれていく傾向があることが気になっていた。Cal lamp の気温による波長変化の可能性は考えにくいので、地球の公転＋自転成分による影響を検討したところ、ほぼ説明できることが判明し、これを補正するモードも設けることにした。地球運動の影響を計算するには、観測地点の位置情報と観測日・観測方向を合わせることで計算できるはずだが、iraf の RVSAO パッケージ 内の bcvcorr の fortran ソース bcv.f があったので、これを debian の f2c パッケージで C に変換して、KOOLS 画像のヘッダ情報から Δz が計算できるようにした。
./bcv 画像.fits
観測地点の位置情報はせいめい望遠鏡の位置情報を用い、観測日(MJD)と天体の座標は fits header から読み取って以下の出力をする。
Modified Julian date (days) : 60451.5126725 Heliocentric Julian date (days) : 2460452.0126725 Geodetic longitude (degrees, west positive) : -133.596722 Geodetic latitude (degrees) : 34.576944 Altitude above sea level (meters) : 341.0 Right ascension of star (hours) : 12.175722 Declination of star (degrees) : 39.405722 Mean equator and equinox for coordinates : 2000 Barycentric correction (km/s) (returned) : -22.8994 Heliocentric correction (km/s) (returned) : -22.8892 Geocentric correction (km/s) (returned) : -0.0760665 Redshift drift correction (Δz) : -7.638411e-05
一番最後の Δz は Barycentric correction の値から計算したもの。隣同士の pixel 情報をできるだけ混ぜないというポリシーで行くと、この情報を適用するのは klwarp の段階となる。すなわち、上記の補正を Cal lamp の輝線の波長を修正することで画像変形マップに反映させる。klwarp を cal file を参照しながら走らせて変形マップを計算する際、cal file に続いて上記 Δz を続けて引数に与えることで、予めこの情報が反映された変形マップを生成する。
./klwarp flatcal0.fits flatcal.fits klwarp_dx.fits klwarp_dy.fits cal_VB.dat -7.638411e-05
自動計算モードでは
./klwarp flatcal0.fits flatcal.fits cal_VB.dat -7.638411e-05
で走らせることになる。注意すべき点は、変形マップの中に Δz 情報が反映されてしまうため、ターゲット毎に変形マップを作り直す必要があるということだ(標準星は地球の大気吸収を補正するものなので、ターゲットと共通の Δz で解析すべき)。明るい天体の長期分光モニタの際には、これを適用することが必須だと思われる。以下に、この補正の適用前と適用後の Δz の季節変化の様子を示す。桃: KOOLS 観測 z と NED で調べた z との差、緑: bcv.c で計算した Δz の値、赤: 適当に平行移動して重ねたもの、水: bcv.c で計算された Δz で修正した変形マップを用いて解析したもの、青: 適当に平行移動したもの。 下図左が VPH495 の結果、下図右が VPH-blue の結果で、VPH-blue でも補正した方が良いことがわかる。
●分光器内のフォーカスずれの影響
-
分光器のカメラレンズは、グリズム切り替えの度にフォーカス位置を変えるが、この調整機構のネジが緩んでグラグラになっていた。これは VPH495 の flatcal の klwarp による形状判断の精度が落ちてきて判明したことだが、どういう状況だったかの備忘録も兼ねてここに掲載する。
- klwarp でかなりの Trace error が表示される
- 波長決定に使われる輝線の数が減る(klwarp.dat で複数輝線コメントアウトされる)
- 波長決定精度が悪い(klwarp.ps が ±0.2pix 以上にばらける)
検出器のフォーカスが悪くなると発生する現象は以下の通り。
上記のどれかに気がついたら flatcal_V*0.fits を確認する。ピンぼけしている可能性が高いので、これを発見したら装置担当者に報告しておく。以下の図は調整前と後の状態比較。
●解析手順まとめ
-
以下、"OBJECT" => 観測天体名で置換、"STAR" => 標準星名で置換
VPH_red 使用の場合 VB => VR, VPH-blue => VPH-red
VPH683 使用の場合 VB => V683, VPH-blue => VPH683 です。
(VPH-blue 以外は利用頻度が少ないのでバグ存在の可能性が若干上がります)
各コマンドで行われることの詳細はここまでに書かれている説明を読んで下さい。
iraf => ds9 だと WCS をちゃんと表示してくれないので、変形後の画像でポインタ位置の波長情報が欲しい場合は ds9 から open でファイルを開いて表示すると WCS 欄に波長(nm)が表示されます。
- IRAF インストール
root で
apt install astro-iraf xgterm saods9
vi /etc/bash.bashrc
alias cl='irafcl' 追加
vi /etc/iraf/login.cl
以下２行のコメントアウトを外し、以下のようにしておく
set stdimage = imt2048
set imtype = "fits"
ユーザーアカウントで解析 directory に移動し、
cp /etc/iraf/login.cl . - 解析 directory に全 c ソースとデータファイルを持って来てコンパイル。
旧ファイバーバンドル(ファイバー本数 127本)のデータ解析の場合は、先に
klwarp.c, mkflat.c, rmcr.c, extract.c, calib.c の#define NFIBER 117
を
#define NFIBER 127に変更してからコンパイルすること。
gcc coros.c -o coros -lm -lcfitsio -Wall
gcc klwarp.c -o klwarp -lm -lcfitsio -Wall
gcc mkflat.c -o mkflat -lm -lcfitsio -Wall
gcc rmcr.c -o rmcr -lm -lcfitsio -Wall
gcc extract.c -o extract -lm -lcfitsio -Wall
gcc calib.c -o calib -lm -lcfitsio -Wall
gcc mkbias.c -o mkbias -lm -lcfitsio -Wall
gcc modbias.c -o modbias -lcfitsio -Wall
gcc headerinfo.c -o headerinfo -lcfitsio -Wall
gcc obj_sky.c -o obj_sky -lcfitsio -Wall
gcc modapt.c -o modapt -lm -Wall
gcc bkgsub.c -o bkgsub -lcfitsio -lm -Wall
gcc bcv.c -o bcv -lcfitsio -lm -Wall
gcc replnoise.c -o replnoise -lm -Wall - 観測日 directory 作成、データ転送
mkdir 2024XXXX
scp ... 2024XXXX/ - iraf & ds9 起動(xgterm を使うとマウスコピーで複数行を一度に送れます。以下、cl シェル内での作業)
cl
!ds9 &
cd 2024XXXX - over scan 差し引きと bad column 補正
!ls kif*.fits | sed -e 's/^/..\/coros /g' > tmp1
!ls kif*.fits | sed -e 's/kif/KLS/g' > tmp2
!paste -d" " tmp1 tmp2 > tmp.sh
!sh tmp.sh - Sky を OBJECT とは別に撮っている場合
#Sky が OBJECT と同じ名前になっている場合は観測ログファイル修正
!grep OBJECT obs.log
!../headerinfo ./ | grep OBJECT
#座標から OBJECT と Sky の番号を確認
!emacs obs.log &
#Sky は天体名に sky をくっつけて OBJECTsky としておく - 観測ログファイルとファイル名の統合
!ls KLS*.fits > tmp1
!cat obs????????.log | sort -n -k 3,3 > tmp2
!paste -d" " tmp1 tmp2 > obs.log
!cat obs.log
#KLS で始まる file 名の番号と右側のログの file 番号が揃っていることを確認 - bias 画像処理と差し引き(縞テンプレート生成も行なっています)
!grep BIAS obs.log | cut -d" " -f1 > bias
imstat @bias
!../mkbias bias bias0.fits
# imstat の値が変なものなどある場合は
# mkbias_*.fits(平均差し引き & 2σクリップ画像) を確認する方がいい
# 高空間周波数成分を bias1800.fits と置き換えて bias.fits を生成する
!../modbias bias0.fits ../bias1800.fits bias.fits
!grep -v BIAS obs.log | cut -d" " -f1 > all
imarith @all - bias @all - flat_VB0 準備
!grep FLAT obs.log | grep VPH-blue | cut -d" " -f1 > flat_VB0
imdel flat_VB0
imstat @flat_VB0
!rm tmpmean
imstat @flat_VB0 fields=mean format- > tmpmean
imarith @flat_VB0 / @tmpmean @flat_VB0
imstat @flat_VB0
imcomb @flat_VB0 flat_VB0 combine=average reject=minmax nlow=0 nhigh=1
displ flat_VB0 1#twilight もある場合(ログに NONE と書かれている場合もあるので注意)、
!grep -i twilight obs.log | grep VPH-blue | cut -d" " -f1 > twilight_VB0
imdel twilight_VB0
imstat @twilight_VB0
#平均が 500 以下のものと最大値が60000以上のものは削除かな...
#!emacs -nw twilight_VB0
!rm tmpmean
imstat @twilight_VB0 fields=mean format- > tmpmean
imarith @twilight_VB0 / @tmpmean @twilight_VB0
imstat @twilight_VB0
imcomb @twilight_VB0 twilight_VB0 combine=average reject=minmax nlow=0 nhigh=1
displ twilight_VB0 1 - cal_VB 準備
!grep Hg obs.log | grep VPH-blue | cut -d" " -f1 > Hg_VB
imdel Hg_VB
imstat @Hg_VB
#imexam @Hg_VB 1
#n キーで順次チェック、良くないものがあればリストから削除
imcomb @Hg_VB Hg_VB combine=average reject=minmax nlow=0 nhigh=1
displ Hg_VB 1
!grep Ne obs.log | grep VPH-blue | cut -d" " -f1 > Ne_VB
imdel Ne_VB
imstat @Ne_VB
#imexam @Ne_VB 1
#n キーで順次チェック、良くないものがあればリストから削除
imcomb @Ne_VB Ne_VB combine=average reject=minmax nlow=0 nhigh=1
displ Ne_VB 1
!grep Xe obs.log | grep VPH-blue | cut -d" " -f1 > Xe_VB
imdel Xe_VB
imstat @Xe_VB
#imexam @Xe_VB 1
#n キーで順次チェック、良くないものがあればリストから削除
imcomb @Xe_VB Xe_VB combine=average reject=minmax nlow=0 nhigh=1
displ Xe_VB 1
imdel cal_VB
imcomb Hg_VB,Ne_VB,Xe_VB cal_VB combine=average reject=none
!../bkgsub cal_VB.fits
!rm tmpmean
imstat cal_VB fields=mean format- > tmpmean
imarith cal_VB / @tmpmean cal_VB
displ cal_VB 1 - flatcal_VB0 & 変形マップ準備
imdel test*,flatcal_VB0
imarith flat_VB0 * 1000 test1
imarith cal_VB * 1000 test2
imarith test1 + test2 flatcal_VB0
displ flatcal_VB0 1
#地球運動による Δz を補正する場合
#!grep OBJECT obs.log
#観測の真中辺りのファイル番号を確認
#../bcv kif[ファイル番号].fits
#一番最後に表示される Δz の値を↓の klwarp コマンドの最後の引数に追加して走らせる
#このモードを使う場合はターゲット毎にここからやり直す必要があることに注意
#もちろん、通常モードで最後まで解析してその後で上記 Δz を補正するのでも可
!../klwarp flatcal_VB0.fits flatcal_VB.fits ../cal_VB.dat
# Loop ends @dxav(x エラー平均値) の値は 0.5 以下になるはず
# Loop ends の表示が無い場合は異常終了していますので岩室まで連絡下さい
!mv klwarp_dx.fits dx_VB.fits
!mv klwarp_dy.fits dy_VB.fits
!../klwarp flatcal_VB0.fits flatcal_VB.fits dx_VB.fits dy_VB.fits ../cal_VB.dat
# 最終マップで上記と同じ dxav が出ることを確認する↑
displ flatcal_VB 1
# ↑左右端近くで輝線の並びが曲がっていたら岩室まで連絡下さい
cat klwarp.dat
# ↑一番右列の値が１を超えている場合は同定ミスをしている可能性あり
# 該当する cal_*.dat で輝線の先頭に # を付けて同定リストから外してやり直し
!gnuplot ../klwarp.plt &
# ↑fitting 残差をグラフで確認(klwarp.dat の 4,5 列目を plot しただけ)
!../klwarp flat_VB0.fits flat_VB.fits dx_VB.fits dy_VB.fits
!../mkflat flat_VB.fits flat_VB2.fits
#twilight がある場合(↑ の結果は上書きされる)
!../klwarp twilight_VB0.fits twilight_VB.fits dx_VB.fits dy_VB.fits
!../mkflat flat_VB.fits flat_VB2.fits twilight_VB.fits
displ flat_VB2 1 - 縞除去用ファイル準備
!echo "sin1.fits" > tmp1 && echo "cos1.fits" >> tmp1
!echo "sin2.fits" >> tmp1 && echo "cos2.fits" >> tmp1
!sed -e 's/^/..\/klwarp /g' tmp1 > tmp2
!sed -e 's/.fits/_VB.fits/g' tmp1 > tmp3
!sed -e 's/$/ dx_VB.fits dy_VB.fits/g' tmp3 > tmp4
!paste -d" " tmp2 tmp4 > tmp.sh
!sh tmp.sh
imdel tmp3
imarith @tmp3 / flat_VB2 @tmp3 - 標準星処理
!grep STAR obs.log
!grep STAR obs.log | grep VPH-blue | cut -d" " -f1 > STAR0
imstat @STAR0
#imexam @STAR0 1
#リストには露出時間が同一で天体が写っている画像だけ残す
!sed -e 's/^/..\/klwarp /g' STAR0 > tmp1
!sed -e 's/KLS/STAR_/g' STAR0 | sed -e 's/$/ dx_VB.fits dy_VB.fits/g' > tmp2
!paste -d" " tmp1 tmp2 > tmp.sh
!sh tmp.sh
!ls STAR_*.fits > STAR
imarith @STAR / flat_VB2 @STAR
#imexam @STAR 1
!sed -e 's/^/..\/rmcr /g' STAR > tmp1
!sed -e 's/STAR/STARc/g' STAR > tmp2
#閾値を２倍にする場合
#!sed -e 's/STAR/STARc/g' STAR | sed -e 's/$/ 2/g' > tmp2
!paste -d" " tmp1 tmp2 > tmp.sh
!sh tmp.sh
!ls STARc_*.fits > STARc
imexam @STARc 1
#CR 誤検出チェック(スペクトルを引っ掛けていたら↑閾値２倍でやり直し)
!sed -e 's/STAR/tmplist/g' STAR > tmplist
imarith @STAR - @STARc @tmplist
imexam @tmplist 1
imdel @tmplist
!sed -e 's/^/..\/extract /g' STARc > tmp1
!sed -e 's/STARc/STARe/g' STARc > tmp2
!paste -d" " tmp1 tmp2 > tmp.sh
!sh tmp.sh
!ls STARe_*.fits > STARe
imexam @STARe 1
#以下は半径 3" アパーチャにする場合の例
#!cp tmp.sh tmp3
#!sed -e 's/^/..\/modapt /g' STARe | sed -e 's/fits/apt 3 B >/g' > tmp1
#!sed -e 's/fits/apt2/g' STARe > tmp2
#!paste -d" " tmp1 tmp2 > tmp.sh
#!sh tmp.sh
#!paste -d" " tmp3 tmp2 > tmp.sh
#!sh tmp.sh - 観測天体処理
#Sky はターゲット名に "sky" がくっついているものとする
!grep OBJECTsky obs.log | grep VPH-blue | cut -d" " -f1 > OBJECTsky0
imstat @OBJECTsky0
#imexam @OBJECTsky0 1
#n キーで順次チェック、リストには Sky 画像だけ残す
!sed -e 's/^/..\/klwarp /g' OBJECTsky0 > tmp1
!sed -e 's/KLS/OBJECTsky_/g' OBJECTsky0 | sed -e 's/$/ dx_VB.fits dy_VB.fits/g' > tmp2
!paste -d" " tmp1 tmp2 > tmp.sh
!sh tmp.sh
!ls OBJECTsky_*.fits > OBJECTsky
imarith @OBJECTsky / flat_VB2 @OBJECTsky
#imexam @OBJECTsky 1
!sed -e 's/^/..\/rmcr /g' OBJECTsky > tmp1
!sed -e 's/OBJECTsky/OBJECTskyc/g' OBJECTsky > tmp2
!paste -d" " tmp1 tmp2 > tmp.sh
!sh tmp.sh
!ls OBJECTskyc_*.fits > OBJECTskyc
#imexam @OBJECTskyc 1
#Sky 画像が無い場合はここから
!grep OBJECT obs.log | grep -v OBJECTsky | grep VPH-blue | cut -d" " -f1 > OBJECT0
imstat @OBJECT0
#imexam @OBJECT0 1
#n キーで順次チェック、良くないものがあればリストから削除
!sed -e 's/^/..\/klwarp /g' OBJECT0 > tmp1
!sed -e 's/KLS/OBJECT_/g' OBJECT0 | sed -e 's/$/ dx_VB.fits dy_VB.fits/g' > tmp2
!paste -d" " tmp1 tmp2 > tmp.sh
!sh tmp.sh
!ls OBJECT_*.fits > OBJECT
imarith @OBJECT / flat_VB2 @OBJECT
#imexam @OBJECT 1
#Sky 画像があれば以下４行追加
#!../obj_sky OBJECT OBJECTskyc
#!ls OBJECTs_*.fits > OBJECTs
#imarith @OBJECT - @OBJECTs @OBJECT
#imexam @OBJECT 1
!sed -e 's/^/..\/rmcr /g' OBJECT > tmp1
!sed -e 's/OBJECT/OBJECTc/g' OBJECT > tmp2
#閾値を半分にする場合
#!sed -e 's/OBJECT/OBJECTc/g' OBJECT | sed -e 's/$/ 0.5/g' > tmp2
!paste -d" " tmp1 tmp2 > tmp.sh
!sh tmp.sh
!ls OBJECTc_*.fits > OBJECTc
imexam @OBJECTc 1
#↑CR が残っていたら閾値を半分にしてやり直し
#CR 誤検出チェック
!sed -e 's/OBJECT/tmplist/g' OBJECT > tmplist
imarith @OBJECT - @OBJECTc @tmplist
imexam @tmplist 1
imdel @tmplist
!sed -e 's/^/..\/extract /g' OBJECTc > tmp1
#512 1536 は OBJECT が明るい波長帯を xmin xmax [pixel] で指定する
#600 60 など第２引数を小さくした場合は 600±60nm と解釈される
!sed -e 's/OBJECTc/OBJECTe/g' OBJECTc | sed -e 's/$/ 512 1536/g' > tmp2
!paste -d" " tmp1 tmp2 > tmp.sh
!sh tmp.sh
!ls OBJECTe_*.fits > OBJECTe
imexam @OBJECTe 1
#OBJECT が写っていないものをリストから削除するならここで以下を入れておく
#例えば 3-5行目を除く場合は(同じものを c のリストからも削除)
#!ls OBJECTe_*.fits | sed '3,5d' > OBJECTe
#!ls OBJECTc_*.fits | sed '3,5d' > OBJECTc
#Sky frame を先に差し引いた場合は apt ファイルから sky を全削除して再度 extract を走らせる
#(Sky は取得していないがアパーチャのみ変更したい場合は標準星の最後の部分参照)
#以下は、半径2" アパーチャでスカイ無しにする場合の例
#!sed -e 's/^/..\/extract /g' OBJECTc > tmp1
#!sed -e 's/^/..\/modapt /g' OBJECTe | sed -e 's/fits/apt 2 B - >/g' > tmp2
#!sed -e 's/fits/apt2/g' OBJECTe > tmp3
#!paste -d" " tmp2 tmp3 > tmp.sh
#!sh tmp.sh
#!sed -e 's/OBJECTc/OBJECTe/g' OBJECTc > tmp2
#!paste -d" " tmp1 tmp2 tmp3 > tmp.sh
#!sh tmp.sh
#imexam @OBJECTe 1
#標準星のカタログスペクトルを STAR.dat として保存
#STAR.dat の第１列は波長(Å or nm or μm)、第２列は flux density であること
#カタログスペクトル中の星の吸収部分は / で、大気の吸収部分は # でコメントアウト
!../calib OBJECTe STARe STAR.dat > OBJECT.dat
#標準星を他の標準星で解析する場合に、出力ファイル名がカタログ
#データのファイル名と同一にならないように注意して下さい。#OBJECT.dat のフォーマットは以下の通り
#column 1: 波長 nm
#column 2: 観測天体の観測スペクトル
#column 3: 観測天体の観測スペクトル中のノイズレベル
#column 4: 標準星の観測スペクトル
#column 5: 標準星のカタログスペクトル
#column 6: 標準星の観測スペクトル(10nm 範囲内のメジアンスムージング後)
#column 7: 標準星のカタログスペクトル(10nm 範囲内のメジアンスムージング後)
#column 8: 観測天体の最終スペクトル (Flux の単位は与えた標準星スペクトルのものと同じ)
#column 9: 観測天体の最終スペクトル中の相対ノイズレベル#Sky を先に引いた場合は相対ノイズレベルの値は正しいものではなくなるので以下で追加
!sed -e 's/OBJECT/tmplist/g' OBJECT > tmplist
imarith @OBJECTc + @OBJECTs @tmplist
imarith @tmplist + @OBJECTs @tmplist
!sed -e 's/^/..\/extract /g' tmplist > tmp1
!sed -e 's/tmplist/tmpliste/g' tmplist | sed -e 's/$/ 512 1536/g' > tmp2
!ls OBJECTe_*.apt2 > tmp3
!paste -d" " tmp1 tmp2 tmp3 > tmp.sh
!sh tmp.sh
!sed -e 's/OBJECTe/tmpliste/g' OBJECTe > tmpliste
!../calib tmpliste STARe STAR.dat > tmplist.dat
!../replnoise OBJECT.dat tmplist.dat
!rm tmplist*#これにより、
#column 9: 観測天体の最終スペクトル中の相対ノイズレベル
#が正しい値(差し引き Sky 情報込み)に差し替えられる。