地学概論

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

恒星などの天体について、 天体から届く放射照度とみかけの等級、天体の放射光度と絶対等級、 天体の有効放射温度(光球の表面温度)とスペクトル型が、それぞれ対応していることを知っておきたい。

HR (ヘルツスプルング Hertzsprung - ラッセル Russell) 図は、 絶対等級と、スペクトル型・色指数・有効放射温度とを2軸にとった図で、 星を「主系列星」「赤色巨星」などに分類して見るのに使われる。

ここでは、HR図を作ることをひとつの目標として、 それにあらわれる量に慣れるような実習を、2回かけてやることにする。 第1回は、少量のデータを読みこんで、単純な操作で、かならずしもきれいではない図にしてみるところまで。

参考教材例

HR図を作成する高校レベルの教材は、あちこちでつくられている。

西浦 慎悟さん(大学の教育学部の天文の教員) のものは、 複数の波長帯でそれぞれとられた天体の画像 (写真のようなもの) から、あかるさを読みとるところから始めている。 天体観測そのものではないが、その次の段階から、研究者がするのに近い体験をすることができる。 しかし、残念ながら、わたしにはそれを指導するのがむずかしい。

Ueno Satoshi さん (小学校教員らしい) は、 Hipparcos星表 【ヒッパルコスは古代ギリシャの天文学者で星座をきめた人だとされるが、 ここでは現代のヨーロッパ宇宙機構 (ESA) による恒星位置測定用衛星の名まえになっており、 Hipparcos星表はその衛星を使った研究プロジェクトの成果である】の データを出発点として、ExcelのVBAで加工し作図している。 有効放射温度にもとづいて実際の星の色に近い色を図示することもしている。 しかし、Uenoさんは、専門的内容には自信がないとものべている。 Hipparcos星表を使う部分は、次のOkamotoさんのやりかたをまねたものである。

Okamoto Yoshio (岡本 義雄) さん (高校の地学の教員) は、 Hipparcos星表のデータを awk で加工し、 GMT (Generic Mapping Tools)というUnix系OSで動く作図ソフトウェアで作図している。 (ただしそのサイトは 2020年6月には見つからなくなってしまった。 おそらく退職とともにサイトがなくなったのだろう。) 「色指数」にもとづいて実際の星の色に近い色を図示することもしている。 ただし色のきめかたの根拠は書いてない。(リンクさきのページも見つからなくなった。) わたしはこれにならって、GMT のかわりに R で作図する教材を 2019年度に作成し、2020年度に更新した。暫定的にそのページ [2020年度実習: HR図をつくる (2)] をリンクしておく。 次回に、それを改訂した教材による実習を予定している。


みかけの等級、年周視差、有効放射温度からの、HR図の初歩的な作図

材料データ

理科年表 2021』 の「天49」ページ「恒星の物理的諸量(2)」の 「物理的諸量の実例」の表の内容を、編集してここに置いた。

【これは、『理科年表 2019』の表を書きぬいたものをもとに、 『理科年表 2020』を参照して変更点をなおし、さらに 『理科年表 2021』を参照して変更点をなおしたものである。 『理科年表 2022』が出たら (2021年11月末ごろの予定) それと照合して、変更点があれば、ファイル名を変えて追加することにしたい。】

最初の行に入れた各列の見出しは R で変数名として使われるので、 (Rは日本語文字にも対応しているが、文字コード化けの心配があるので) 内容を、アルファベット・数字・下線文字だけで構成した。 見出しの項目の、教材用データと理科年表の表との対応は、次のとおり。

(補足説明)


作業準備

自分が作業するパソコンに、この実習のための 作業用ディレクトリ (フォルダ) を用意しよう。 【ディレクトリ名は、わたしの授業に関するかぎり、日本語文字がはいっていてもよいが、 R Studio をつかうならば、アルファベット・数字などのASCII文字にかぎったほうがよいと聞いている。】

rika2021ten49.txtを、作業用ディレクトリ にダウンロードしよう。


Rによる作業手順

上記のデータをRプログラムの次のような文で、「ten49df」という「データフレーム」に読みこむことができる。

ten49df <- read.table("rika2021ten49.txt", header=TRUE)

読みこまれると、表の各列は、「ten49df$mass」のようなベクトル変数として参照することができる。 Rのコンソール画面に変数名を書けば、その変数の値が表示される。

