Phoebe2を使ってみる（データ入力と表示編）

　　　　　　　　　　　　　　　　　　　　　　　　　　　 　　　　　　　　　　　　　　　　　(2021.11.06)

はじめに

　自分で観測した食連星系の光度曲線を解析して連星系の物理パラメータを決定するというのは、醍醐味のあることです。これまで必要なときは、Phoebe1（Version0.3）を使ってこの作業を行っていましたが、何分Phoebe1は、かなり良いパラメーターを与えないと簡単にプログラムが異常終了してしまうので、設定ファイルを頻繁に保存しておかないと、振り出しに戻り、困り物のソフトでもありました。そこでかなり安定して高機能になったPhoebe2に早く移行したかったのですが、ユーザーインターフェースがこれまでとはガラリと変わり、さっぱり扱い方がわかりませんでした。

　悪戦苦闘して、サイト上のドキュメント（主にTutorial）を読みながらPCを操作して来ましたが、途中で諦めようかと思うくらい理解が進まず、これだけのことを行えるようになるまでに1ヶ月半もかかりました。しかし、やっと観測データをモデル計算のデーターセットに与える方法がわかりましたので、その中間報告を書きます。

　 　そのポイントは、Phoebe２の仕組みを理解することでした。その仕組みは、モデル計算に使うデーターセットに、観測データのフラックスと時刻を与えるということでした。理解が進まなかったのは、私が間違ってその仕組を予想していたという側面が大きいです。間違って予想していた仕組みとは、データセットには、観測のデータセットと計算モデルのデータセットの2つがあって、その2つのデータセットを比べて計算モデルを改良して観測データに合わせる、というものでした。

　 　以下に書くのは、観測データをPhoebe2に読み込ませ、連星系の大雑把な初期値から計算した光度曲線と同時にグラフ表示して見比べるところまでです（モデル計算と観測は全く合っていません）。連星系のパラメータを改良して最適解を見つけるにはどうしたらよいかということは、これからまだ手探りしながら調べていきます。

参考にしたオンライン資料へのリンク

　Docments

　Tutorials

　PHOEBE Workshop(June 2021) Virtual Workshopスライドの動画も見れてわかりやすい

Phoebe2の原理についての論文

　Physics of Eclipsing Binaries. II. Toward the Increased Model Fidelity

　Physics of Eclipsing Binaries. III. Spin–Orbit Misalignment

　Physics of Eclipsing Binaries. IV. The Impact of Interstellar Extinction on the Light Curves of Eclipsing Binaries

　Physics of Eclipsing Binaries. V. General Framework for Solving the Inverse Problem　その日本語訳　ぜひ読むことをお薦めします

　　以下は、Miniconda（またはAnaconda）でPython3,Phoebe2,その他必要な依存関係パッケージをインストールし、jupyter notebook上でphoebe2が動く環境での説明です。

　　作業ディレクトリで 　　$ conda activate phoebe 　　としておく。

　　持っている測光観測データから次のことを準備する。

　（１）極大等級として位相0.25か0.75付近の等級m_maxを調べる. 　 　（例）m_max=1.245

　（２）観測データから、時刻、極大に対する相対光度（相対フラックス）、誤差の3項目が並んだデータファイル（textファイル）を作る。

　　　Phoebe2が膨大な計算に要する時間を短くしたいので、観測データ数はできれば100個に押さえたい。そこで、観測データを位相0.00から0.99までの100個のビンに入れて、その平均値と標準偏差を求める＜付録B参照＞。

（例）ファイルのフォーマット
     #time    f_rel  sigma
     0.00000  0.854  0.008
     0.00332　0.872  0.006
      ‥‥

　　これを.例えばファイル名"AB_And21,dat"としてpath名"/mnt/c/ab_and/"に保存するとする。（以下、W UMa型=接触型食連星の例）

データの準備ができたら、いよいよPhoebe2を動かしてみましょう。 まず、次の4行は、Phoebe2を動かすためのいわゆる「おまじない」と思って黙って実行します。第4行目は、記録に残すレベルを指定するコマンドで、'WARNING'だと警告レベルから表示されます。もう少し情報を減らしたい時は、'ERROR'にします。初めての人は'WARNING'を指定して、丁寧にそのメッセージを読むと参考になります。英語のメッセージでは意味がわからんという人は、コピーしてGoogleなりDeepLなりで翻訳してもらうと、それなりに意味が取れて、結構よいヒントになります。

import phoebe
from phoebe import u # units
import numpy as np
logger = phoebe.logger(clevel='WARNING')

次に、連星系のオブジェクトのインスタンスを作ります。（オブジェクトもインスタンスも、オブジェクト指向言語の技術用語です。天文学では、オブジェクトと言えば天体のことですが、混同しないように文脈で捉えてください）

b = phoebe.default_binary(contact_binary=True)

