この前は国立天文台で直接撮像法で系外惑星を探す研究について勉強する機会があって、色々勉強できたのでこの記事で纏めてみたいと思います。

以下は直接撮像法で系外惑星を探す流れを手短に説明します。

興味があったら誰でもデータをサイトからダウンロードして説明の通り行えると思います。

作業は大体pythonのコードで済ませます。

ここで説明すること

直接撮像法について説明
必要なソフトとモジュールのインストール
データのダウンロード
vlt-sphereモジュールとesorexでデータ解析
ADI手法でデータを合成する

直接撮像法とは

系外惑星を探す方法は色々あります。
- 直接観測法 (imaging)
- トランジット法 (transit)
- 視線速度法 (radial velocity)
- 重力マイクロレンズ法 (microlensing)
- その他

これらの方法についてこの前の記事で説明してあります。
https://qiita.com/phyblas/items/ec29c49302c7b0a60905

その中で直感的で一般の人にはわかりやすいのは直接観測法です。つまり、直接系外惑星の写真を撮ることです。

ウィキペディアとかで直接観測法で発見された系外惑星を検索したら、例えば

これらのページでは実際に望遠鏡によって撮られた写真が載っているのです。

勿論惑星はたった小さな点にしか見えないのですが、それでもちゃんと写っているのです。

『百聞は一見に如かず』と言われているのだから、直接観測法は一般人にとって印象が深い方法だと思えます。

ただし実際に、惑星は主星と比べたらとても小さくて暗いものなので、惑星を撮るのは簡単なことではないのです。この方法で見つかるのは主に随分主星から遠く離れた巨大惑星です。

一般に、コロナグラフで主星の光を遮って、ADIという手法を使います。これについてあとで説明します。

使うデータ

ここで例として使ってみるのはVLT（Very Large Telescope）という望遠鏡によって撮られた写真データです。

VLTはチリののアタカマ砂漠の中のパラナル天文台に建設された望遠鏡で、経緯度は24°37′38″S 70°24′15″W，海抜2635メートル。

VLTには色々な装置が入っています。それぞれは違う目的のため準備されたのです。その中で、直接撮像法による系外惑星探索を目的として作られたSPHERE（Spectro Polarimetric High contrast Exoplanet REsearch）という装置が含まれています。

SPHEREの中でも主に3つの部分からなされています。

IRDIS（Infrared Dual-band Imager and Spectrograph）同時に2つの違う波長で写真を撮る
IFS（Integral Field Spectrograph）光を違う波長に分別し3次元の画像を作る（空間2次元+波長1次元=3次元）
ZIMPOL（Zurich Imaging Polarimeter）偏光の機能を持っている

ここで一番簡単であるIRDISのデータを使います。

データ解析について

普段は撮影された天文写真は使う前にまず補正を行わなければなりません。これは天体データ解析という段階です。

これについて、日本語の資料はすでに充実しているので、は深く触れません。

例えばこれらのサイトを参考に

補正をするためには、対象の写真と共に必ずダーク画像やフラット画像などを撮る必要があります。

ESOサイトからダウンロードされた画像データはまだ処理されていない生データですが、補正に使うダーク画像やフラット画像も準備されているので、すぐ一緒にダウンロードしてデータ解析できます。

データ解析の方法は望遠鏡や装置によって違うので、普段は解析しやすいようにパイプラインというものが準備されています。

装置によってパイプラインは多少違いますが、共通する部分も多いです。

VLT-SPHEREのパイプラインについてはここで説明してあります。
http://www.eso.org/sci/software/pipelines/sphere

ただ、もっと簡単のためにpythonで自動的にVLT-SPHEREのパイプラインを操作するモジュールが準備されています。

それはvlt-spereというモジュールです。

これを使うと、あまりデータ解析について詳しくわからなくても解析できるほど簡単になるとも言えると思います。

vlt-spereモジュールを使うにはまずesorexというソフトが必要です。