絶対等級 [absolute magnitude] M は、見かけの等級 m と、 距離 d [ただし単位はパーセクとする] または 年周視差 p [単位は秒角] と から、次のように計算することができる。 これは、等級の近代的定義、絶対等級の定義と、放射照度が距離の2乗に反比例することから導かれる。

M = m + 5 - 5 log10 d

M = m + 5 + 5 log10 p

(『図表』 1章 「恒星の明るさと色」38ページ [3] 「絶対等級」に式がある。 ただしそこでは、距離 d の単位を光年としている。 『理科年表 2021』では、天47ページ「天体の距離と絶対等級」に式がある。 ただしそこでは、年周視差を π [パイ] の字であらわしている。円周率でないことに注意。)

Rでは、ベクトル変数の演算式で、一連のデータをいっせいに変換することができる。 (結果の amag は ten49df の部分ではなく別のベクトル変数になる。)

amag <- ten49df$vmag + 5 + 5*log10(ten49df$prlx / 1e3)

【理科年表からとった表では、年周視差の数値が、10-3 秒角 を単位としてあらわされていて、 上の M の式の p の単位は秒角なので、log10の中に「/1e3」(「10の3乗で割る」)を入れている。】

横軸に有効放射温度、縦軸に絶対等級をとってプロットしてみる。

plot(ten49df$temp_eff, amag)

縦軸、横軸それぞれの向きを標準的HR図に合わせるように、図のx、yの範囲指定をする。

plot(ten49df$temp_eff, amag, xlim=c(30000, 0), ylim=c(15,-5))

plot の中で、図の表題は main="..." 、x軸のラベルは xlab="..." 、 y軸のラベルは ylab="..." で指定することができる。 次のかぎかっこ内の文字列を入れてみよう。

([注] Windows 7での R で日本語文字を入れるには、文字コードは「Shift_JIS」にする。 他のOSでは「UTF-8」にする必要があるかもしれない。)

目盛りや格子線の入れかたをいろいろくふうすることができる。

星の大きさ (半径) を表現する

plot の中で、cex に値をあたえることによって、記号(〇)の大きさを変えることができる。 星の半径によって変えてみることにする。 Okamoto さんの「天文教材作成備忘録(パート2)」[2019-10-25 閲覧] を参考にして (しかし、そのままではなく)、 cex にあたえる値を、星の半径 (太陽を1とした相対値) の 0.04倍、 ただし、それが 1未満のときは 1、としてみた。

symsize0 <- ten49df$radius*0.04

【授業時には「ten40df」となっていましたが正しくは「ten49df」です。授業後訂正しました。】

symsize <- ifelse(symsize0 > 1, symsize0, 1)

plot(ten49df$temp_eff, amag, xlim=c(30000, 0), ylim=c(15,-5), cex=symsize)

星の名まえを表示する

上記の plot のあとで、つぎのようにすれば、 それぞれの星の名まえを、星をあらわす記号の付近に文字列で表示することができる。 文字の色を青 (「col=4」で指定) にしてみた。

text(ten49df$temp_eff, amag, ten49df$name, cex=0.7, col=4)

文字列どうしがかさなるのをさけて、文字列が読めるようにするために、 ひとまず、文字の大きさを cex=0.7 として指定したが、まだかさなって読めないところがある。 読みやすくするには、さらにくふうが必要である。


表計算ソフトウェアによる作業手順

読みこみ

今回用意した rika2021ten49.txt は、「スペースくぎりのテキストファイル」である。 ただし、複数のスペース(空白)がつづいても、ひとつのくぎりとみなす。 表計算ソフトウェアは、このようなファイルのデータを表として読みこむことができるが、 その手順は、ソフトウェアのバージョンごとにちがっている可能性がある。 ここでは、Windows 10 上で、 Microsoft Excel 2013と、LibreOffice Calc version 6 でやってみた手順をしめす。 rika2021ten49.txt がダウンロードずみであることを前提とする。

Microsoft Excel 2013 での読みこみ

まず、ブックを作成し、空白のワークシートを開いた状態にする。

上のメニューの「ファイル」から「開く」を選択する。 「ファイルを開く」のメニューが出る。 「ファイル名」の右のファイルの種類を指定するところを「すべてのファイル」にしておく (これが「Excelファイル」などだと、".txt" でおわるファイル名のものは見えない)。 フォルダーのツリーをたどって、rika2021ten49.txt を選択する。

「テキストファイルウィザード」が開く。

「次へ」を押して、くぎられかたを確認し、「完了」を押す。 ワークシートの A1 から J32 までがうまれば、読みこみ成功。