（bがインスタンスで、phoebe.default_binary()がオブジェクトのクラスです‥‥プログラミング用語） これで、１つの連星系が計算機のメモリ上に作られました。この時点では、その各種パラメータは適当な仮のデフォルト値が代入されています。 ここから手順を追ってphoebe2に指示を出していくわけですが、その大まかな手順は＜付録C　Phoebe2を動かす手順の概要＞を参照。

観測データファイルを読み込みます。

times, fluxes, sigmas = np.loadtxt('/mnt/c/ab_and/AB_And21,dat',comments='#',unpack=True)
b.add_dataset('lc',times=times,fluxes=fluxes,sigmas=sigmas,
            dataset='lc01',overwrite=True)
b.set_value('pblum_mode@lc01','dataset-scaled')  #モデル計算するフラックスを観測データの相対フラックスにスケールを合わせる
b['period@binary']=0.3318912

　先行研究の論文などから既知のパラメータがあれば入力しておく（注）。

b['passband@lc01'] = 'TESS:T'  #バンドパス名
b['incl@binary']=85.45
b['teff@primary']=5890
b['teff@secondary']=5873
b['gravb_bol@primary']=0.32
b['gravb_bol@secondary']=0.32
b['l3_mode@lc01']='fraction'  #3rd light
b['l3_frac@lc01']=0.28

（注）入力できるパラメータは、どれでも自由にできるとは限らず、初期設定のままで自由に入力できるのは次の表の左端の量だけである。右端のパラメータは入力できないが、初期設定の制約を外して右端にあるパラメータと左端のそれを入れ替える事(Phoebe2では"flipping constraints"と呼ぶ)もできる。 TutorialsのConstraintsを参照のこと

初期設定のパラメータ

b.run_compute(irrad_method='none', model='tess')
afig, mplfig = b.plot(x='phases',ylabel='flux_rel', show=True)

ここまでで、測光観測データと大雑把な初期モデルによる光度曲線を見ることができた。次はここから出発して、食連星系のパラメータを改良して満足のいくモデルを構築することになる。

　　まずOSとしては、Phoebe2はLinux, MacOS,Windowsの３つのOS上で動くそうですが、Windowsは現時点（2021年秋）でまだ動作が不完全な可能性があるので、前2者がよいでしょう。普段はWindows10/11を使っている人は、WSL2でLinux（Ubuntu20.04LTS）を入れるのがお勧めです。私もWindows10+WSL2+Ubuntu20.04環境です。Windowsのホルダー中のファイルがシームレスに使えるので便利です。

　　 　　MinicondaでPython3,Phoebe2,その他必要な依存関係パッケージをインストルするのが手っ取り早いです。 　http://phoebe-project.org/install/latest/linux/auto

なお、windows10へのWSL2、Ubuntu20.04のインストールについては、Web検索するとたくさん出てきます。なるべく新しいもので丁寧な記述のサイトを選びましょう。

例

https://www.kkaneko.jp/tools/wsl/wsl2.html

まずCondaの導入方法:

