ヒッパルコス星表で太陽近傍恒星のHR図を描こう

ヒッパルコス星表からエクセルを使って太陽近傍恒星のHR図(ヘルツシュプルング・ラッセル図)を作図します

Satoshi Ueno提供サイト

あなたの名前の世界旅行…「世界の文字で遊ぼう」
小学校６年理科「土地のつくりと変化」デジタル教材…「発見！地球号～滋賀の大地」
かんたんに星座を見つけたい人は…「星座早見を使おう」
滋賀県の市町名を楽しく学びたい人は…「滋賀県地図パズルで遊ぼう」

ＨＲ図の作図

(1)グラフ挿入（散布図）で作成

B-V色指数と視差角について、標準誤差から精度の低いデータをはじきます。（今回は誤差10％未満）

これにより、恒星データ数約20000となります。

「B-Vの列」と「絶対光度の列」をグラフ範囲に指定して、グラフの挿入を選択します。XY軸の範囲やマーカーの種類等を適宜変更します。単純に恒星をプロットするだけならこれで完成です。

(2)VBAを利用して作成（恒星の半径と色を反映）

恒星の半径と色を反映させたＨＲ図を描くためにVBAを利用します。ワークシート上に「B-V色指数」と「絶対光度」によりオートシェイプの「円/楕円」を挿入します。「円/楕円」の大きさと色に「太陽との半径比R/Rs」と「表面温度Tから求めたRGB値」を指定します。

※今回使用したPC環境では、恒星20,000個のプロットに分単位の時間を要しました。

恒星の半径と色を反映させたＨＲ図のためのVBA の例

(1)エクセルの入力例

恒星をオートシェープで描画する際に、恒星同士の重なりが分かりやすいように「線」と「塗り」で輝度を変えてみました。さらに、表面温度（最大波長）からの色の計算に２種類の関数を作成して比較してみたところ、結果はほぼ同じでした。（aネットで見つけた近似式利用,bプランクの放射則とCIEの等色関数から積分により計算）

恒星データシート名：HIPデータ

HR図の作成マクロ：drawHR2()

恒星色の推定関数：getRGB_vba(表面温度,輝度)※ネット上で見つけた「xy色度図上の黒体軌跡の近似式」を利用

恒星色の推定関数：getRGB_cie(表面温度,輝度)※黒体放射の分光分布とCIEの等色関数の積を可視波長域で積分して計算

エクセルシートの例（HR図描画マクロで使用するデータ部分）
	A	T	AB(28)	AC(29)	AD(30)	AE	AF(32)	AG(33)	AH(34)	AI(35)
1	HIP番号	表面温度T	B-V色指数	絶対等級M	半径比R/Rs	色hex	色(塗り)	色(線)	色(塗り)	色(線)
2	7	5660.377	0.74	5.884768	0.648951	=getRGB_hex(T2)	=getRGB_vba(T2,0.8)	=getRGB_vba(T2,0.4)	=getRGB_cie(T2,0.8)	=getRGB_cie(T2,0.4)

※列番号（A～）に添えた（数字）は、マクロ内の Cells() メソッドの引数として指定の値です。

(2)恒星の半径、色を反映させたHR図を作成するマクロ

グラフ用シート名：HR図2

一等星名シート名：一等星