esorex（ESO Recipe Execution Tool）はESOによって準備された天文写真処理用のソフトで、主にESOの各種の望遠鏡によって撮影された写真を処理する機能を持っています。

必要なソフトのインストール

主にpythonを使うので、必要なモジュールをインストール必要があります。

それと、それとESOのサイトからダウンロードされた画像データを解析するために使うesorex。

esorex

esorexはSPHEREパイプラインの中に同梱されています。

macを使うこともできますが、linuxの方が問題なくインストールできるようです。

インストールファイルはこのリンクでダウンロード >> spher-kit-0.40.0.tar.gz

そして凍結してインストールします。

tar xzf spher-kit-0.40.0.tar.gz # 解凍
cd spher-kit-0.40.0 # フォルダに入る
./install_pipeline # インストール

どこかにインストール先フォルダを指定する必要があります。

インストールの時、gasganoというソフトを一緒にインストールするかって聞かれます。ここでは使う必要がないのでインストールしなくてもいいのです。

時間が経ったら、リンクが使えないとか、もっと新しいバージョンが出る可能性があるので、その場合はSPHEREパイプラインの説明のページを調べたらいいです。

インストールが終わったら、ターミナルとかでesorexというコマンドを実行できるように、PATHの設定が必要です。

~/.bashrcファイルにこれを追加

export PATH="インストール先フォルダ/bin:$PATH"

そして改めてターミナルとか開いてesorexというコマンドを実行して、問題がなければインストール終了となります。

esorex

     ***** ESO Recipe Execution Tool, version 3.13.2  *****

Libraries used: CPL = 7.1.1, CFITSIO = 3.45, WCSLIB = 4.24, FFTW (normal precision) = 3.3.4-sse2, FFTW (single precision) = 3.3.4-sse2, CPL FLOP counting is unavailable, enable with -DCPL_ADD_FLOPS, OPENMP = 201511

pythonモジュール

astropy: 天文学的計算
astroquery: 天文データのダウンロード
skimage: 画像処理
vlt-sphere: SPHEREパイプラインの操作
vip_hci: ADIの実装

それと、言うまでもないかもしれませんが、numpy、pandas、matplotlibも必要です。

全部pipを使って簡単にインストールできます。

pip install astropy
pip install --pre astroquery
pip install skimage
pip install vlt-sphere
pip install vip_hci

データのダウンロード

データをダウンロードする方法について説明します。

ESOサイト

VLTはESO（ヨーロッパ南天天文台）によって運営されている望遠鏡です。ESOの観測データは基本的にESOのサイトで公開されているから、全世界の研究者や学生は簡単にダウンロードして、研究や勉強に使うことができます。

ただし一部の比較的に新しいデータは関係者が優先に使用するものだから非公開とされています。ダウンロードできるのは数ヶ月前か一年前からのデータです。

このリンクでデータを検索してダウンロードすることもできます。
http://archive.eso.org/eso/eso_archive_main.html

ただ、pythonを使うのならもっと便利な方法があります。それは、astroqueryを使って検索するという方法です。

データをダウンロードするためにはまずESOのサイトでユーザー登録をしてログインする必要があります。
https://www.eso.org/UserPortal/pre-authenticationSecure/registration.eso

astroqueryでデータをダウンロードする

astroqueryは色んなサイトからの天文データをダウンロードするために使える便利なモジュールです。

ESOのデータはastroquery.esoというサブモジュールで取得できます。

まずはSPHEREのIRDISによって撮られた星のデータの一覧を見ます。

ここでHR8799という星を対象として選びます。

HR8799は今まで直接撮像法で惑星が4つも見つかった星なので有名です。

query_instrumentというメソッドを使ったらどのようなデータがあるか調べられます。

from astroquery.eso import Eso

eso = Eso()
eso.ROW_LIMIT = 10000
cf = {
    'target': 'HR8799', # ほしい天体の名前
    'dp_type': 'OBJECT', # 補正のためのデータ（ダークやフラット）ではなく、対象自体
    'dp_tech': 'IMAGE%',
    'seq_arm': 'IRDIS'
}
df = eso.query_instrument('sphere',column_filters=cf).to_pandas()
df = df.set_index('DP.ID') # DP.IDをindexに設定
print(df)

