地学 II 関連の計算機実習: HR図をつくる

恒星のHR図をつくる。手作業では困難な多数のデータを、Python で図示してみる。

【[2021年度の他大学の授業のために R でつくった教材] をもとにしている。】

つぎのウェブページを参考にした。

データ: ヒッパルコス星表

ヒッパルコスは古代ギリシャの天文学者だが、 その名まえが、現代のヨーロッパ宇宙機構(ESA)による恒星位置測定用衛星の名まえに採用された。 ここでいう「ヒッパルコス星表」はその衛星を使った研究プロジェクトの成果である。

このデータは公開されていて、日本の国立天文台の次のところに(も)ある。

この hip_main.dat.gz というファイルは、gzip というソフトウェアで圧縮されていて、 GNUソフトウェアにふくまれる (Linuxではかならずある) gunzip で展開すると、 hip_main.dat というテキストファイルになる。

これを、Okamoto さんと Ueno さんにならって、awk [オーク] というプログラム言語で加工し、 Ueno さんが抽出したのとだいたい同じ、次の項目を抽出した。 【ただし、天球上の星のみかけの位置をあらわす 赤経・赤緯は、 Ueno さんは、天文観測の慣例にしたがった表現 (赤経は時間の時・分・秒、赤緯は角度の度・分・秒) によっているが、 hip_main.dat には角度の「度」単位の小数で表現された項目もあるので、わたしはそちらをえらんだ。】

スペクトル型も書きだしたほうがよいかもしれないが、今回は省略した。

hip_main.dat には、数値項目の値がはいっておらず、空白になっているところがある。 上記 2～8 のいずれかの項目が空白である行は、採用しない。

年周視差の値が負になっている行も、採用しない。

誤差の大きいデータを採用しないという観点から、 年周視差の相対誤差(「年周視差の標準誤差」を年周視差でわったもの) が 0.1 以上のものと、 色指数B-Vの誤差が 0.1 以上であるものは、採用しないことにした。 【この条件は、Okamoto さんにならったものだが、Okamoto さんは、色指数についても相対誤差による基準をつかっていた。 色指数は 0 にちかい値をとりうる量なので、相対誤差によるよりも標準誤差による基準がよいと判断した。】

その結果を CSV (コンマくぎりのテキストファイル) で書きだした。(ファイル名の「hip」のあとに下線文字「_」がある。)

1行が1つの星のデータをあらわし、20542 行ある。

実習にはつかわないが、抽出につかった awk プログラムもリンクしておく。 (Linux上で「awk -f extract.awk hip_main.dat > hip_work.csv」として実行した。)

残念ながら、星の名まえは、hip_main.dat に ふくまれていないので、hip_work.csv にもない。

つぎのサイトの「ヒッパルコス星表」のページの「恒星名データ」のところに、 3215個の星についての、ヒッパルコス星表 ID と 星座のメンバーという形の星の名 との対応表 hip_name.csv (文字コード UTF-8) がある。

つぎのサイトには、みかけの明るい星についての、星の固有名とヒッパルコス星表IDとの対照表がある。

数量の計算

絶対等級をもとめる

みかけの等級 vmag と 年周視差 prlx から計算する。標準的な公式はつぎのようなものである。

絶対等級 = みかけの等級 + 5 + 5 × log10(年周視差) ただし、年周視差の単位は角度の秒。

Python のプログラムでは次の形になる 【実際につかうプログラム中では pandas のデータフレームの要素 (たとえば df.vmag[krec]) を操作するが、 このページでは数値の演算式をしめすため、単独の変数 (たとえば vmag ) の操作にしておく。】

ここでは、年周視差が、角度のミリ秒の単位であたえられているので、値を10の3乗で割ってからつかう。

amag = vmag + 5 + 5*math.log10(prlx/1.0e3)

相当黒体温度をもとめる

色指数(B-V)から相当黒体温度を概算する。色指数が小さい (青い) ほど温度が高い。 ここでは、Wikipedia 英語版 [[Color index]] の記事に採用されている Ballesteros (2012) の式を採用する。