Sub drawHR2()
  Dim  i  As Long
  Dim  x, y  As Integer					‘オートシェイプの挿入座標
  Dim  r_rs  As Single					‘半径比
  Dim  x0, x1, x2, y0, y1, y2, xm, ym  As Integer	‘オートシェイプの座標計算用
  Dim  rgb_f,rgb_l  As Long				‘オートシェイプの塗り色、線色（ｒｇｂ値）
  Dim  name  As String
  Dim  oShape  As Shape					‘オートシェイプの一時保管用
 'セルをHR図の軸の目盛として利用するためにセル幅、高さを調整する
 '標準フォントを設定（セル幅の設定単位とptの関係を一定にするため）
  Application.StandardFont = "ＭＳ ゴシック"
  Application.StandardFontSize = 12
  xm = 180: ym = 24					‘ｘ軸、ｙ軸の表示倍率　１→(pt)
  x1 = 20: x2 = 20: y1 = 50: y2 = 30			‘1,2行、1,2列のセル幅、セル高の調整用(pt)
  x0 = x1 + x2 + xm * 0.5: y0 = y1 + y2 + ym * 10	‘グラフの原点の位置
 'セル高をｙ軸の目盛単位（5等級）に合わせて設定
  Worksheets("HR図2").Rows(1).RowHeight = y1
  Worksheets("HR図2").Rows(2).RowHeight = y2
  For i = 3 To 8
    Worksheets("HR図2").Rows(i).RowHeight = ym * 5
  Next
 'セル幅をｘ軸の目盛単位（0.5）に合わせて設定
 '（○-3.75）/6の部分は、ptをセル幅独特の単位に変換するため（excelの設定によっては変更する必要）
  Worksheets("HR図2").Columns(1).ColumnWidth = (x1 - 3.75) / 6
  Worksheets("HR図2").Columns(2).ColumnWidth = (x2 - 3.75) / 6
  For i = 3 To 8
    Worksheets("HR図2").Columns(i).ColumnWidth = (xm * 0.5 - 3.75) / 6
  Next
 'オートシェイプの円と一等星名テキストボックスをすべて削除
  For Each oShape In Worksheets("HR図2").Shapes
      If oShape.AutoShapeType = msoShapeOval Or (oShape.Type = msoTextBox And oShape.name <> "表題") Then oShape.Delete
  Next
 '恒星のプロット
  For i = 2 To 20178
    r_rs = Worksheets("HIPデータ").Cells(i, 30).Value + 1
    rgb_f = Worksheets("HIPデータ").Cells(i, 32).Value
    rgb_l = Worksheets("HIPデータ").Cells(i, 33).Value
    x = Worksheets("HIPデータ").Cells(i, 28).Value * xm + x0
    y = Worksheets("HIPデータ").Cells(i, 29).Value * ym + y0
    Worksheets("HR図2").Shapes.AddShape(msoShapeOval, x - r_rs / 2, y - r_rs / 2, r_rs, r_rs).Select
    Selection.ShapeRange.Fill.Solid
    Selection.ShapeRange.Fill.ForeColor.rgb = rgb_f
    Selection.ShapeRange.Line.Weight = 0.5
    Selection.ShapeRange.Line.ForeColor.rgb = rgb_l
  Next
 '一等星名（太陽）の表示
 For i = 2 To 23
    x = Worksheets("一等星").Cells(i, 3).Value * xm + x0
    y = Worksheets("一等星").Cells(i, 4).Value * ym + y0
    name = Worksheets("一等星").Cells(i, 2).Value
    Worksheets("HR図2").Shapes.AddTextbox(msoTextOrientationHorizontal, x - 6, y - 6, 12, 12).Select
    Selection.Characters.Text = "・" + name
    Selection.ShapeRange.Fill.Visible = msoFalse
    Selection.ShapeRange.Line.Visible = msoFalse
    Selection.AutoSize = True
  Next
‘予めシート上に「表題」という名前のテキストボックスを配置しておく
  Worksheets("HR図2").Shapes("表題").ZOrder msoBringToFront　　
End Sub

(3)表面温度ｔからｒｇｂを求める関数の例(a)

※ｘｙ色度図上の黒体軌跡の近似式を利用

'VBAのRGB()関数と同じ形式の値を返す関数
Function getRGB_VBA(r,Y)　　　
   Dim rgb_arry()
   rgb_arry = getRGB(r,Y)
   getRGB_VBA = rgb_arry(0) + rgb_arry(1)*256 + rgb_arry(2)*256*256
End Sub

'rgbの各値を配列として返す関数（getRGB_VBA()関数から呼び出し）
Function getRGB(t,Y)
'単純に表面温度（最大波長）のみで色を推定
  Dim  i  As Integer
  Dim  x0, y0, X, Z  As Double
  Dim  rgb(3),  mtrx(9)
'sRGB(CIE 61966-2-1)
  mtrx(0) = 3.2410       :   mtrx(1) = -1.5374    :  mtrx(2) = -0.4986
  mtrx(3) = -0.9692      :   mtrx(4) = 1.8760     :  mtrx(5) = 0.0416
  mtrx(6) = 0.05563      :   mtrx(7) = -0.2040    :  mtrx(8) = 1.0570