結果

                              Release Date       ...       DIMM Seeing-avg
DP.ID                                            ...                      
SPHER.2014-07-13T08:42:55.582   2017-04-28       ...                   N/A
SPHER.2014-07-13T08:50:37.101   2017-04-28       ...                   N/A
SPHER.2014-07-13T08:54:08.958   2017-04-28       ...                   N/A
...
（長いので省略）
...
SPHER.2018-08-20T02:59:58.413   2019-08-20       ...           0.36 [0.01]
SPHER.2018-08-20T03:00:19.458   2019-08-20       ...           0.36 [0.01]
SPHER.2018-08-20T03:00:39.605   2019-08-20       ...           0.34 [0.01]

[115 rows x 33 columns]

データの情報を調べて、自分がダウンロードものを選べます。

ほしいデータのDP.IDを指定してretrieve_dataというメソッドを使ったらデータをダウンロードできます。

with_calibというキーワードを'raw'にしたらダークやフラットなど補正するための画像データも一緒にダウンロードされます。

普通は、同じ日のデータは同じセットのデータです。with_calibを'raw'にしたら同じセットのデータは全部一緒にダウンロードできます。

今回は2015年09月28日のデータをダウンロードしてみます。

ダウンロードするにはまずloginメソッドでログインして、パスワードを入力したらダウンロードは始められます。

eso.login('ユーザー名', store_password=True)
data = eso.retrieve_data(['SPHER.2015-09-28T02:23:34.953'],destination='保存するフォルダ',with_calib='raw')

ただ、普通はこうするとエラーが出る場合が多いです。その場合、直接サイトにログインしてダウンロードするしかないのです。

普通はretrieve_dataを使ったとたんESOサイトからメールは送ってきます。そのメールの中のリンクをクリックしてログインしたらデータをダウンロードすることができます。

それとも、このリンクでDP.IDを入力してリクエストをしてダウンロードすることもできます。
http://archive.eso.org/cms/eso-data/eso-data-direct-retrieval.html

ダウンロードされたデータを確認

ダウンロードされたファイルはいっぱいあって、その中には対象自体のデータとダークやフラットなど補正するためのデータが含まれています。

どのファイルが何のデータなのか、ヘッダの中の『hierarch eso dpr type』というキーを調べればわかります。

ダウンロードされたファイルは圧縮されたzipファイルなので、まずは解凍する必要があります。解凍したら.fitsファイルになります。pythonでは.fitsファイルはastropyを使って開けられます。

全てのファイルを開けてヘッダを調べてみます。

from astropy.io import fits
from glob import glob
import pandas as pd
import os

fd = 'データの保存されたフォルダ'
dpr_type = pd.Series()
for f in sorted(glob(os.path.join(fd,'SPH*.fits'))):
    dpr_type[f.split('/')[-1]] = fits.getheader(f)['hierarch eso dpr type']

print(dpr_type)

結果

