実習: HR図をつくる (その2)

[その1]も参照してほしい。

ひきつづき恒星のHR図をつくるのだが、今回は、手作業では困難な多数のデータを図示してみる。

その1の「参考教材例」でふれた

• Okamoto Yoshio (岡本 義雄)さん(高校の地学の教員)による教材ページ 「衛星データの加工とH-R図」 http://www.tennoji-h.oku.ed.jp/tennoji/yossi/Astro/ (2019-10-17 閲覧、現在はネット上にない) での GMT (Generic Mapping Tools) というソフトウェアでの作図
• Ueno Satoshi さんによるページ 「ヒッパルコス星表で太陽近傍恒星のHR図を描こう」 http://uenosato.net/hr_diagram/ での Excel VBA での作図

を参考にして、R で作図してみる。

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

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

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

• ftp://dbc.nao.ac.jp/DBC/NASAADC/catalogs/1/1239/hip_main.dat.gz ... 星表データ
• ftp://dbc.nao.ac.jp/DBC/NASAADC/catalogs/1/1239/ReadMe ... その説明文書(英語) [ReadMe の コピー]

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

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

• 列番号、変数名、説明
• 1. hid : ヒッパルコス星表での星のID。1以上の整数。いまあるデータでは最大6けた。
• 2. radeg: 赤経 (right ascension)。 単位は角度の「度」。赤道ひとまわりを360とし、0～360のあいだの値。
• 3. dedeg: 赤緯 (declination)。 単位は角度の「度」。北を正として、-90～90のあいだの値。
• 4. vmag : みかけの等級 (visible magnitude)。等級の数値。整数にかぎらない。正・負のあたいをとりうる。
• 5. prlx : 年周視差 (parallax)。単位は角度のミリ秒 (1度の60分の1の60分の1の千分の1)。
• 6. e_prlx : 年周視差の標準誤差 (error)。単位は年周視差と同じ。
• 7. bminusv : 色指数 B-V。単位は等級と同じ。
• 8. e_bminusv : B-Vの標準誤差。単位は等級と同じ。

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

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

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

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

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

• hip_work.csv (日本語文字をふくまず ASCIIだけ)

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

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

• extract.awk

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

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

• [外] Astro Commons: [ヒッパルコス星表] http://astronomy.webcrow.jp/hip/

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

• [外] 川口市立科学館: [明るい星一覧（ヒッパルコスカタログによる）] http://www.kawaguchi.science.museum/astro_room/document/brightstar50hipparcos.html

数量の計算

絶対等級をもとめる

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

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

R のプログラムでは次の形になる。 ここでは、年周視差が、角度のミリ秒の単位であたえられているので、1.0e3 がかかっている。 (あとにつけるプログラム例では、さらに、df$prlx が数値であることを明示するために as.numeric(df$prlx) としているが、 いまの hip_work.csv に対しては、そうしなくてもよい。)

amag <- df$vmag + 5 + 5*log10(df$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) )

Rのプログラムでは次の形になる。有効放射温度 (単位 K) を temp_eff とした。

temp_eff <- 4600 * ( 1/(0.92*df$bminusv+1.7) + 1/(0.92*df$bminuxv+0.62) )

文献

• F. J. Ballesteros, 2012: New insights into black bodies. EPL (Europhysics Letters), 97: 34008. doi:10.1209/0295-5075/97/34008, arXiv:1201.1809 .

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

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

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

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

rel_luminos <- 10^((4.83-amag)*0.4)

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

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