以下の説明は、磯貝 瑞希さんの「PyRAF ミニ講習会」資料(2020 MAR 12 国立天文台・天文データセンター）を

2つのdistribution(Miniconda/Anaconda)があり、それぞれに Python 2.7系 (=Miniconda2/Anaconda2) と Python 3.x系 (=Miniconda3/Anaconda3) が提供されている :

- Miniconda: 必要最小限のConda 管理環境を提供

- Anaconda: 完全なConda 管理環境 + 数百の有用なツール・ライブラリを(デフォルトで)提供

　軽量な Miniconda(=Miniconda3)を採用

ブラウザーを起動し、本家サイトのdownloadページ

URL: https://docs.conda.io/en/latest/miniconda.html

より、 Linux installers の Python 3.x, 「Miniconda Linux 64-bit」 を選択、表示されるダイアログで「Save File」を選択し、 OKを押す

ダウンロードしたファイルのハッシュ値計算と確認: $ sha256sum /home/student/Downloads/Miniconda3-latest-Linux-x86_64.sh → 本家サイトのハッシュ値と一致すればOK.

ダウンロードしたファイルを bash で実行: (=ユーザ権限でのインストール) $ bash /home/student/Downloads/Miniconda3-latest-Linux-x86_64.sh

→ Enterを押すとライセンス条項が表示される。最後まで進んで、「yes」を入力

→インストール先[/home/student/miniconda3」はそのままとし、 Enterを押す

→ Miniconda3 の初期化の実施: yes

→ ~/miniconda3 ディレクトリが作成される。ディレクトリのサイズは約300MB

確認:

$ source ~/.bashrc

→ condaのbase環境がactivateされる

→ コマンドプロンプトの前に、「 (base) 」 と表示されるようになる デフォルトの設定では、次回以降、ログイン時に自動起動

参考: 自動起動設定の解除

$ conda config –set auto_activate_base false
$ which conda

→ ~/miniconda3/bin/conda

Phoebe2 のインストール (ここでは環境名: phoebe とする )

$ conda create -n phoebe python=3.8
$ conda activate phoebe

→ Proceed([y]/n)?: [Enter] → 問題なければ5-6分程度で終了 (環境=通信帯域幅に強く依存)

$ conda deactivate
$ conda install -c conda-forge notebook
$ conda install -c conda-forge nb_conda_kernels
$ sudo add-apt-repository ppa:ubuntu-toolchain-r/test
$ sudo apt-get install g++-5 gcc-5 libstdc++-5-dev
$ export CXX=g++-5
$ sudo apt-get install python3-dev

一緒にnotebookもインストール

$ conda install -c conda-forge notebook
$ conda install -c conda-forge nb_conda_kernels

PHOEBE2環境のインストール

$ conda create -n phoebe python=3.8
$ conda activate phoebe
$ conda install -c anaconda ipykernel
$ python -m ipykernel install --user --name phoebe_env --display-name "PHOEBE env"
$ pip install numpy phoebe

Phoebe2環境に入るには

$ conda activate phoebe

Phoebe2環境から抜けるには

$ conda deactivate

　最終的にデータファイルに入れる時刻は、各位相0.00から0.99（か0.995）に対応した実時間（単位はday）とする。すなわち時刻=位相×周期で与える。 　

　相対光度（相対フラックス）は、極大等級に対する相対値で、まずその時刻における値f_relを次式で計算する。その時刻の等級をm、極大等級をm_maxとすると、

 f_rel=10^-0.4(m-m_max) .

ここで使う等級は、一貫していれば、絶対値でも相対値でも良い

位相0.00から0.99までの各位相ごとにf_rlの平均値と標準偏差（=誤差とする）を求める。ファイルには位相でなく時刻に直した値を入力する。

（例）
     #time    f_rel  sigma
     0.00000  0.854  0.008
     0.00332　0.872  0.006
      ‥‥

以下は、概要ですが、詳しくはTutorialsのGeneral Concepts: The PHOEBE Bundleを読んでください。

０．「おまじない」を実行

１．計算機の中に「連星系」を作る

b = phoebe.default_binary(contact_binary=True)

２．その「連星系」のデータセットに、観測データを与える（観測データの代わりに、Phoebe2にデータセットを与えさせることもできます。Tutorialsでは後者で話を進めています）

b.add_dataset()

３．データセット中のパラメータで既知の値があれば、できるだけ多く与えておく。（ただし制約を受けるパラメータに注意。上記本編３．の（注）を参照）

b.set_value(qualifier='teff', component='primary', value=5890)

または、次の書き方でも同じ

b['teff@primary']=5890

４．連星系モデルを実行して、ここまでに与えたデータセットから定義される観測変数(observables)から成る合成モデルを作る。

b.run_compute()

　これを実行すると、生まれたての連星系モデルにアクセスが可能になり任意のパラメータの値を表示させたり、代入したりできます。

print(b)

　グラフにプロットすることもできます。

afig, mplfig = b.plot(show=True)

上記本編で扱っているのは、ここまでです。以下は連星系のパラメータを改良して観測とうまく合わせる（逆問題を解く）までの手順の概要です。

５．逆問題を解く

　（１）パラメータの可能な範囲distributionsの追加

　（２）観測変数の推定を行う（estimators）

　（３）最適化（optimizers）、または最適化のためのサンプラー（samplers）を実行する 　

１．コマンドプロンプトからwsl2を起動しておく

> wsl

　　"wsl2"の２は不要というより、2をつけると起動しない。

　　私の場合は、Ubuntu20.04 LTSをインストールした。

　　WSL2の起動は、Powershellからでもよいが、軽いのでわたしはcmdを使っている.

２．windows10からX-termソフトを立ち上げ、そこからphoebe2の環境に入る

  $ conda activate phoebe

　　X-termが起動したら、もとのコマンドプロンプトはexitで終了させても、wsl2環境は引き続きX-term上から利用できます。

　　私の場合はMobaXterm、軽くて、いろいろと便利な機能があります。計算結果待ち時間にソリティアゲームもできるし。

３．使用する作業ディレクトリへ移り、notebookを起動する

  $ cd /mnt/c/etc/TZ_Boo/
  $ jupyter notebook

　firefoxが立ち上がりnotebook画面になるまでしばらく待つ

　なお、Phoebe2は必ずしもnotebook環境で作業する必要はなく、コマンドラインからpython scriptを実行する形でも動きますし、その方が実行時間は２割ほど早いようです。しかし、慣れないうちは、notebook環境上でPhoebe2に触った方が扱いやすいし、自動で保存され記録が残るなど、何かと便利なように軟弱な私は思います。

４．そして存分にPhoebe2を使用する。 　慣れないうちは大変です。覚悟の上でしっかりDocments類に目を通して（大いにDeepL翻訳を利用して）、恐れずに触りまくって慣れましょう。

５．終了は、X-term上で、"ctrl-c"。５秒以内に"y"を押すとnotebook 環境が閉じます。続いてX-term自体もexitで終了させます。

OTO-Ohshima Tamashima Observatory-トップへ戻る