SPHER.2015-09-28T02:23:34.953.fits             OBJECT
SPHER.2015-09-28T02:28:24.730.fits             OBJECT
SPHER.2015-09-28T02:33:14.199.fits             OBJECT
SPHER.2015-09-28T02:38:03.790.fits             OBJECT
SPHER.2015-09-28T02:42:53.412.fits             OBJECT
SPHER.2015-09-28T02:47:43.182.fits             OBJECT
SPHER.2015-09-28T02:52:32.769.fits             OBJECT
SPHER.2015-09-28T02:57:22.273.fits             OBJECT
SPHER.2015-09-28T03:02:11.572.fits             OBJECT
SPHER.2015-09-28T03:07:01.002.fits             OBJECT
SPHER.2015-09-28T03:11:50.616.fits             OBJECT
SPHER.2015-09-28T03:16:39.785.fits             OBJECT
SPHER.2015-09-28T03:21:29.073.fits             OBJECT
SPHER.2015-09-28T03:26:18.236.fits             OBJECT
SPHER.2015-09-28T03:31:07.473.fits             OBJECT
SPHER.2015-09-28T03:35:56.912.fits             OBJECT
SPHER.2015-09-28T03:41:18.749.fits      OBJECT,CENTER
SPHER.2015-09-28T03:42:13.969.fits        OBJECT,FLUX
SPHER.2015-09-28T03:44:20.557.fits                SKY
SPHER.2015-09-28T03:45:06.576.fits                SKY
SPHER.2015-09-28T03:45:52.561.fits                SKY
SPHER.2015-09-28T10:36:15.112.fits               DARK
SPHER.2015-09-28T11:59:24.391.fits               DARK
SPHER.2015-09-28T12:57:39.906.fits    DARK,BACKGROUND
SPHER.2015-09-28T13:12:49.106.fits    DARK,BACKGROUND
SPHER.2015-09-28T14:01:51.749.fits    DARK,BACKGROUND
SPHER.2015-09-28T14:17:13.269.fits          FLAT,LAMP
SPHER.2015-09-28T14:17:44.993.fits          FLAT,LAMP
SPHER.2015-09-28T14:21:25.267.fits          FLAT,LAMP
SPHER.2015-09-28T14:22:38.309.fits          FLAT,LAMP
SPHER.2015-09-28T14:26:01.484.fits          FLAT,LAMP
SPHER.2015-09-28T14:37:57.083.fits          FLAT,LAMP
SPHER.2015-09-28T14:38:19.050.fits          FLAT,LAMP
SPHER.2015-09-28T14:41:48.845.fits          FLAT,LAMP
SPHER.2015-09-28T14:42:30.701.fits          FLAT,LAMP
SPHER.2015-09-28T14:46:24.616.fits          FLAT,LAMP
SPHER.2015-09-28T15:04:51.724.fits       LAMP,DISTORT
dtype: object

dpr typeがOBJECTであるファイルはほしい対象自体の生データです。その他は全部補正用のデータです。

matplotlibを使ってそれぞれのファイルを見てみます。

描く時は重要な部分を強調するためにastropyのZScaleIntervalでスケールします。

zscaleについてはこの記事で書いてあります https://qiita.com/phyblas/items/87667c250b29e195f7fb

import matplotlib as mpl
import matplotlib.pyplot as plt
from astropy.visualization import ZScaleInterval,ImageNormalize

for dt in set(dpr_type): # タイプ毎一枚だけ取り出す
    obj = fits.getdata(os.path.join(fd,dpr_type.index[dpr_type==dt][0])) # 画像読み込み
    plt.figure(figsize=[7,7]) # 画像

    # logscale
    plt.subplot(211,facecolor='#000000')
    lognorm = mpl.colors.LogNorm()
    plt.imshow(obj[0],origin='bottom',norm=lognorm,cmap='plasma')
    plt.colorbar(pad=0.01,aspect=50)

    # zscale
    plt.subplot(212,facecolor='#000000')
    norm = ImageNormalize(obj[0],interval=ZScaleInterval())
    plt.imshow(obj[0],origin='bottom',norm=norm,cmap='plasma')
    plt.colorbar(pad=0.01)

    plt.suptitle(dt)
    plt.tight_layout(0)
    plt.savefig(dt.lower()+'.png') # 保存
    plt.close()

IRDISのデータは2つの波長が一緒に含まれているので、画像は左と右2つの部分に分けられます。左と右の部分は違う波長の写真です。

データ解析

データのファイルが揃ったら次はvlt-sphereモジュールを使ってデータ解析を行います。つまり、生データを補正するのです。

vlt-sphereモジュールを使う

まずは一つフォルダを作成して、その中でrawというフォルダを作成して生データと補正用のデータを全部そのrawのフォルダに置きます。

例えばホームにHR8799_20150928というフォルダを作成したらデータは

~/HR8799_20150928/raw/

に置くのです。