半径(相対値) = ((太陽の有効放射温度/有効放射温度)² × sqrt(光度(相対値))

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

rel_radius <- ((5800/temp_eff)^2) * sqrt(rel_luminos)

Rによる実習

プログラム例

• 注意: 日本語文字コードは、 OS によって UTF-8 か Shift_JIS をえらんでつかってほしい。
• hrhipw1.R
  • hrhipw1_utf8.R (文字コード UTF-8)
  • hrhipw1_sjis.R (文字コード Shift_JIS)
  • hip_work.csv のデータを R の「データフレーム」に読みこむところまで。 データフレームの名まえは簡単に「df」とした。
  • 【11月15日の授業中の時点でこのページに置かれたプログラムは、 読むファイル名がまちがっていました。授業後にさしかえました。】
• hrhipw2.R
  • hrhipw2_utf8.R (文字コード UTF-8)
  • hrhipw2_sjis.R (文字コード Shift_JIS)
  • 横軸に色指数(B-V)、縦軸に絶対等級をとり、H-R図を作図する。
  • どの星も同じ大きさの黒いまるで表示する。 (記号の大きさを、plotで「cex=0.2」と指定することによって、標準の0.2倍にした。)
  • [結果の図]
• hrhipw3.R
  • hrhipw3_utf8.R (文字コード UTF-8)
  • hrhipw3_sjis.R (文字コード Shift_JIS)
  • 横軸に有効放射温度、縦軸に絶対等級をとり、H-R図を作図する。 横軸のちがい以外は hrhipw2.R と同じ。
  • [結果の図]
• hrhipw4.R
  • hrhipw4_utf8.R (文字コード UTF-8)
  • hrhipw4_sjis.R (文字コード Shift_JIS)
  • 横軸に有効放射温度、縦軸に光度(太陽を1とした相対値、対数目盛り)をとり、H-R図を作成する。
  • [結果の図]
  • plot のパラメータ yaxp の対数目盛りのための機能が期待どおり動作せず、 値が 1 (=10⁰）以上のところにだけ目盛りがうたれている。 yaxp をはずせば、目盛りの位置は既定値になるが、縦軸全体に目盛りがうたれる。
• hrhipw5.R
  • hrhipw5_utf8.R (文字コード UTF-8)
  • hrhipw5_sjis.R (文字コード Shift_JIS)
  • 星の半径によって、記号の大きさを変える。 Okamotoさんにならって、図の記号の大きさを、rel_radius * 0.04 + 0.04 としてみた。 (これが適切かどうかの検討は まだできていない。)
  • Rのプログラムでは、plotにcexをベクトルであたえることができないので、 ひとつの星ごとに points をよびだすことにした。
  • [結果の図]
• hrhipw6.R
  • hrhipw6_utf8.R (文字コード UTF-8)
  • hrhipw6_sjis.R (文字コード Shift_JIS)
  • B-Vによって、記号の色を変えて、実際の星の色に近い色をあたえる。 色のRGB値は、Okamotoさんの「教材作成備忘録 2」によった。 その数値の根拠をわたしは理解できていない。 また、Okamotoさんは、RGB値を B-V の値の 0.1 きざみであたえたうえで、 GMTの機能をつかって 0.1よりこまかい中間段階の色の階調も表現しているのだが、 わたしは、B-Vをまるめて 0.1 きざみにしてつかい、それよりもこまかい階調は表現しない。
  • 背景は黒、文字や目盛りの線は白にする。
  • [結果の図]
• rgbtable.txt
  • hrhiw6.R でつかうRGB値の表。Okamotoさんの「教材作成備忘録 2」から切り出したもの。 スペースくぎり、見出し行なし。8列あるが、実際に使うのは4列で、 第1列が星のB-V値、第2, 3, 4列がそれぞれ色のR(赤)、G(緑)、B(青)の 0..255の範囲の値。
• testrgb.R
  • testrgb_utf8.R (文字コード UTF-8)
  • testrgb_sjis.R (文字コード Shift_JIS)
  • rgbtable.txt の色を確認するための簡単なプログラム。
  • [結果の図]

プログラムの実行手順

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

[R言語を使ってみる]のページと基本的に同じなので、くわしくはそちらを見てほしい。

• まず、作業用ディレクトリ(フォルダともいう)をつくっておく。
• ブラウザから、Rプログラム( ".R" でおわるファイル名) と 入力データファイル("hip_work.csv") を ダウンロードする。行き先を作業用ディレクトリにする。
• Rを起動し、上のメニューの「ファイル」から「ディレクトリの変更」で、作業用ディレクトリを指定する。
• 「ファイル」から「スクリプトを開く」でRプログラムを読みこむ。 編集画面が開いてプログラムが見えるはず。
• 編集画面を選択したうえで、上のメニューの「編集」から「全て実行」。 グラフをつくる例文プログラムならば、もうひとつ画面が開いてグラフが表示できるはず。
• グラフを画像ファイルに保存したいばあいには、 編集画面で、プログラムの「plot」よりもまえに「png("xxxx.png")」(xxxxは自分がつけたいファイル名)、 あとの「def.off()」を入れる。 (例文プログラムでは、行の頭に「##」を入れてコメント行にしてあるので、「##」をはずせばよい。) そして「全て実行」をやりなおす。作業ディレクトリに xxxx.png ができるはず。

Linux上で実行するばあいに変更するところ

• R言語処理系はエディタをふくまないので、別のエディタとくみあわせてつかう必要がある。
• 作業用ディレクトリの変更は、setwd("ディレクトリ名")
• R プログラムを読みこみ すぐ実行するという命令は source("Rプログラム名")

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

なにかの条件にあてはまる星だけを表示すること。 手順の紹介を用意できなかったが、Rの機能に関心のある人は試みてほしい。

• 条件の例
  • われわれに近い星だけを表示したいならば、年周視差 prlx の値にしきい値をつけてそれ以上のものをえらべばよい。
  • 半径の大きい星だけ、半径の小さい星だけを表示したいならば、 hrhipw5.R にあるように半径 (太陽を1とした相対値) を計算しておき、それの値でえらべばよい。
• 計算手順
  • rhipw5.R, hrhipw6.R のように、星をひとつひとつ points で表示するプログラムならば、 forループの中に points の呼び出しがあるのを、 forループの中に ifによる条件つき実行、その中に points の呼び出しがある形に書きかえる。
  • hrhipw2.R, hrhipw3.R, hrhipw4.R のように、データの配列を plot にわたすプログラムならば、 あらかじめデータを選別し、条件にあうデータの配列をつくっておけばよい。 R のデータフレームでは、条件にあてはまる行だけをとりだして、 もうひとつのデータフレームとしてあつかうことが簡単にできるはずだ (わたしはおぼえておらず、R の教科書的な本の記述をさがす必要があるが)。

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

• 記号の色の変えかたは、ひとまず Okamoto さんの B-V と色の関係を採用したが、 わたしにはその根拠が理解できていない。 Ueno さんのページには、有効放射温度から黒体放射のスペクトルにもとづいて色に変換する式が説明されている。 色についてはそれを採用したほうがよさそうなので、今後プログラミングしてみたい (しかし、いつまでにと約束できない)。
• Rのplotでの対数目盛りで、軸の目盛りのつけかたがうまくいっていない。 対策がわかったら補足したい。
• Rのplotでの、画面の寸法や、画面に対するグラフの相対寸法を変えたい。 Rの機能でできることはわかっているが、具体的に採用することがまだできていない。
• このデータの絶対等級は、実視等級から計算したものなので、実視(V)の波長帯で評価している。 星の物理を考えるうえでは、全波長の放射で評価した「放射絶対等級」にしたほうがよい。 それには、スペクトル型ごとにちがった補正をする。 補正値の概略は『理科年表 2021』 天48ページに「輻射補正」として書かれている。

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

Excelなどの表計算ソフトウェアでも、「散布図」の機能をつかえば、H-R図をつくることができる。 過去にこころみた手順が [(2020年度春学期) : HR図をつくる (その2)]の ページにある。 ただし、そのときとは、入力データの構成を変えてしまったので、 今回の入力データでは作業手順をいくらか変える必要がある。 とくに、実視等級、有効放射温度などの計算をくみこむ必要がある。 Ueno さんは計算を Excel の中で VBA 言語によってやっているが、 わたしは、上記の R プログラムでしているのと同じ計算ならば、 計算式をセルに書きこむ形のほうがよいと思う。 (Ueno さんの色のきめかたを採用したばあいは、やや複雑になるので、プログラムを組む必要があるかもしれない。)

たとえば、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という関数があることを前提としている。)