KOOLS 画像解析
http://www.kusastro.kyoto-u.ac.jp/~iwamuro/KOOLS/index.html

岩室 史英 (京大宇物)


●重要なお知らせ

    2020.10.30
    暗い天体を 2020.10.30 より前のソースで解析した場合、ファイバー間の効率ムラによる sky の凹凸の影響で正しくアパーチャが選択できていないことがわかりました。apt ファイルでアパーチャ形状がかなり不規則になっている場合は、新ソースでの1からの再解析をお勧めします。

    2021.12.27
    波長同定に Xe の暗い輝線まで加えて7次関数で fit するのが最良の波長決定方法であることが判明しました。cal 画像に Xe まで含まれている場合は一番下の解析手順の10の klwarp から、Xe を入れてない場合は9の後半(Xe 画像を作る所)から再解析お願いします(もちろん klwarp と 3つの cal*.dat を更新してから作業して下さい)。旧ファイバーのデータ再解析の際は klwarp.c の NFIBER 値を 127 に変更してからコンパイルすることに注意して下さい。また、これまでの参照波長データが空気中での値であったため、真空波長よりも 0.2nm 波長が短くなっていたことも判明しました...

●ソース(随時更新)

  • coros.c
  • klwarp.c
  • cal_VB.dat
  • cal_VR.dat
  • cal_V683.dat
  • cal_V495.dat
  • mkflat.c
  • rmcr.c
  • extract.c
  • calib.c
  • models.dat
  • mkbias.c
  • modbias.c
  • bkgsub.c
  • klwarp.plt
  • bias1800.fits
  • headerinfo.c
  • obj_sky.c
  • modapt.c

    改変履歴:
    2019.05.16: とりあえず完成
    2019.05.18: mkstripe.c 追加、extract.c を縞パターン除去に対応するように修正
    2019.05.19: extract.c, calib.c を最終 flux がより正しくなるよう修正
    2019.05.19: 縞除去確認用画像(extract.fits) を生成するよう extract.c を修正
    2019.05.20: 大気吸収のある3つの波長帯でメジアンスムージングなしとするよう calib.c を修正
    2019.05.20: 一次元化したイメージの WCS が波長(nm)となるよう extract.c を修正
    2019.05.22: ノイズレベルを最終スペクトルの実測値に合うよう calib.c を修正
    2019.05.22: CCD の正しい読み出し向きに合うよう mkstripe.c を修正
    2019.05.27: mkbias.c 追加、rmcr.c に縞除去機能追加
    2019.05.28: rmcr.c の CR 除去閾値を修正
    2019.05.29: 標準星データがないときに calib.dat を参照するよう calib.c を修正
    2019.05.30: klwarp.c の形状トレースアルゴリズム修正、cal_VB.dat に1輝線追加
    2019.05.30: klwarp.c を VPH_red にも対応、cal_VR.dat 作成
    2019.05.31: calib.dat を calib_VB.dat に名称変更、klwarp,rmcr,extract,calib fits入出力修正
    2019.05.31: rmcr.c を2成分 fit に改良
    2019.06.05: header に CD1_1 CD2_2 追加(imexam 時の warning が消えるだけで結果は同じ)
    2019.06.06: rmcr.c に手動で Threshold を与える場合の第3引数を追加
    2019.06.08: rmcr.c で波長情報により縞テンプレートの VB/VR を自動選択するよう修正
    2019.06.12: calib.c で、標準星カタログスペクトルの分解能が 0.25nm よりも悪い場合にも対応
    2019.09.09: coros.c で x=866 の列の処理をマイナー修正
    2019.09.11: calib.c で、標準星カタログスペクトルの修正(かなり下の方参照)に対応
    2019.10.08: coros.c で 一番右と一番左の quadrant のクロストークをある程度補正
    2019.10.08: mkbias.c で M パターン傾斜のサーチ範囲拡大と収束条件修正
    2019.10.11: rmcr.c に1成分追加、閾値を各行・各列の平方根カウントの平均値の積で決定
    2019.10.13: extract.c で apt ファイル(aperture 制御用)の入出力に対応
    2019.10.13: 上記による extract.c の出力仕様の追加に対応
    2019.10.17: coros.c で 2019/10/17 以前のデータにのみ修正を適用するよう修正。
    2019.10.19: extract.c で apt ファイル内の重み変更に対応
    2019.10.30: rmcr.c で差し引きスペクトルの生成方法を微修正
    2019.10.31: rmcr.c で第3引数を閾値の値から default 閾値に対する factor に変更
    2019.11.06: calib.c と extract.c で apt ファイル使用時にノイズが0となるバグを修正
    2019.11.06: klwarp.c でフラットに必要な最小カウントを 3000 => 1500 に変更
    2019.11.09: rmcr.c でサチるほど明るいスペクトルを CR と間違わないよう微修正
    2019.11.15: extract.c で夜光輝線強度を合わせる factor を apt file で指定できるよう修正
    2020.03.30: klwarp.c で、連続光がぼけてトレースしづらい場合のトレース性能を強化
    2020.05.02: extract.c に大気分散補正の機能を追加
    2020.08.20: klwarp.c で、変形マップの最適解探しを自動化、cal_VB.dat に Xe を追加
    2020.10.23: 新ファイバー用に klwarp.c, mkflat.c, rmcr.c, extract.c, calib.c を準備
    2020.10.27: extract.c, mkflat.c で各 fiber の実効効率の影響を反映
    2020.10.28: extract.c で sky 差し引き factor を 各 fiber の実効効率からの算出に変更
    2020.10.29: extract.c で大気分散補正のアルゴリズムを各 fiber の実効効率を考慮して修正
    2020.10.30: coros.c で 2019/10/17 以前のデータのクロストーク補正は露出61秒以下のみに修正
    2020.10.30: 新ファイバー用のソースを旧ファイバー用のソースと共通化
    2020.11.13: extract.c で大気分散補正強度に ALT 値に基づいた上限を設定
    2020.11.13: mkflat.c で twilight 情報を用いて flat 低周波成分を補正するモードを追加
    2020.11.15: mkflat.c, extract.c 間の伝達ファイル名に観測モード追加(mkflat_VB/VR.dat)
    2020.11.15: rmcr.c にマイナー修正(ファイバーグループ間ギャップのメジアン引きを最後に)
    2020.11.15: extract.c で、新バンドルでは大気分散の方向と強度を与えるのをデフォルトに
    2020.11.16: extract.c で、旧バンドルで ImR OUT を header で判断して NOR モード自動適用
    2020.11.24: klwarp.c で、新バンドル sky 用ファイバーに光が入らない場合にも対応
    2020.11.24: calib.c で、大気分散補正確認用自己相対スペクトルを calib.dat に出力
    2020.11.28: cal_V683.dat 追加、klwarp,mkflat,rmcr,extract,calib を V683 に対応
    2020.12.02: calib.c で、OBJ/REF 露出時間比を header の EXPTIME の重み平均から自動計算
    2021.03.26: 解析処理手順で、新ファイバーの sky fiber 部分を extract で無視する方法を追加
    2021.08.07: extract.c で大気分散強度の計算方法を少し変更
    2021.09.21: mkbias.c で、bias パターン異常条件を追加
    2021.12.26: klwarp.c で波長 fit 次数を7次に上げ同定輝線リスト3種(cal*.dat)にも輝線を追加
    2021.12.28: mkflat.c を微修正(スペクトル両端の収差が大きくなってきたのでそれに対応)
    2022.01.27: coros.c で一番上3本分の処理を微修正(2021.11-この日までのデータは再解析が無難)
    2022.01.28: klwarp 後の波長エラー確認用 gnuplot script として klwarp.plt 追加
    2022.03.08: mkbias.c で bias の縞パターンの多様性を許容するための微修正
    2022.03.12: mkbias.c で bias の縞パターンの多様性に対応するための大改修
    2022.03.12: rmcr.c で微小領域毎に独立して縞テンプレート fit するよう大改修
    2022.05.01: rmcr.c で画像端で NaN が発生し停止する事例に対応
    2022.05.19: 1000枚の画像から生成した bias1000.fits を追加
    2022.06.06: 800枚の bias 情報を追加した bias1800.fits で差し替え
    2022.06.06: 観測当日の bias を template と融合する modbias.c を追加
    2022.06.28: extract.c で 2022/06/27 以降の大気分散方向にロテータの 90°原点移動を反映
    2022.08.04: klwarp.c で CAL lamp がかなりサチっている場合の異常終了に対応
    2023.01.31: klwarp.c で不要な Warning を表示しないように微修正
    2023.01.31: calib.c でファイルリストの入力方法を微修正
    2023.09.26: mkbias.c で非常に低い確率で発生するバグ修正
    2023.09.27: cal_V495.dat 追加、klwarp,rmcr,extract,calib を V495 に対応
    2023.10.02: calib.c で 823nm の大気吸収を正しく補正するよう修正
    2023.12.10: mkflat.c で mkflat_**.dat に nan が入ることがないように修正
    2023.12.22: headerinfo.c を追加 (解析には関係ありません)
    2023.12.24: sky 有観測の際に obj-sky ペアを決定する obj_sky.c を追加
    2023.12.25: calib.c で標準星との airmass 差で大気吸収に違いが出るのを補正
    2023.12.26: apt ファイルを指定のアパーチャ形状に修正する modapt.c を追加
    2023.12.28: coros.c で、quadrant 4 の bit 7 不安定に対処
    2024.01.12: coros.c で、quadrant 1,2 の bit 6 不安定に対処
    2024.01.19: coros.c で、quadrant 1,2 の bit 6 不安定の対処を改良
    2024.01.23: coros.c で、bit 不安定の対処を再改良
    2024.01.28: coros.c で、bit 不安定の対処を再改良
    2024.02.04: coros.c で、bit 不安定の対処を再改良
    2024.02.16: coros.c で、bit 不安定の対処を大改良(多分これが最後)
    2024.02.23: coros.c で、Over Scan 領域の扱いを変更(3σclip + 1次 fit)
    2024.02.27: coros.c で、bit 不安定の対処を微修正
    2024.02.28: extract.c で、xmin > xmax で入力時に wlcen dwl 入力と解釈するよう変更
    2024.02.28: extract.c で、夜光輝線差し引きの際の位置の微調整を追加
    2024.03.03: Cal 画像から背景散乱光を差し引く bkgsub.c を追加、解析手順も修正
    2024.03.13: klwarp.c の軽微なバグを修正

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