そしてpythonを使ってこのコードを実行するだけでデータ解析が始まります。

import sphere.IRDIS as IRDIS
path = '~/HR8799_20150928/'
reduction = IRDIS.ImagingReduction(path)
reduction.full_reduction()

注意：pathはホームの代わりに~を使うこともできますが、相対パスを使うことができません。

しばらく待って、間違いがなければデータ解析は完了します。

その後productsというフォルダに解析完了の画像データが置かれます。

~/HR8799_20150928/products/

解析完了の画像データを確認

解析したあとのデータはscience_cube.fitsというファイルにあります。

データは4次元

第1次元は違う波長を示す
第2次元は違う時間を示す
第3-4次元は画像の空間

2つの波長はwavelength.fitsに書いてあります。単位はnm（ナノメートル）

例として最初と最後のフレームの画像を見てみます。2つの波長は別々に表示します。

import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.visualization import ZScaleInterval,ImageNormalize

sc = fits.getdata('HR8799_20150928/products/science_cube.fits') # 画像
wl = fits.getdata('HR8799_20150928/products/wavelength.fits')/1000 # 波長
print(sc.shape) # (2, 512, 800, 800) 4次元データ

plt.figure(figsize=[7,6])
for i in [0,1]: # 最初と最後のフレームを別々
    # 第1波長
    plt.subplot(221+i*2)
    norm = ImageNormalize(sc[0,-i],interval=ZScaleInterval())
    plt.imshow(sc[0,-i],norm=norm,origin='bottom',cmap='Blues_r')
    plt.colorbar(pad=0.01)
    if(i==0):
        plt.title(r'$\lambda=%s\mu m$'%wl[0])

    # 第2波長
    plt.subplot(222+i*2)
    norm = ImageNormalize(sc[1,-i],interval=ZScaleInterval())
    plt.imshow(sc[1,-i],norm=norm,origin='bottom',cmap='afmhot')
    plt.colorbar(pad=0.01)
    if(i==0):
        plt.title(r'$\lambda=%s\mu m$'%wl[1])

plt.tight_layout()
plt.show()

ここで違う波長だとわかりやすいように、短い方の波長を青に、長い方の波長を赤にするが、別に特に必要っいうわけではなく、グレースケールを使ってもいいです。

データから見ると512フレーム（違う時間で撮られた写真）があります。

このらの画像はノイズがいっぱいあるので、惑星を見つけることができなくて、全て合わせてノイズを消された一つの画像にする必要があります。

惑星が見えるように、次はこの数百枚の画像を合成して一つの画像にします。

角度差分撮像（ADI）で画像合成

VLT-SPHEREでは角度差分撮像（angular differential imaging）、通称ADI（Marois et al. 2006）という方法が採用されています。この方法によって「PSF」と「スペックルノイズ」というものを消してはっきりと惑星が見えるようにするのです。

スペックルノイズとは、微妙に歪んだ波面の干渉によって発生する粒状であるノイズです。地上望遠鏡は大気の影響があるためスペックルノイズは特に激しいです。

それと、望遠鏡の中の回折によって、星の見かけの姿は単なる小さな点ではなく、その点から広がった明るい領域になります。その姿を説明する関数はPSF（point spread function = 点拡がり関数）と呼ばれます。惑星は主星に近いしずっと暗いため、主星のPSFによって埋もれやすくて観測しにくいです。

PSFの形は大体エアリーディスクですが、スペックルノイズなどのノイズの影響も受けて単純なエアリーディスクではなく、姿は複雑になります。

このようなノイズは光学系に由来する準静的なノイズであり、視野に固定されて写るものなので、望遠鏡の視野が回転しても同じような形をします。

近年建設された天文台は経緯台式架台を採用することが多いです。普段経緯台式架台は対象の星を眺める時に視野はその星の周りに回転していくのです。
もしその星の周りに回転する惑星が存在したら、その惑星の位置は時間によって変わります。
その同時に主星のPSFとスペックルは視野に対して固定されるので、そのノイズと惑星の姿を区別することができます。