LibreOffice Calc version 6 での読みこみ

上のメニューの「ファイル」から「開く」を選択。 フォルダーのツリーをたどって、rika2021ten49.txt を選択する。

「テキストのインポート」のメニューが出る。

これで「OK」を押す。ワークシートの A1 から J32 までがうまれば、読みこみ成功。

絶対等級を計算する (Excel, LibreOfficeに共通)

セル K2 に、式「=D2+5+5*LOG10(E2/1e3)」を入れる。 (D列がRの例文の「vmag」、E列が「prlx」にあたる。)

【理科年表からとった表では、年周視差の数値が、10-3 秒角 を単位としてあらわされていて、 上の M の式の p の単位は秒角なので、log10の中に「/1e3」(「10の3乗で割る」)を入れている。】

セル K2 の内容を「コピー」し、K3:K32 に「はりつける」。 ただし、数値でなく計算式が複写されるようにする。 【ここでの計算式の中でのセル番号の参照のしかたは「相対参照」なので、 複写とともに自動的に、ちがったセルをさすように書きかえられる、 ここでは意識的にその機能をつかう。】

この結果、セルK2には計算エラーが出るはず。 これは、太陽の年周視差の値が定義できないので、 セルE2 に数値でなく「NA」がはいっているからである。 今回の課題では、太陽をふくめないで作図してみることにする。 もし太陽をふくめたければ、太陽の絶対等級の値をしらべて、 セル K2 に(計算式でなく) その数値を書きこんでしまえばよい。

作図 (暫定)

(考えかたはつぎのとおり。 実際にマウスやキーボードをどう操作すればよいかは、わたしにはまだ説明できない。)

「グラフ」の「散布図」の機能を利用する。

横軸に、I3:I32 (太陽の絶対等級を追加したばあいは I2:I32) をとる。 左にいくほど大きい値になるようにする。

縦軸に、K3:K32 (太陽の絶対等級を追加したばあいは K2:K32) をとる。 上にいくほど小さい値になるようにする。


発展: 縦軸の変数を絶対等級から光度にかえる

縦軸を「光度」(luminosity)の太陽を1とした相対値 (変数名の案: rel_luminos )に変換して表示することを考えてみよう。 絶対等級は光度の対数に比例し、絶対等級の差 5 が、光度の100倍に相当する。 太陽の絶対等級は 4.83 とされている。

Rでは、次のようにすればよさそうだ。(ただし、まだよく検算していない。)

rel_luminos <- 10^((4.83-amag)*2/5)

ここで、記号「^」を使った「a^b」のような演算式は「aのb乗]を意味する。 また、2が出てくるのは 100が10の2乗だからである。 グラフを縦軸だけ対数尺度にするには、plot のパラメータとして log="y" とする。 なお、ylim も適切に変更する。

表計算でも同様に、絶対等級が K 列にはいっているならば、 セル L2 につぎのような計算式を書いて、下にコピーペーストすればよさそうだ。ただし検算していない。

=10^((4.83-K2)*2/5)


発展: データの追加

次回にもっとたくさんの恒星のデータをあつかう予定ですが、 そのまえに事例をふやしてみたい人のために:

[(*) 図表 41ページ [4] 「恒星の諸量」] に収録された恒星は、 rika2021ten49.txt とかさなるものが多いが、そこにふくまれていないものもある。 ふくまれていなかったものを追加して、データ件数をふやしてみよう。 ただし、表計算ワークシートをつかうばあい、セルに数値がはいっているか、計算式がはいっているか、意識して、 とりちがえないように注意してください。


残された問題: 等級をきめる波長帯の選択

ここで使った見かけの等級は実視等級なので、 今回表示した絶対等級も「実視」の波長帯で決めたもの (実視絶対等級) になっている。 天体自体を論じるときには、放射絶対等級で論じたほうがよいだろう。 実視絶対等級から放射絶対等級を推定するには、スペクトル型ごとにちがった補正値をくわえればよい。 『理科年表 2021』の「天48」ページ「恒星の物理的諸量(1)」の表に、「輻射補正」として示されている。 これを考慮に入れることは、もうすこし進んだレベルの課題になり、わたしの理解はまだ不じゅうぶんである。


R言語を使ってみる

ギリシャ文字一覧

文献

参考サイト


(2019-10-28, 2020-06-19 の教材をもとに改訂) 2021-11-07
増田 耕一 (MASUDA Kooiti)
[「地学概論」日程と各回の内容] [今学期教材入口ページ] [地学の教材目次ページ]