●ポリシー

  • スペクトルを水平かつ間隔一定、波長は縦で揃えて分散一定となるよう変形(画像変形はこの1回のみ)
  • flat は空間周波数の高い成分のみを利用
  • 天体寄与の大きいスペクトルと小さいスペクトルを結合させて天体+背景光を差し引き CR 検知
  • bias 取得時に検出された縞ノイズ成分は CR 除去後にテンプレートで fit して差し引き
  • スペクトルは S/N 最大となるようカウントで重みをつけて積算(マニュアルで重み修正可)
  • ALT,INR の値から大気分散の量と方向を決定し、アパーチャをずらして積算

    概要はこちら (旧バージョンは こちら)

●下準備

  • over scan 差し引きと bad column 補正
    ./coros 生画像.fits 新画像.fits

  • bias 画像差し引き
    全 bias 画像を2σクリッピングで平均 combine、bias 以外の全画像から差し引く
    => mkbias で同時に縞パターン解析とテンプレート生成を行う
    (「その他気づいたこと」2番目の項目参照)

  • flat0 準備
    各 flat の平均を1に規格化
    2σ クリッピングをかけながら平均 combine (flat0.fits)

  • flatcal0 準備
    Hg,Ne の平均 cal 画像と上記 flat0(x1000 程度) を足しあわせ (flatcal0.fits)

●変形マップ作成

  • 画像変形の概要はこちら

  • 回転+歪曲・弓形補正マップ作製
    ./klwarp flatcal0.fits flatcal.fits
    (補正マップ初期値 klwarp_dx.fits klwarp_dy.fits が生成される)

  • 補正マップに3次修正項追加 (fiber 間隔 = 7pix)
    ./klwarp flatcal0.fits flatcal.fits klwarp_dx.fits klwarp_dy.fits cal_VB.dat

  • 補正マップに3次修正項再追加・波長方向3次修正 (1nm = 4pix)
    ./klwarp flatcal0.fits flatcal.fits klwarp_dx.fits klwarp_dy.fits cal_VB.dat

  • あと4回程度で一番良い補正マップになる
    ./klwarp flatcal0.fits flatcal.fits klwarp_dx.fits klwarp_dy.fits cal_VB.dat
    ./klwarp flatcal0.fits flatcal.fits klwarp_dx.fits klwarp_dy.fits cal_VB.dat
    ./klwarp flatcal0.fits flatcal.fits klwarp_dx.fits klwarp_dy.fits cal_VB.dat
    ./klwarp flatcal0.fits flatcal.fits klwarp_dx.fits klwarp_dy.fits cal_VB.dat

    Spec_dy: dyav=0.113393 dymx=0.204449 スペクトルの曲がりの残差平均と最大(pix)
    Spec_dx: dxav=0.137171 dxmx=0.350624 輝線スリット像の横ずれ残差平均と最大(pix)

    数値が悪化してしまったら klwarp を1回目からやり直していいところの1回前で止める
    (前回マップで変形 => 結果測定 => マップ修正と書き出し、なのでベスト値が表示されたらその次のマップが書き出されてしまい手遅れ)。

  • 上記作業を以下で自動的に行うようにした
    ./klwarp flatcal0.fits flatcal.fits cal_VB.dat

  • マップができたら別名で保存
    mv klwarp_dx.fits dx_VB.fits
    mv klwarp_dy.fits dy_VB.fits
    最終マップでベスト値が出ることを確認
    ./klwarp input.fits output.fits dx_VB.fits dy_VB.fits cal_VB.dat

    以降は任意の画像をこのマップを用いて形状補正できる
    ./klwarp input.fits output.fits dx_VB.fits dy_VB.fits


元画像 (flatcal0.fits)     形状補正後 (flatcal.fits)


dx 成分マップ (dx_VB.fits)   dy 成分マップ (dy_VB.fits)


pixel - 波長関係(左: VPH_blue, 右: VPH_red) (NIST のデータを使用)

●flat 準備

  • 形状補正
    ./klwarp flat0.fits flat.fits dx_VB.fits dy_VB.fits

  • 光源の波長特性とファイバー PSF 特性除去
    ./mkflat flat.fits flat2.fits

    mkflat のしていること

    1. 全ファイバーのスペクトルを平均一次元化し、波長方向に幅50pix でスムージング
    2. 1. の結果で各行を割る(波長特性除去)
    3. 各行をそれぞれの平均値で割り、平均を1に規格化(PSF 特性除去)

flat 画像 (flat.fits)     特性処理後 (flat2.fits)
表示レベル 黒: 0, 白: 5    表示レベル 黒: 0, 白: 2