'ｘｙ空間の黒体軌跡の近似式（Wikipedia「Planckian locus」Approximationによる）
  If t < 4000 Then
     x0 = -0.2661239 * (10 ^ 9 / t ^ 3)
     x0 = x0 - 0.234358 * (10 ^ 6 / t ^ 2)
     x0 = x0 + 0.8776956 * (10 ^ 3 / t)
     x0 = x0 + 0.17991
  Else
     x0 = -3.0258469 * (10 ^ 9 / t ^ 3)
     x0 = x0 + 2.1070379 * (10 ^ 6 / t ^ 2)
     x0 = x0 + 0.2226347 * (10 ^ 3 / t)
     x0 = x0 + 0.24039
  End If
  If t < 2222 Then
     y0 = -1.1063814 * (x0 ^ 3)
     y0 = y0 - 1.3481102 * (x0 ^ 2)
     y0 = y0 + 2.18555832 * x0
     y0 = y0 - 0.20219683
  ElseIf t >= 2222 And t < 4000 Then
     y0 = -0.9549476 * (x0 ^ 3)
     y0 = y0 - 1.37418593 * (x0 ^ 2)
     y0 = y0 + 2.09137015 * x0
     y0 = y0 - 0.16748867
  Else
     y0 = 3.081758 * (x0 ^ 3)
     y0 = y0 - 5.8733867 * (x0 ^ 2)
     y0 = y0 + 3.75112997 * x0
     y0 = y0 - 0.37001483
  End If
'Yxy→XYZ
  X = (Y / y0) * x0
  Z = (Y / y0) * (1 - x0 - y0)
'XYZ→sRGB
  rgb(0) = mtrx(0) * X + mtrx(1) * Y + mtrx(2) * Z
  rgb(1) = mtrx(3) * X + mtrx(4) * Y + mtrx(5) * Z
  rgb(2) = mtrx(6) * X + mtrx(7) * Y + mtrx(8) * Z
  For i = 0 To 2
     If rgb(i) < 0 Then　rgb(i) = 0
'ガンマ補正γ＝2.2
     rgb(i) = rgb(i) ^ (1 / 2.2)	
     rgb(i) =　Round(rgb(i) * 255)
     If rgb(i) > 255 Then　rgb(i) = 255
  Next
  getRGB = rgb
End Function

(4)表面温度ｔからｒｇｂを求める関数の例(b)

※プランクの放射則とCIEの等色関数を利用

Function getRGB_CIE(t, YY)
   Dim  X,  Y,  Z,  M ,c2　 As Double
   Dim  x0, y0  As Single
   Dim  i  As Integer
   Dim  CIE_x(94),  CIE_y(94),  CIE_z(94),  rmd(94),  rgb(3),  mtrx(9)
'sRGB(CIE 61966-2-1)
   mtrx(0) = 3.2406:  mtrx(1) = -1.5372:  mtrx(2) = -0.4986
   mtrx(3) = -0.9689:  mtrx(4) = 1.8758:  mtrx(5) = 0.0415
   mtrx(6) = 0.0557:  mtrx(7) = -0.204:  mtrx(8) = 1.057
'∫Me（λ）・x_y_z_（λ）dλ
   c2 = 1.438775 * 10 ^ 7
   For  i = 0 To 94 Step 2
      rmd(i) = Worksheets("等色関数").Cells(i + 2, 1).Value
      CIE_x(i) = Worksheets("等色関数").Cells(i + 2, 2).Value
      CIE_y(i) = Worksheets("等色関数").Cells(i + 2, 3).Value
      CIE_z(i) = Worksheets("等色関数").Cells(i + 2, 4).Value
      M = (rmd(i) ^ 5) * (Exp(c2 / rmd(i) / t) - 1)
      X = X + CIE_x(i) / M
      Y = Y + CIE_y(i) / M
      Z = Z + CIE_z(i) / M
   Next
   x0 = X / (X + Y + Z)
   y0 = Y / (X + Y + Z)
' xyY→XYZ
   X = (YY / y0) * x0
   Z = (YY / y0) * (1 - x0 - y0)
' XYZ→sRGB
   rgb(0) = mtrx(0) * X + mtrx(1) * YY + mtrx(2) * Z
   rgb(1) = mtrx(3) * X + mtrx(4) * YY + mtrx(5) * Z
   rgb(2) = mtrx(6) * X + mtrx(7) * YY + mtrx(8) * Z
   For i = 0 To 2
       If rgb(i) < 0 Then  rgb(i) = 0
'ガンマ補正（γ=2.2）
       rgb(i) = rgb(i) ^ (1 / 2.2)
       rgb(i) = Round(rgb(i) * 255)
       If  rgb(i) > 255 Then rgb(i) = 255
   Next
   getRGB_CIE = rgb(0) + rgb(1) * 256 + rgb(2) * 256 * 256
End Function

←前のステップへ|次のステップへ→