T = 4600 K × ( 1/(0.92×(B-V) + 1.7) + 1/(0.92×(B-V) + 0.62) )

Python のプログラムでは次の形になる。相当黒体温度 (単位 K) を temp_eff とした。

temp_eff = 4600 * ( 1/(0.92*bminusv+1.7) + 1/(0.92*bminusv+0.62) )

文献

光度 (太陽を1とした相対値) をもとめる

絶対等級 amag から、つぎの式でもとめることができる。 ただし、太陽の絶対等級 4.83 をつかう。 0.4 は、等級が5だけふえると光度が100分の1になることを10の指数でまとめたもの。

光度(相対値) = 10 (4.83 - 絶対等級) × 0.4

Python のプログラムでは次の形になる。光度(相対値)を rel_luminos とした (relative luminosity の略)。

rel_luminos = 10**((4.83-amag)*0.4)

半径 (太陽を1とした相対値) を概算する

放射輝度は相当黒体温度の4乗に比例する。光度は放射輝度に表面積をかけたものである。 表面積は半径の2乗に比例する。 以上のことから、光度 (相対値) と相当黒体温度から、半径 (相対値) を概算することができる。 ただし、太陽の相当黒体温度 5800 K をつかう。平方根を便宜上「sqrt()」と書いておく。

半径(相対値) = ((太陽の相当黒体温度/相当黒体温度)2 × sqrt(光度(相対値))

Python のプログラムでは次の形になる。半径(相対値)を rel_radius とした。

rel_radius = ((5800/temp_eff)**2) * math.sqrt(rel_luminos)

Python による プログラム例

プログラムの実行手順

教材準備のつごうで、Windows 上の Anaconda Power Shell での実行のばあいについて説明する。

発展としてできそうなこと

なにかの条件にあてはまる星だけを表示すること。 手順の紹介を用意できなかったが、Pythonプログラミングに関心のある人は試みてほしい。 Python の pandas のデータフレームの機能をつかい、条件にあてはまる行だけをとりだし、 もうひとつのデータフレームとし、それ以後はそちらをつかうのがよいだろう。

できればやりたいが、まだ用意できていないこと

表計算ソフトウェアによる作図

Excelなどの表計算ソフトウェアでも、「散布図」の機能をつかえば、H-R図をつくることができる。

実視等級、相当黒体温度などの計算を、 Ueno さんは計算を Excel の中で VBA 言語によってやっているが、 わたしは、計算式をセルに書きこむ形のほうがよいと思う。 ただし色をきめるところは、やや複雑になるので、VBA プログラムを組む必要があるかもしれない。

たとえば、hip_work.csv の内容をそのまま、ワークシートの A..H 列に読みこんだとしよう。 たぶん、つぎのようにすれば各変数の値を計算することができると思う。ただし、まだ動作確認していない。

I列で絶対等級を計算する。
セル「I1」に 「=D1 + 5 + 5*log10(E1/1.0e3)」を入れ、その内容をI列の各行にコピーする。 (log10という関数があることを前提としている。) コピーすると、参照されたセル番号「D1」「E1」は、その位置に応じて書きかわっていく。それを見こんだ計算式である。 以下同様。

J列で相当黒体温度を概算する。
セル「J1」に「=4600 * ( 1/(0.92*G1+1.7) + 1/(0.92*G1+0.62) )」を入れ、その内容をJ列の各行にコピーする。

K列で光度 (太陽を1とした相対値) を計算する。
セル「K1」に「=10^((4.83-I1)*0.4)」を入れ、その内容をK列の各行にコピーする。 (「^」(山形)が べき乗の演算子であることを前提としている。)

L列で半径 (太陽を1とした相対値) を概算する。
セル「L1」に「=((5800/J1)^2) * sqrt(K1)」を入れ、その内容をK列の各行にコピーする。 (「^」(山形)が べき乗の演算子であること、sqrtという関数があることを前提としている。)

2023-05-03
増田 耕一 (MASUDA Kooiti)
[地学 II の授業教材ページ] [立正大学むけの教材目次ページ] [教材目次ページ]