●標準星解析

  • Sky 差し引き(あれば)
    星の写っていない Sky を取得している場合は、星が写っている場合と同様に CR 除去までを行ない、星が写っている画像の CR 除去前に差し引きする。

  • 形状補正
    画像毎に天体が写っているファイバーが変化するので、1枚ずつ処理する
    ./klwarp star.fits STAR.fits dx_VB.fits dy_VB.fits

  • flat 補正
    単に flat2.fits で割るだけ

  • Cosmic ray と Hot pixel 除去、縞パターン補正
    ./rmcr STAR.fits STARc.fits [threshold_factor]

    rmcr でやっていることは、

    1. 画像内の127本全スペクトルを(突発的なカウントの増加部分はとりあえず適当に切り落として)平均化した平均スペクトルを明るいもの40本とそれ以外のもの87本に分けて最大値を1つだけ除いて(これで残っている CR 1発は除ける)残りを平均、各line 毎に2種類を混ぜ合わせて残差最小となるように差し引く(127x7 line) => その後、もう1段階差し引きを追加
    2. スペクトル外の領域は、幅15pix のメジアンフィルター画像差し引きで強制的に平坦化
    3. 差し引いた画像でスペクトルが写っている領域の各 line (2048pix) の大きい方から 40 番目のカウントをリストし、更にそれら(127x7 個の数値)の大きい方から 40番目のカウントの2倍+10 を閾値として、閾値を超えている部分とその外側 2pix まででマスク作成(第3引数を追加した場合はそちらを優先) => その後、閾値の決め方も変更
    4. このマスクを用いてスペクトル差し引きからやり直し、閾値を決め直してマスク再作成(これを2回繰り返す)
    5. マスク部分を横方向マスク両端外側各3pixずつの median 値を用いて内挿して置き替え
    6. 2種類の縞パターンを用いて背景の縞を除去 (詳細は下の「その他気づいたこと」の2番目の項目参照)
    7. 確認用画像を rmcr.fits として書き出した後、初めに引いたスペクトルを足して元に戻す

    非常に輝線強度の強い明るい天体に適用する場合、天体の輝線が cosmic ray として処理されてしまう可能性があります。その可能性があると予想される場合は処理前後の画像で差し引き確認し、天体輝線部分が処理されている場合は第3引数として閾値に掛け合わせる値を与えて下さい。


標準星スペクトル (STAR.fits)     CR 処理後 (STARc.fits)

●切り出し

  • 1次元化と Sky 差し引き
    ./extract STARc.fits STARe.fits [xmin xmax] [APT_file]

    fiber 間隔は 7pix となっているので、{0.0, 0.4, 0.8, 1.0, 0.8, 0.4, 0.0} という重みをつけて足し合わせ。 各スペクトルの左から3つのデータは各スペクトルの代表値で以下の役割を持つ

    新ファイバーでの各 column, line の意味はこちらを参照

    column 1: 各スペクトルの重み (0: 一番暗いスペクトル、1: 一番明るいスペクトル)
         (各スペクトルの x=[xmin:xmax] 間での明るさ代表値で default は [101:1800])
    column 2: ノイズファクター (fiber 幅 7pix の重みの2乗和の√で、1.612)
    column 3: 規格化ファクター (ここでは column 1 と同じ値)

切り出し後 (STARe.fits)
    2048x323 pix の画像が出力される(上図の上端に3行追加された)。一番下の line から順に、

    line 1-127: ファイバー 1-127 を一次元化したもの
    line 128 : ファイバー 1-127 の均等な重みでの平均値
    line 129-255: 各スペクトルの重み(一番左の column 値)の大きい順に下から並べたもの
    line 256 : 空行
    line 257-264: 下から重み >0.2 >0.3 >0.4 >0.5 >0.6 >0.7 >0.8 >0.9
           のスペクトルを各スペクトルの重みを付けて平均化したもの。
           column 1: 重みの合計、column 2: ノイズファクター、column 3: 規格化ファクター
           (ノイズファクターは各スペクトルのノイズファクター×重みの2乗和の平方根、
           規格化ファクターは重みの重み平均値で、この値で各 line を割る事で重み1の
           スペクトルのカウントに合わせることができる値)。
    line 265-272: 下から重み <0.4 <0.35 <0.3 <0.25 <0.2 <0.15 <0.1 <0.05
           のスペクトルを均等な重みで平均化したもの。
           column 1: 積算スペクトル数、column 2: ノイズファクター、column 3: 1 と同じ
    line 273-280: line 257-264 に対し、ノイズファクターが半分以下となる sky スペクトルを
           line 265-272 のできるだけ上のものから選んで、夜光輝線強度が合うように
           微修正して引き算したもの。
           column 1: 引き算したスペクトルの平均値を最終的なノイズファクターで割った値
           column 2: ノイズファクター、column 3: 測光ファクター(line 257-264 の
           重み÷規格化ファクターで、この値をスペクトルにかけることで対応する aperture
           内での全カウント合計値に変換できる)
           column 1 の値は S/N に大体比例する値なので、この値の大きいものを最終結果と
           すべき

    line 281-320: スペクトルを 16nm 毎に切り分けて IFU 上のイメージにしたもの。
           一番左は明るさ評価をする波長範囲内(この次の項目参照)での積分値で、ピーク値が
           スペクトルの最大値と一致するように規格化したもの(単に見やすくするため)。
    line 321-323: 下から、S/N 最大となる Object+Sky, Sky, Object のスペクトル。sky は夜光輝線
           強度を合わせるための factor をかけたものとなっている。
    スペクトルが明るいと、分光器内部の散乱による隣のファイバーからの光の漏れこみで合成像にスパイクが入ったように見える場合があります。元画像のスペクトルのピークがサチっていないか確認して下さい。

    つづく

●観測天体の処理

  • 1次元化までは標準星の解析と全く同じ

観測天体スペクトル (OBJECT.fits)   CR 処理後 (OBJECTc.fits)
  • 1次元化と 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 が存在すれば自動で読み込んで適用する)

    ここまでくると、非常に弱い縞パターンが気になる。人工的に特定周期の縞画像(4つの 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 値)で重みをかけて足し合わせる。その際、対応するスカイ差し引き前のスペクトルとスカイのみのスペクトルもそれぞれのノイズレベルの2乗を掛けて足し合わせから平方根を取ることで、スペクトル内の相対ノイズレベル(最小値を1として規格化)として出力する。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 ダウンロードできない場合はこちらを参照)。第1列は波長(Å or nm or μm で自動判定)、第2列は 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% レベルでの事なので、通常の観測では特に問題とはならないと思われる。それよりも、補正できるはずの大気吸収がメジアンスムージングによって補正できなくなる方が気になる。試しに、異なる2つの標準星のデータを用いて、片方の結果をもう1方でキャリブレーションしたのが以下のグラフ。


標準星同士でのキャリブレーション  大気吸収部分修正後

    傾きが急なので縦軸を対数軸にし、上述の通りメジアンスムージングをかけてキャリブレーションしたものが、メジアンスムージングをかけずにキャリブレーションしたものがとなっている。大気吸収のある5つの波長帯の内、590, 628 はのままで問題無さそうだが、687,720,761 の3つの波長帯は明らかにの方が正しそうだ(761 は吸収が深いため airmass の違いで完全には補正できていないようだが)。この3つの波長付近では、距離 5nm 以内でメジアンスムージングなし、その外側距離 10nm まではメジアンスムージングなし→ありに線形変化するようにした(上右図線)。
    → その後、VPH-red で 823nm の大気吸収も修正が必要であることが判明し、これも加えて補正する波長帯は4つとなった。

    最後に、ノイズレベルを正しく出力するため、最終スペクトルを幅20nm で部分的に2次関数で fit し、3σクリップで外れ値を除去、再度残りで 2次関数 fit を3回繰り返してその位置でのノイズ値として最終スペクトルからのノイズ実測値を計算。このノイズ実測値を上記の画像から評価したノイズレベルで割り、その値を再度主要波長範囲で2次関数 fit、3σクリップで外れ値を除去、再度 fit を3回繰り返して、ノイズレベルと実測値の比を 2次関数として評価、この関数でノイズレベルを規格化することで、これまで最小値で適当に規格化していたノイズレベルが正しい値となった(即ち、大局的なノイズレベルは天体スペクトルから実測、局所的な変化はカウントの2乗積算値から出した相対ノイズレベルで決まっている状態)。下左グラフの線はこの方法で計算された 3σノイズレベルを表している。


観測スペクトル(大気吸収修正後) 標準星のスペクトル(大気吸収修正後)