これについて、アニメーションでもっと詳しい説明はこの記事で https://qiita.com/phyblas/items/d3fdca7c42b798e65b89

ADIという方法は、何時間も連続で何枚の写真を撮って、それらの写真からPSFのパターンを探して、そして写真からPSFをひきます。

一番簡単な方法は全ての写真の中央値を取って、その中央値でひき算するのです。

各枚の写真の中の惑星の位置はぐるぐる廻るので、中央値に影響はないはずです。なので、中央値でひいたらPSFの影響は消えて、惑星が残ってはっきり見えます。

これについてもっと詳しくはこれらの論文を参考に
http://www.asj.or.jp/geppou/archive_open/2012_105_03/105_148.pdf
https://www.wakusei.jp/book/pp/2013/2013-4/2013-4-255.pdf
http://www.astro-wakate.org/ss2015/web/file/shuroku/star_c15.pdf
https://www.jstage.jst.go.jp/article/yuseijin/22/4/22_KJ00008993174/_pdf/-char/ja

その他に、SPHEREのような装置ではADIと共に波長差分撮像（wavelength differential imagingかspectral differential imaging、SDI）や偏光差分撮像（polarimetric differential imaging，PDI）という方法が使えます。つまり、偏光や波長の違いも加えて分析することです。

視差角

対象の天体を中心として天頂の方向と北の方向にそれぞれ線を描けば、その線の間に成す角度は視差角（parallactic angle）と呼びます。

何時間も撮影し続けると視差角は少しずつ変わります。

視差角によって写真を回転する必要があるので、ADIを使う時には視差角を知る必要があります。

vlt-sphereモジュールで解析されたデータの視差角はproductsフォルダ内のscience_parang.fitsというファイルに書いてあります。

このファイルを開けてグラフを描いてみます。

parang = fits.getdata('HR8799_20150928/products/science_parang.fits')
plt.plot(parang,'#BB8799')
plt.ylabel('parang (°)')
plt.show()

ADIを実装する

次はADIを実装するコードです。流れは以下の通り

画像と視差角を収めている.fitsファイルを開ける
全ての画像の中央値を取って、この値で全ての画像からひく
画像を視差角によって回転させる
回転したあとの画像の中央値を取る
結果の画像を描く

import matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits
from astropy.visualization import ZScaleInterval,ImageNormalize
from skimage.transform import rotate

path_sc = 'HR8799_20150928/products/science_cube.fits'
path_parang = 'HR8799_20150928/products/science_parang.fits'
path_wl = 'HR8799_20150928/products/wavelength.fits'
sc = fits.getdata(path_sc) # 解析完了の画像
parang = fits.getdata(path_parang) # 視差角
wl = fits.getdata(path_wl) # 波長
sc = sc-np.median(sc,1)[:,None,:,:] # 中央値をひく

# 回転させる
for j in [0,1]: # 2つの違う波長を別々
    for i,parang_i in enumerate(parang): # 違うフレームの画像
        sc[j,i] = rotate(sc[j,i],-parang_i,cval=np.nan)

res = np.median(sc,1) # 中央値
res = res[:,200:-200,200:-200]

# 画像を表示してみる
## 第1波長
plt.figure(figsize=[7,6])
norm = ImageNormalize(res[0],interval=ZScaleInterval())
plt.imshow(res[0],norm=norm,origin='bottom',cmap='Blues_r')
plt.colorbar(pad=0.01)
plt.title(r'$\lambda=%s\mu m$'%wl[0]) # 波長を記す
plt.tight_layout()

## 第2波長
plt.figure(figsize=[7,6])
norm = ImageNormalize(res[1],interval=ZScaleInterval())
plt.imshow(res[1],norm=norm,origin='bottom',cmap='afmhot')
plt.colorbar(pad=0.01)
plt.title(r'$\lambda=%s\mu m$'%wl[1])
plt.tight_layout()

plt.show()

結果