●その他気づいたこと

  • 3つのファイバーバンドルの透過率のカラー特性が、標準星観測時とフラット取得時で結構異なる変化をしているようだ。Az の回転によるファイバーバンドルのねじれの変化が関係している可能性があるので、様々な Az 位置でフラットを取得することを2~3セット繰り返し、安定性をチェックしてみるべきだと思う。(Az 軸中心でのファイバーバンドル固定方法を色々試しているようですが、改善はしていません)

  • bias 画像に関して
    3.5925pix 周期の縞は 10枚の bias 画像平均を個別の bias 画像から差し引いたものでも確認できるが、縞の強弱が単なる唸りでは説明できそうにないので、10枚の bias の標準偏差フレームを作ることで縞が出やすい部分と出にくい部分を調べてみた。すると下図右のような画像になり、各 quadrant の左と右が繋がらないのでこれでは唸りモデルで再現できない。なぜこれほどまでに左端と右端が不連続になるのかの原因は全く不明だが、とにかく、唸りは使わずに強度変化を決めてやる必要がある。


  bias 画像          bias の標準偏差画像

    縞の成分には垂直な 3.6pix 周期のもの(下左図)と、傾斜した 7.8pix 周期のもの(下右図)の2種類があるので、この2成分(sin,cos の4枚)を10枚の bias 画像から自動で抽出するソフトを作った。

    ./mkbias biaslist bias.fits

    出力される画像は以下の通り。
    bias.fits: 入力 bias リストの2σクリッピング平均画像
    mkbias_[番号].fits: 2σクリップ後の各入力画像から上記平均を引いた画像
    mkbias.fits: 4つに畳んだ標準偏差画像に M パターンに沿った boxcar スムージングをかけたもの
    sin1.fits, cos1.fits: 3.6pix 成分縞テンプレート
    sin2.fits, cos2.fits: 7.8pix 成分縞テンプレート
    mpat.fits: M パターンマスク
    unity.fits: 変形後の M パターンマスクのカウント修正に使う全面1の画像

    上右図の bias 標準偏差に出ている "M パターン" を数値化した "M パターンマスク" を作り、この値が同じものに沿って cosmic ray 除去のタイミングで2種類の縞成分を除くことにした(実際は M パターンに沿った幅21pix で2成分最小2乗計算をし、結果を中央1列にのみ適用)。とりあえず、バイアス画像で動作確認 (以下は変形・フラット処理後の画像)。


3.6pix 周期成分縞画像      7.8pix 周期成分縞画像


平均引き・変形・フラット割り後のバイアス  縞軽減後      


3.6pix 周期成分振幅分布     7.8pix 周期成分振幅分布

    目で見ると 7.8pix 周期のほうが目立って見えるが、実際は 3.6pix 周期のものの方が倍以上振幅が大きい。細かい縞が見づらいこともあって、処理前後の画像の違いはまあ、こんなもんかなといった感じ。天体が写っていても大丈夫か確認してみる。


    上記縞処理方法で処理後の観測天体スペクトル(1次元化後)

    以前の処理法とほぼ遜色ない感じ。念のため、以前の手法も重ねて適用してどれだけ追加で除けるか調べてみた画像がこれ(上半分が適用前、下半分が適用後)。数値的にも zsdav: 4.87853 => 4.87477 でσ値がほとんど変化しないので、ほとんど縞が除けない状況になっている事が確認できた。最後に、できるだけ多くの bias 画像を用いて、縞の状態が過去どのように変化しているのか調べてみる。

    以下、5/5-5/16 の bias 画像の縞変化(傾斜は縦軸に対する傾き)

    日付M 傾斜3.6成分周期3.6成分傾斜7.8成分周期7.8成分傾斜
    5/52.285713.60469650.00000007.77968490.4799947
    5/61.896403.5149401-0.00060247.45974990.4853059
    5/71.822583.60400590.00419107.75789230.4756021
    5/71.929413.60427700.00335057.77907000.4692092
    5/82.100483.57915950.00000007.78715610.4905279
    5/914.40003.29055460.00000007.71218720.5625702
    5/101.547533.6028866-0.00085057.78238220.4559666
    5/102.098523.6034753-0.00084977.77695990.4723620
    5/112.246843.5559916-0.00269497.64233810.5436668
    5/143.044943.5781561-0.00157387.62215450.5893466
    5/1514.31253.6045019-0.00299587.77733890.4740607
    5/161.466423.6048799-0.00298787.77151490.4722291

    最も重要なのは 3.6成分周期の変化で、同じ日にこれが変化すると始めの縞除去のやり方の try & error での縞除去しかなくなり、自動化が難しくなる。今の所、同じ日なら大丈夫そうだが、観測の始めと終わりに10枚ずつ取得しておくほうが安全そうだ。上表で黄色の行は bias 取得枚数が5枚しかなく、2σクリッピングで異常値が除けない場合があるもの。また、5/15 の M 傾斜の異常値は、bias 標準偏差画像(4つにたたんでスムージングをかけた画像)に暗い横縞成分が入っているため。

  • 標準星のデータがない状態での calib
    観測中の同時解析などの際、通常は標準星を撮るのは後回しになるので、そういう場合のために標準的な "効率曲線" を calib から参照してやる必要がある。今の所、VPH_blue で取得された2つ分のデータしかないので、それらを平均したものを calib_VB.dat として書き出しておき、標準星が無い場合にはそれを参照して結果を出すように calib.c を修正した (VPH_red の場合はヘッダー内の波長情報から calib_VR.dat を自動選択するが、VPH_red で取得された標準星データがまだないので準備できていない)。


2つの標準星のカウント => Flux density 変換係数

  • rmcr での2成分 fit による差し引き
    非常に明るい天体を観測した場合、天体の写っているファイバーと sky のみが写っているファイバーのスペクトルが大きく異なるため、単に中間的な明るさのものを平均して スケーリングして差し引きするだけでは奇麗にスペクトルを消すことができない。127本のファイバーの内、明るいものから 40本とそれ以外の 87本でそれぞれ平均を取り、この2つのスペクトルの一次結合で各ファイバーのスペクトルを合わせてやることで、明るい天体の場合でも奇麗に消すことができた(←上下端の縞除去を行わない領域と比較することで、同時に処理した2つの縞成分も消えている事が確認できる)。


明るい天体のスペクトル (STARX.fits)  CR 処理後 (STARXc.fits)   

  • rmcr でのスペクトル外迷光の処理
    上図左上端にも薄く写っているが、ファイバースリット端から分光器室内の赤色光(780-820nm)が漏れ込んでいるようだ。分散素子を VPH-blue から VPH-red にすると中央付近に移動して更に明るくなり、CR 除去や縞パターン除去に影響を及ぼすので、単純な幅15pix (3.6 周期の約4倍で 7.8 周期の約2倍)のメジアンフィルター自己画像差し引きで、できるだけ縞情報を残して強制的に平坦に矯正する事にした。


  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 1と4の間にクロストークがある事が判明。
    以下は、明るい標準星の生データと 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. 左側3つの平均値 - 16 よりも小さい
    3. 左側2つの外挿値 - 8 よりも小さい
    4. 右側3つの平均値 - 16 よりも小さい
    5. 右側2つの外挿値 - 8 よりも小さい
    6. 左右2つの平均値 - 8 よりも小さい

    まず、左から上記 1,2,3,6 の条件全てにあてはまるものに 32 を足し、次に右から 1,4,5,6 の条件全てにあてはまるものでまだ 32 を足していないものに 32 を足す。これを2回繰り返す(最大で64加算される)。このやり方で修正したものが下図右、適用前のものが下図左。

  • 2019/10/17 に CCD 読み出しクロックの修正でこの問題が解決したため、ヘッダ情報を見てこれよりも古いデータにのみ上記を適用するように修正。但し、1分以上の露出の場合、宇宙線イベントが増えて(なぜか宇宙線イベントはクロストークが出ない)この補正が逆効果になるので、長時間露出には適用しないことにした。

●rmcr に1成分追加+閾値2次元化

  • CR 除去ルーチン(rmcr)では、明るいスペクトルのグループ(object+sky)と暗いスペクトルのグループ(sky)の2つの平均スペクトルを最小二乗条件で重ねて観測スペクトルから差し引く、ということをしているが、天体が非常に明るい場合やフラット取得時と観測時の装置状態の微小変化によるずれなどにより、残差が大きくうねっている事がある。これがあると CR 除去の際に誤って天体スペクトルを引っ掛けてしまう場合があるため、残差のうねりを除くため最後に以下の手順でもう1成分差し引くことにした。

    1. 横方向のσが最も大きい行を探す
    2. 他の行を上記σ最大行と順次掛けあわせそれぞれの行の平均を計算
    3. 上記の値が負のものは符号を反転させて最大値を1つだけ除き重ね合わせて平均する
    4. 上記残差テンプレートを用いて最小二乗条件で各行に合わせて差し引く

    これよる残差画像の変化が下図。暗い天体ではほぼ変化はないが、特に明るい標準星などの場合は以下のように改善される。

  • これまで、残差画像に対して定数の閾値で CR 判断を行なっていたが、夜光輝線のある場所や天体の輝線が強い場合には輝線の先端を CR と誤判定する確率が高くなる。それを避けるため、全体の閾値を上げるのもつまらないので、以下の手順で場所により閾値を変えることにした。

    1. CR 除去後の残差画像に差し引いたスペクトル+残差うねりモデルを加算して元に戻す
    2. スペクトルが写っている部分の各行と各列のカウントの平方根の平均を別々に出す
    3. 各列の平均値(横方向の変化)で全体の平均が1となるよう規格化後、1以下の部分を1にする
    4. 各行の平均値(縦方向の変化)で全体の平均が1となるよう規格化(こちらは1以下でもそのまま)
    5. 上記2つの値の積に√(スペクトルカウント平均)x30を掛けたものをその場所の閾値とする
    6. これにより CR 除去マスクを作り直し、1. から繰り返す(3回まわす)

    これにより、輝線などが縦に並んで写っている部分は上記3の情報で閾値が引き上げられ、天体が入っているスペクトルが明るい部分は上記4の情報で閾値が引き上げられることになる。但し、スペクトル両端部ではかなりカウントが下がり、標準星スペクトルでの誤検出を増やしてしまうため、1以下の部分を1にすることで波長方向への変化を抑えた。以下は非常に明るい天体(下左)と中程度の明るさ(下右)の天体の場合の閾値分布の例。

  • CR 除去効果の確認が毎回必要であることがわかったので、CR 除去後の画像の確認の他、CR 除去前後の画像を引き算して目で画像チェックする(以下「解析手順まとめ」で #CR 誤検出チェック として2箇所に追加済)。以下の例では、CR 除去後の画像(下左上)で default 閾値では黄色○部分に CR 引き残しがあることがわかる。前後の差し引き画像では夜光輝線の先端などが入っていない事がわかるので、閾値を default 値(rmcr を走らせた際に Threshold = で表示される値)の半分にするため rmcr の第3引数に 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 部分の重みは常に1なので変更しても意味なし、天体部分の数値は 99 とする必要はないが2桁だと上書きモードで修正できるので楽)。一番最後の夜光輝線強度を合わせるための 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 の一番上3行をこのアパーチャで計算したもので置き換える(下図)。この3行の意味は下から Object+Sky, Sky, Object に対応し、257-264,265-272,273-280行に相当するもの(先頭3列の意味も同じ)となっている(但し Sky は夜光輝線強度を合わせる factor をかけた後のもの)。calib.c では、この3行が追加されている場合はこの部分だけを参照するように修正している。

  • 夜光輝線が 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 ファイルにパラメータを書き出すことにする。手順は、

    1. 明るさ評価をする波長範囲内(中心波長 wpsf)での合成イメージにファイバー1本分の幅でスムージングをかけたものを psf とし、apt ファイルで "+" (または "x") で選択されている部分で切り取ったものを重み付きアパーチャとする。
    2. 適当な大気分散ベクトルを仮定し、上記 d(w)-d(wpsf) を掛け合わせて波長 w での大気分散によるアパーチャずれ量とする
    3. 上記の量で波長ごとにずらしたアパーチャを波長別の合成イメージにかけ合わせて空間的に積分し、全波長分足し合わせた際に最大値を与える大気分散ベクトルを探す

    この手法で高度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 で光が来ないファイバーにも対応(但し、3つのバンドルグループの端が消える場合には未対応)
  • 波長 fitting を3次から4次に変更


3次           4次

  • 少し収束が早くなるようにアルゴリズムを修正

    新ファイバーでの dome, twilight, dome/twilight の flat_VB2 画像(表示レンジ 0.7-1.3, 比は log で ±0.3 表示)。twilight でも Sky 用 fiber は透過波長特性が異なるようだ(光の当たり方が異なるのであれば問題だが...)。dome flat は twilight よりも更に差が激しく、やはり使えない感じだ。

    一次元化した dome, twilight, dome/twilight は以下のような感じ(比は log で ±0.5 表示)。視野中央の横並び3本は良くない。視野右端はケラレがあり、Sky用 ファイバーは flat lamp の光がちゃんと当たっていない感じだ。実は twilight flat も Sky fiber にケラレがあるようで、対処が必要。

●新ファイバー用の mkflat

  • sky 用ファイバーの flat に関して
    sky ファイバーは天体用ファイバーの位置からかなり離れたところにあり、拡散板を用いてもランプによる flat ではファイバーに対する光の入射角度に違いが出て、sky 用ファイバー部分は明るさだけでなく波長特性も異なる flat となりこのままでは使うことができない。twilight flat は、光の当たり方は同一になるが大気の吸収パターンや短波長側での光量不足など、波長により S/N が不足する状態となる。これを解決するため、mkflat で通常の flat の他に twilight の画像も読み込み、sky ファイバーの光量の違いと波長特性の違いを補正するモードを設けた。

    ./mkflat flat.fits flat2.fits twilight.fits

    mkflat を3引数で用いた場合には、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

  • 各ファイバーの重み付けに関して
    新ファイバーでは、直前のレンズによるケラレ等でファイバー毎にかなり実効効率が異なるため、以下の問題に対処する必要が出てきた。

    • ある程度より暗い光しか入らないファイバーは解析に加えない方がいい
    • 各ファイバーの平均カウントだけで天体がどの程度入っているかを判断することができない

    これに対処するため、flat (できれば twilight で取得したもの) を mkflat で処理する際に各ファイバーの相対効率を出力し(mkflat_VB/VR.dat)、その値を extract で参照することにした。

    各ファイバーの相対効率(平均1):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 値の最小値を0、最大値を1にスケーリングしたもの)
    各ファイバーの実効天体光含有率: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 の特性代表値を1つ増やし column 1~4 に以下の値を置くことにした。

    column 1: 各スペクトルの重み、すなわち実効天体光含有率
         (0: object 光が無いと思われるもの、1: 最も object 光が入っていると思われるもの)
    column 2: ノイズファクター (fiber 幅 7pix の重みの2乗和の√で、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: 相対効率の重み平均値(ノイズファクターは各スペクトルのノイズファ
           クター×重みの2乗和の平方根、規格化ファクターは重みの重み平均値で、この値で
           各 line を割る事で重み1のスペクトルのカウントに合わせることができる値)。
    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 ファイルで与えられた重み分布に対しファイバー1本分の幅でスムージングをかけてから指定されたアパーチャを適用して計算している。そのため、事前計算で算出された最大 S/N を若干下回る S/N となる。

  • 大気分散補正に関して
    大気分散補正は天体光の大半が視野内に入っていることが前提で、はみ出ている場合は大抵の場合大気分散0と判定される。その場合、大気分散による傾斜がスペクトルに残ることになるが、新バンドルでは視野が狭いためこの状況に陥りやすい。そのため、データごとに大気分散量を fit するのではなく、ALT, INR の角度から予想される値を用いて、視野外に出ている割合を波長毎に補正するようにした。

    以下は、従来の方法と上記測定値を用いて視野外に出ている割合を補正する方法の比較。同じ晩に南天にある同じ標準星を ALT〜(45°,50°,35°) で観測した結果を2種類の方法で解析し、それぞれの平均値で規格化した自己相対スペクトル(下図左&中:時間の経過とともに天候が悪化している)。ある程度は補正が働いている事が確認できるが、大気分散方向で大きく視野外に出ていると補正が困難になる。下図右は別の日に同じ天体を視野中央に入れて新方式で解析した結果。少なくとも高度角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: ↑ とカタログスペクトルから算出した効率(%)で、大体カウントに
           比例する重みで積算されているため、真の効率の6割程度の値であることに注意

    旧バンドルでの大気分散の大きさと方向を確認するため、2020/4〜9 のいくつかの標準星データを用いて確認してみた。大気分散の強度や方向の fit 結果にはかなりのばらつきがあるが(気圧や湿度の違いによるばらつきがどの程度効くのかは不明)、旧バンドルと新バンドルとの見かけの分散強度比が 1/1.16:1/0.81 であることも考慮すると以下のようになる(天頂角<30° の点は大半が分散0に集中)。ImR を使わない場合の NOR を付けて extract を走らせるモードでは、これまでは方向のみを決めて強度は個々のデータの fitting で決めていたが、方向・強度どちらも下図緑線で決めたものを使うことにする。即ち、旧バンドルでの NOR モードと新バンドルの通常モードはどちらも大気分散の方向と強度として緑線で示された測定値を用いることにする。また、ImR の IN/OUT を header で判断して、OUT の場合には extract が NOR モード相当の動作をするよう改訂した。逆に fitting したい場合には FIT と入れることで大気分散の方向と強度を fitting で決定する(新旧バンドル共通)。

●calib で観測天体と標準星の露出時間比を自動計算

    これまで、観測天体と標準星の露出時間比を手動で与えていたが、露出時間の異なるデータが混在する場合にも簡単に対応させるため、EXPTIME 情報を header から読み取り、スペクトル重ね合わせの際の重みを用いて平均し、露出時間の比を自動評価するようにした。以後、第3引数を省略した状態で走らせる使い方を 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、波長分解能の違いは2倍強程度だが、V683 の方は 650nm 以下はほとんど使えない感じ。ほぼ同じ露出時間でもかなり S/N に差があり、V683 はかなり損な印象だ。

●露出時間が短い場合のカラー安定性に注意

  • 2つの標準星で比較
    Feige110 (高度角50°, 露出時間60秒)と HR153 (高度角70°, 露出時間1秒) の比較。どちらも7枚ずつ取得し、それぞれの平均スペクトルで規格化して plot したものが下左図(Feige110, HR153)。高度角が高い HR153 の方がカラーが安定していないことがわかる(3つ上の項目中のグラフよりも縦軸レンジが狭いことに注意)。

    それぞれの31枚の IFU 上の波長別合成イメージをそれぞれの波長帯での平均カウントで規格化し、一番長波長側の1枚を除いた残りの30枚を10枚ずつに分け、RGB 画像にしてみた。下左が HR153, 下右が Feige110。やはり、露出時間が短いと PSF が安定せず、それがスペクトルの color 変化を引き起こしている可能性が高い(PSF がいい瞬間が最もスペクトルが青い)。こういう場合は S/N 最大条件でアパーチャを決めるのではなく、できるだけ重み一様で大きく広げたアパーチャにすべきとは思うが、自動判定はなかなか難しそうだ...(試しに手動で重み一様かつできるだけ広げたアパーチャにしたものが上右図)。

    ところで、高度角70°でも結構大気分散が大きいことが気になった。大気分散は天頂角の tan に比例するので(高度が低くなると tan 3乗の成分も出てくるが)、KOOLS での観測結果に合わせて天頂角の2次関数(下図線)で決めていた大気分散は、やはり天頂角の tan (下図線)で決める方がいい気がしてきた。(下図左は新ファイバー、右は旧ファイバーのデータ)

    この大気分散強度を用いて extract からやり直してみたのが以下。

    まあ、以前のものと似たり寄ったりで良くなった部分もあれば悪くなった部分もある。一応、こっちの方が理論的にも正しいはずなので、こちらのソースで差し替え(再解析が必要なほどではないと思う)。そもそも、HR153 の方は、ALT+INR の値では、青い光が4時の方向に伸びていないといけないはずが、合成イメージでは6時位の方向になっていてちょっと大気分散の方向が変(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-445nm445-505nm505-565nm565-625nm
Hg
Ne
Xe
625-685nm685-745nm745-805nm805-865nm
Hg
Ne
Xe
    一番長波長側の Xe の強い2本の輝線のピークが見えていないが、この2本は各段に強度が強く完全にサチっているので無視。それにしても Xe は Hg,Ne と相性が悪い。どうしてこうなるのか...

    ここで、NIST のデータは default では空気中の波長が表示されている事に今頃気が付いた... う~ん、ここはじっくり考えないといけないところ...

  • fitting 次数の変更

    これまで、通常の分散素子の分散は3次程度で歪曲の次数も高くても5次だろうという予想から最大で6次までしか確認していなかったが、大塚さんより7次で fit できるという情報をもらったので klwarp.c の大改修を行った。fitting 次数を上げると少しの揺らぎでも大きく影響が出るので、元となる輝線データの精査・追加と検出数が足りない場合に次数を5次に落とす等、様々な例外処理を加えて何とか安定して動くようになった。また、NIST も理科年表も空気中での波長で表示されるため、同定する輝線リストは空気中での波長を参照し、読み込んだ直後に klwarp.c 内部で真空波長に変換している。まず、NIST で使われている空気の屈折率を確認するため、default で表示される空気中での波長と、設定変更して真空中の波長で表示した場合の比(下図線)から NIST で使われている空気の屈折率を確認した。Owens 1967 の式で、気圧 1012.5hPa の乾燥空気で計算(下図線)すると大体合うことが確認できた。

    とにかく何とか 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 のもので、上左図の青線。この日の4次 fit の結果はちょっと他とは異なる振る舞いをしているが、旧ソフトでは自動では同定できない輝線があったようで、その影響で輝線がない 480nm,800nm 付近での振る舞いが変わるようだ。

    全体が上っているのは、NIST のデータを真空中の数値に変換した影響(これまでのものは波長が短くなってしまっていた)で、これまでのデータの波長を 0.2nm 増やせば、波長帯中央付近に関しては大体正しくなる。また、実際のデータが短波長側で折れ曲がっているのは、fitting 次数を7次に上げたことにより参照した較正用輝線の外側では急速に値が変化するため、波長補正を最後の輝線の少し外側で打ち切っていることによるもの。

    新しい klwarp を走らせた際に出力される波長較正確認用の klwarp.dat の内容は以下の通り。

    1行目:中心波長、Δλ/pix、有効波長範囲の画像上での右端の位置
    2行目以降:
     1列目:輝線を画像上で探す際の初期値 [pix]
     2列目:空気中での波長 [nm] (内部では全て真空波長で計算している)
     3列目:輝線の種類
     4列目:中央のファイバーで、この輝線の変形後の画像(flatcal_VB.fits など)上での位置 [pix]
     5列目:中央のファイバーで、この輝線の目標位置(1025-(波長-620)/0.25)とのずれ [pix]
     6列目:この輝線の輝線位置の横方向のばらつき [pix] (0.3pix 以上のものは較正には用いない)
     7列目:この輝線を検出できたファイバー数 (ファイバー数の半分以下のものは較正には用いない)
     8列目:この輝線の目標位置からのずれの標準偏差 [pix] (1を超える場合は同定ミスの可能性大)

    1,4,5 列目の値が nan となっているものは、中央のファイバーではこの輝線のずれが大きく、クリッピングされてしまったことを表している。較正に用いなかった輝線は、先頭に # が付いてコメントアウトされた状態になる。この klwarp.dat の4列目と2列目、5列目×0.25 を plot したものが、上の pix - 波長関係図。

    また、この作業の際、CAL 画像に強い宇宙線がスペクトルに近い角度で入っていると、宇宙線につられて初期のトレースで失敗する事例を発見したので、CAL 画像の重ね合わせを median ではなく minmax rejection (最大値1つを削除) での平均の計算とした。

●縞パターンの変化

    いつ頃からか、検出器に現れる縞のパターンがかなり変化してきた。以前は検出器の読み出しクロック起因の高周波ノイズと思われるほぼ垂直な細かい縞が出ていたが、それがなくなり、左右45°位の角度の縞が4つの quadrant 間の強度相関が無いような感じで現れる。今後のことも考え、縞パターンが変化しても対応できるように mkbias と rmcr の大改修を行った。

    まずは、bias の連続画像平均からの差分画像を4つに畳んで、bias 取得中の縞を見やすくする(下図一番左)。水平な細かい縞と左上から右下方向への縞があることがスムージングをかけると確認できる(下図左から2番目)。スムージング前の画像を2次元フーリエ変換してパワースペクトル(の平方根なので振幅スペクトル、下図左から3番目で表示は log スケール)を調べてみる。様々な縞成分があるが、中央上下端にある水平縞成分と、右下がり縞に対応する右上方向に幾つかの周波数で強い成分が存在することがわかる。bias は30枚あるので、それぞれで計算したスペクトルを重ね合わせてみる(下図左から4番目)。すると、上記の成分だけが共通して存在し、他の成分は取得中に変化した感じだ。この情報から2つの縞テンプレートを作成するが、水平方向の縞はスペクトルに影響を与える可能性があるので、画像上で水平から±30°の角度内の縞は無視して、強度の強い順に2つの成分を取り出す。quadrant 間の強度相関が無くなったようなので、これまで mkbias で作成していた強度相関の参照ファイル (mpat.fits と unity.fits) は不要となった。

    取り出した2つの成分が以下の通り。これを画像変形して rmcr での処理で用いる。
    rmcr では、まず天体スペクトル成分を差し引き、宇宙線ノイズを除去(下図右)する。

    この画像に画像変形・フラット割り後の縞テンプレート2種(sin,cos パターンあるので4枚)を 30x30 pix を単位として最小2乗条件で4成分の強度を算出、これをずらしながら繰り返して4成分の強度マップを作る。これから縞成分を合成したもの下図左だが、表示スケールを合わせると非常に薄い縞となり、大して除けていない(というか、そもそも bias 取得時の縞パターンと観測中の縞パターンが異なる...)。まあ、悪化はせず若干の改善という程度だ。縞成分を差し引き後の画像が下図右。上図右と比べても大して変化はしていない...(前後の画像を拡大してブリンクすると、テンプレートパターンの縞が軽減されていることは確認できるが)

    まあ、やっていることに間違いはないと思うので、観測中の縞とバイアス取得時の縞パターンが同じであれば、ちゃんと除けるのではないかと思う。

●標準 bias

    大塚さんが 1000枚の bias を撮ってくれたので、早速それで標準 bias (bias1000.fits) を作ってみた。一応、当日の bias と比較して、形状に変化無ければこちらを用いることにする。以下、これまでの 20枚 bias との比較。もっと早く取得しておくべきだったと思う...

    3週間後に磯貝くんが 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 パターンの縦縞の位置は観測日の検出器の状態により変化するようだ(2つ上の図の左右比較)。低空間周波数成分は観測日のものを反映させる必要がある事がわかったので modbias.c を追加した。

  • 観測当日bias (bias0.fits) と標準bias (bias1800.fits) の融合
    ./modbias bias0.fits bias1800.fits bias.fits

    modbias が行なっていること

    1. 観測当日bias - 標準bias を計算
    2. 1.の結果に 3x3 median smoothing をかける
    3. 2.の結果の4つの quadrant を、平均値を計算しながら左右反転重ねあわせ
    4. 3.の結果を 11x51 範囲で boxcar smoothing (領域端では z 歪みを無くすため幅を減らす)
    5. 各 quadrant の平均を再現させながら4.の結果を4つの quadrant に展開
    6. 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 でファイバーの空間方向の重みを全て1とし、apt ファイルで最も光量の多いファイバーに対し 5% 以上の光が入ったファイバーを全て同じ重みで足しあわせることでできるだけカウントを減らさないようにして評価してみた(conversion factor 1.389)。以下は大気の透過率も含む VPH-blue(4回), VPH-red(3回) での測定値。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 かで番号の小さい方か大きい方かを決定する。

●APT ファイル自動修正

    これまでは、point source の計測を前提にしていたため、アパーチャの形状にはあまり気を遣わずに S/N 最大となる事だけを考えてきたが、apt ファイルを修正して決められた半径での固定アパーチャでの結果が出せるように apt ファイルを自動修正するソフトを準備した。

    ./modapt aptfile 半径 タイプ [-]

    半径の単位は arcsec で、タイプは B(Box)/G(Gauss)/C(Cone) のどれか1つ、"-" は Sky 領域無し(予め Sky が引かれている場合に用いる)となる。Box は半径内が1で外が0、Gauss は半径内が Gaussian の1σまでと同じ形状で外が0、Cone は半径内が1~0.5で直線的に変化し外が0、という重みのアパーチャ。中心は与えられた apt ファイルでの重心で決定する(各ファイバー形状は六角形としている)。指定されたアパーチャに対し、各ファイバー内で三角形グリッド上の約3万点での重み評価を行い、その平均でファイバー1本の重みを決定する。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枚程度の連続露出の平均で、は連続して観測が行われている(但し、は観測開始直後)。連続観測の最中ですら、波長帯両端は全然安定しないことがわかる。

    この原因を調べるために、1次元化したデータの平均に対する比を見ていくと、観測開始付近でかなりの長時間検出器が暴れていて、下左図のような状態になっていることがわかった。一番右の quadrant の非常に暗いカウントが入っている辺りだけが選択的にダークが増えたような状態になっている(再生像がリング状に!)。これは、辿っていくと生画像から入っているため(下下図)、検出器側の問題(ビット不良の一種?)のような感じだ...

    詳しく見ていくと、quadrant 4 の bit 7 だけが、1088ADU (2進数で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 (2進数で0000010001000000) に近づくと、bit 6 が0になるという bit 異常が発生していることが判明した。基本的には上記と同様に対応するが、一度誤った修正が始まると芋づる式にどんどん繋がって伸びてしまうため、完了後に逆側から確認して間違った補正は元に戻すようにした。非常に運悪く cosmic ray がこの地帯の両端に入ると、その間の誤補正が修正できないが、可能性は低いのでこれには特に対処しないこととする。

    気温の低下に伴い bias のカウントも下がってきた感じで、とうとう quadrant 2 の bias のカウントが危険領域に入ってしまった。この状態では、オーバースキャン領域でも bit 異常が発生するため、オーバースキャンの値の利用の仕方によっては注意が必要だ。オーバースキャン領域(quadrant 1,2 と 3,4 の間の部分のみ)では、bit 異常が発生するカウントでかつ左右の pixel のカウントに比べて異常値のジャンプが見られた場合は即修正、ということにすることでほぼ解決できた。

    2つの宇宙線が同一 quadrant の同じ列に入った場合に、その間が誤判定されてしまうことに対処するためにかなり複雑なアルゴリズムで対処した。あまりに複雑でわかりやすく説明できないのでここでは割愛。この問題は、bias レベルを 1200 まで引き上げることで解決したようなので、観測日が 20231001~20241118 の間のデータにのみこのアルゴリズムを適用することにした。

    ...と思ったが、どうやら 1088+128=1216ADU (2進数で0000010011000000)に近づくと bit 異常が発生する場合があるような感じなので、20241118 以降のデータに関してはこのレベルに近づいたら注意ということで bit error を修正することにした。coros を走らせた際、bit error が100個以上検出された場合に8個の整数が表示される。左側の4個が順に quadrant 1〜4 のオーバースキャン領域で補正した bit 異常の数で、右側の4個が画像領域での bit 異常の数を表している。安定している場合は全て0になるので、これら8個の整数は表示されない。

    更に、64 と 32 の bit 異常が混在して現れる現象が確認されたので(結構な頻度で混在することが判明)、これにも対処できるように修正する。幸い、bit 異常はなぜか以下のパターンでしか現れないので、「bit 異常境界値」さえ決定しておけば(2024/1/18 より前は 1088、以降は 1216)そこから±(32~64) の範囲さえ探せば良く、周辺の状況から 32 なのか 64 なのかを判別できる。ただ、32 や 64 の段差は、状況によっては正常な画像であってもたまたま起こりうるため、bit 異常による段差なのか、輝線や連続光スペクトルの裾野や弱い輝線のピークで起こっている事なのかを判別する必要がある。

    また、ほぼ正常な画像には bit 異常修正はしない方が無難なので、第0段階としてどのタイプの bit 異常がどの程度出ているのかを各 quadrant で解析し、優先度の高いものから2タイプだけを修正することとした。優先度は -64,+64,-32,+32 の順。また、ほぼ全ての pixel で bit 異常が出ている場合は、1タイプに統一される傾向があるため、それも考慮して各 quadrant でどのタイプの補正を行うのかを決定してから補正を始める。

    1.5ヶ月にわたり、2023年11月から2024年2月までの手持ちの1000枚以上の画像全てに適用できるようにほぼ毎晩改良していったので、最終版はソースが複雑になりすぎて時間が経ったら自分でも理解不能になりそうなレベルにまで成長した。修正手順の概要は以下の通り。

    1. 各 quadrant にどのタイプの異常がどの程度出ているかを大まかに調査、2タイプまで選択
    2. bit 異常候補を全て一旦修正
    3. 通常のオーバースキャン領域による bias 合わせと quadrant 間の gain 修正
    4. 輝線検出と輝線周辺部で bit 異常修正結果を再修正(輝線周囲は判断アルゴリズムが異なる)
    5. その他の領域で bit 異常修正結果を再修正
    6. 全体で再修正結果を再々修正

    という感じで、初めに可能性のあるものを全部修正し、後から様々な条件でその結果を見直す、という流れとなっている。

    ./coros 生画像.fits 新画像.fits +

    で、coros.fits も生成される。この画像は bit 異常をどのように処理したかを表すファイルで、1 は高優先度で bit 異常処理されたもの、2 は低優先度で bit 異常処理されたもの。10 を超えているものは、一旦処理されたが、後に元の状態に戻されたもの、20 を超えたものは輝線周辺での異なるアルゴリズムで判断したもの、を示している。多分、私以外は見ても良くわからないと思う。

    coros を走らせると入出力ファイル名の後に以下のような表示が出る(実際は改行なしに1行表示)。
    この例では、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 の補正を行うと、上下の輝線の間が繋がってしまう場所が発生するが、スペクトル抽出の際はその部分はほぼ重み0となるため、他の部分の補正が行われることを優先している。

    2024/3/5 に読み出し回路のボード交換が行われ、bit 異常問題は解決した。この日以降に取得されたデータに対しては coros.c では bit 異常処理は行わない。

●Cal lamp 取得時の背景光

    VPH495 を使うと、背景に原因不明の連続光が発生する。slit 以外の部分も光っているが検出器全面ではないので、slit unit 付近からの光の漏れ込み(光源はファイバーバンドル自身?)があるようだ。VPH-blue や VPH-red ではこの成分は弱いのでこれまではあまり気にしなかったが、VPH495 ではこの光が flatcal でのスペクトルと輝線位置のトレースに影響を与えるので、ソフト的に大まかに差し引くことにした。

    ./bkgsub cal画像

    で、輝線が写っている行と写っていない行を自動識別して、上下の輝線が写っていない行の1次 fit を内挿して差し引く。これで、VPH495 でのトレースが安定した。

●解析手順まとめ

    以下、"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)が表示されます

  1. 解析 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 -Wall
    gcc modapt.c -o modapt.c -lm -Wall

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

  3. iraf & ds9 起動(以下、cl シェル内での作業)
    cl
    !ds9 &
    cd 2021XXXX

  4. 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

  5. Sky を OBJECT とは別に撮っている場合
    #Sky が OBJECT と同じ名前になっている場合は観測ログファイル修正
    !grep OBJECT obs.log
    !../headerinfo ./ | grep OBJECT
    #座標から OBJECT と Sky の番号を確認
    !emacs obs.log &
    #Sky は天体名に sky をくっつけて OBJECTsky としておく

  6. 観測ログファイルとファイル名の統合
    !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 番号が揃っていることを確認

  7. 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

  8. 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=sigclip lsigma=2 hsigma=2
    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=sigclip lsigma=2 hsigma=2
    displ twilight_VB0 1

  9. 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

  10. 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
    !../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
    # ↑一番右列の値が1を超えている場合は同定ミスをしている可能性あり
    # 該当する 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

  11. 縞除去用ファイル準備
    !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

  12. 標準星処理
    !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
    #閾値を2倍にする場合
    #!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 誤検出チェック(スペクトルを引っ掛けていたら↑閾値2倍でやり直し)
    !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
    #新バンドルの sky fiber を強制的に無視する場合は、mkflat_VB.dat (twilight or dome flat で
    #得られた相対カウント)の 78-84番の値を 0.1 以下にして下さい。結果は apt ファイルで確認。
    !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

  13. 観測天体処理
    #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 画像があれば以下3行追加
    #!../obj_sky OBJECT OBJECTskyc
    #imarith @OBJECT - @OBJECTskyc @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 など第2引数を小さくした場合は 600±60nm と解釈される
    !sed -e 's/OBJECTc/OBJECTe/g' OBJECTc | sed -e 's/$/ 512 1536/g' > tmp2
    !paste -d" " tmp1 tmp2 > tmp.sh
    #新バンドルの sky fiber を強制的に無視する場合は、mkflat_VB.dat (twilight or dome flat で
    #得られた相対カウント)の 78-84番の値を 0.1 以下にして下さい。結果は apt ファイルで確認。
    !sh tmp.sh
    !ls OBJECTe_*.fits > OBJECTe
    imexam @OBJECTe 1
    #Sky frame を先に差し引いた場合は apt ファイルから sky を全削除して再度 extract を走らせる
    #以下は、半径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/OBJECTe1/g' OBJECTc > tmp2 #!paste -d" " tmp1 tmp2 tmp3 > tmp.sh
    #!sh tmp.sh
    #imexam @OBJECTe 1
    #標準星のカタログスペクトルを STAR.dat として保存
    #STAR.dat の第1列は波長(Å or nm or μm)、第2列は 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 + @OBJECTskyc @tmplist
    imarith @tmplist + @OBJECTskyc @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/OBJ_VBe/tmpliste/g' OBJ_VBe > tmpliste
    !../calib tmpliste STARe STAR.dat | cut -f9 > tmplist.dat
    !paste OBJECT.dat tmplist.dat > tmp.dat && mv tmp.dat OBJECT.dat
    !rm tmplist*

    これにより、
    #column 10: 観測天体の最終スペクトル中の相対ノイズレベル(差し引き Sky 情報込み)
    が追加される。



iwamuro@kusastro.kyoto-u.ac.jp