波長によって回折が違って長い波長の方はPSFの範囲が広いし、スペックルも大きくなるので、違う波長で起こったノイズは別々ですが、それに対し、天体の位置は波長には関係ないので、惑星である点はどんな波長でもちゃんと同じ位置に写るはずです。

なので、2つの波長を合わせて調べたら惑星の点はもっとはっきりと強調されます。

2つの波長を合わせてみたらこうなります。

norm = ImageNormalize(res,interval=ZScaleInterval())
plt.imsave('HR8799.jpg',norm(np.stack((res[1],(res[0]+res[1])/2,res[0]),2))[::-1])

画像をよく見てみたら3つの白い点が見つけます。それはHR8799b、HR8799c、HR8799d

実はもっと中心に近いところにはもう一つHR8799eという惑星が存在するはずなのですが、主星から一番近いため見えにくくて無理です。

発見の歴史によっても、最初は2008年からHR8799bとHR8799cとHR8799d、3つだけは先に発見されて、そのあと一番見えにくいHR8799eは2010年になってから発見を確認されたそうです。

ADIの派生

ADIを基礎として発展された手法は色々あります

例えば

環状ADI（annular mode ADI）
PCA（principal component analysis）とKLIP（Karhunen-Loève image projection） Amara & Sascha (2012), Soummer et al. (2012), Wang et al. (2015)
NMF（non-negative matrix factorization） Ren et al. (2018)
LOCI（locally optinized combination of images） Lafrenière et al. (2007)
LLSG（local low-rank plus sparse plus Gaussian noise decomposition） Gonzalez et al. (2016)
IRS (image rotation and subtraction) Ren et al. (2012), Dou et al. (2015)
ANDROMEDA（angular differential optimal exoplanet detection algorithm） Mugnier et al. (2009), Cantalloube et al. (2015)
などなど

これらの方法を使えば本来のADIよりもっと綺麗に合成できます。

これに対して本来のADIは古典ADI（classical ADI）と呼ばれることもあります。

方法は更に複雑になるので、ここでは説明しませんが、vip_hciというモジュールを使って実装します。

vip_hciによるADI

vip_hciはADIを実装するためのpythonモジュールです。これを使ったら簡単に色んなADIから派生した方法が使えます。

vip_hciを使って本来のADIを含めて5つの方法で画像を合成します

import vip_hci as vip
import matplotlib.pyplot as plt
from astropy.visualization import ZScaleInterval,ImageNormalize
import os

fd = 'HR8799_20150928/products'
sc = vip.fits.open_fits(os.path.join(fd,'science_cube.fits'))
parang = vip.fits.open_fits(os.path.join(fd,'science_parang.fits'))

res = [{},{}]
for j in [0,1]:
    res[j]['classical ADI'] = vip.medsub.median_sub(sc[j],parang)[200:-200,200:-200]
    res[j]['annular ADI'] = vip.medsub.median_sub(sc[j],parang,mode='annular')[200:-200,200:-200]
    res[j]['NMF'] = vip.nmf.nmf(sc[j],parang,ncomp=4)[200:-200,200:-200]
    res[j]['PCA'] = vip.pca.pca(sc[j],parang,ncomp=4)[200:-200,200:-200]
    res[j]['annular PCA'] = vip.pca.pca_annular(sc[j],parang,ncomp=4)[200:-200,200:-200]
    res[j]['LLSG'] = vip.llsg.llsg(sc[j],parang,4)[200:-200,200:-200]
    res[j]['LOCI'] = vip.leastsq.leastsq.xloci(sc[j],parang)[200:-200,200:-200] # LOCIは凄く時間がかかるので注意

for meth in res[0]:
    plt.figure(figsize=[5,8])

    plt.subplot(211)
    norm = ImageNormalize(res[0][meth],interval=ZScaleInterval())
    plt.imshow(res[0][meth],norm=norm,origin='bottom',cmap='PuBu_r')
    plt.colorbar(pad=0.01)

    plt.subplot(212)
    norm = ImageNormalize(res[1][meth],interval=ZScaleInterval())
    plt.imshow(res[1][meth],norm=norm,origin='bottom',cmap='afmhot')
    plt.colorbar(pad=0.01)

    plt.tight_layout()
    plt.savefig('HR8799 '+meth+'.jpg')
    plt.close()

古典ADI

環状ADI

NMF

PCA

環状PCA

LLSG

LOCI

こんな風に、違う方法で合成すると結果は違います。

他の惑星系

ここでHR8799の惑星を例として挙げましたが、その他にも色々直接撮像法で目撃することができる系外惑星があります。

VLT-SPHEREは色々な惑星系の写真を撮っています。その星達のデータを調べてみてもいいのです。

すでに発見された惑星を見たいなら、astroqueryのNasaExoplanetArchiveモジュールを使って、惑星のデータを取得することができます。

これについてはすでにこの記事で説明しました
https://qiita.com/phyblas/items/ec29c49302c7b0a60905

今まで直接撮像法で発見された惑星のデータを見てみましょう。

from astroquery.nasa_exoplanet_archive import NasaExoplanetArchive
tb = NasaExoplanetArchive.get_confirmed_planets_table(cache=0)
df = pd.DataFrame({c: tb[c] for c in tb.columns})
df.set_index('pl_name',inplace=True)
df[df['pl_discmethod']=='Imaging']
print(df[df['pl_discmethod']=='Imaging']['pl_hostname'].sort_values())

今のところでは47の系外惑星の名前が出ます。

pl_name
1RXS J160929.1-210524 b                1RXS J160929.1-210524
2MASS J01225093-2439505 b            2MASS J01225093-2439505
2MASS J02192210-3925225 b            2MASS J02192210-3925225
2MASS J04414489+2301513 b            2MASS J04414489+2301513
2MASS J12073346-3932539 b            2MASS J12073346-3932539
2MASS J21402931+1625183 A b        2MASS J21402931+1625183 A
2MASS J22362452+4751425 b            2MASS J22362452+4751425
51 Eri b                                              51 Eri
AB Pic b                                              AB Pic
CFBDSIR J145829+101343 b              CFBDSIR J145829+101343
CHXR 73 b                                            CHXR 73
CT Cha b                                              CT Cha
DH Tau b                                              DH Tau
FU Tau b                                              FU Tau
Fomalhaut b                                        Fomalhaut
GJ 504 b                                              GJ 504
GQ Lup b                                              GQ Lup
GSC 06214-00210 b                            GSC 06214-00210
GU Psc b                                              GU Psc
HD 100546 b                                        HD 100546
HD 106906 b                                        HD 106906
HD 203030 b                                        HD 203030
HD 95086 b                                          HD 95086
HIP 65426 b                                        HIP 65426
HIP 78530 b                                        HIP 78530
HIP 79098 AB b                                  HIP 79098 AB
HN Peg b                                              HN Peg
HR 2562 b                                            HR 2562
HR 8799 c                                            HR 8799
HR 8799 b                                            HR 8799
HR 8799 e                                            HR 8799
HR 8799 d                                            HR 8799
LkCa 15 c                                            LkCa 15
LkCa 15 b                                            LkCa 15
Oph 11 b                                              Oph 11
PDS 70 b                                              PDS 70
PDS 70 c                                              PDS 70
ROXs 12 b                                            ROXs 12
ROXs 42 B b                                        ROXs 42 B
Ross 458 c                                          Ross 458
SR 12 AB c                                          SR 12 AB
USco CTIO 108 b                                USco CTIO 108
VHS J125601.92-125723.9 b            VHS J125601.92-125723.9
WD 0806-661 b                                    WD 0806-661
WISEP J121756.91+162640.2 A b    WISEP J121756.91+162640.2 A
bet Pic b                                            bet Pic
kap And b                                            kap And
Name: pl_hostname, dtype: object

これらの星の名前で調べればいいです。

HIP65426

もう一つの例として、上述と同じような方法をHIP65426という星に実行してみます。この星には一つ系外惑星が発見されています。

結果は以下