物理のかぎしっぽ CO/gnuplotでプラネタリウム

CO

gnuplotでプラネタリウム †

星空を見るためのプラネタリウムソフトはたくさんありますが、自分で星空を描けたら素敵ですね。 このページに書いてあることをマスターすれば下のような絵を作ることも可能です。(もっと凝った絵も作れます) :)

CO/gnuplotでプラネタリウム/QandA

display01.png

↑

カタログの入手 †

以下から全天の星についての情報が書かれているカタログ(データファイル)を入手します。 カタログ本体は catalog.dat.gz です。

ftp://dbc.nao.ac.jp/DBC/NASAADC/catalogs/5/5050/

これは Bright Star Catalogue 第5版です。 全天の肉眼で見える程度(mag. < 7.1等)の星(約9000個ほど)に関する情報が記載されています。

ダウンロードしたファイルは圧縮されていますので、解凍コマンドを使って解凍します。

gunzip catalog.dat.gz

解凍するとデータファイル catalog.dat ができています。

データファイルの中身は以下のようになっています。

  1          BD+44 4550      3 36042          46           000001.1+444022000509.9+451345114.44-16.88 6.70  +0.07 +0.08         A1Vn               -0.012-0.018      -018      195  4.2  21.6AC   3
  2          BD-01 4525      6128569                       235956.2-010330000503.8-003011 98.33-61.14 6.29  +1.10 +1.02        gG9                 +0.045-0.060      +014V
  3 33    PscBD-06 6357     281285721002I         Var?     000013.0-061601000520.1-054227 93.75-65.93 4.61  +1.04 +0.89 +0.54   K0IIIbCN-0.5       -0.009+0.089 +.014-006SB1O < 17  2.5   0.0     3*
  4 86    PegBD+12 5063     87 917012004                   000033.8+125023000542.0+132346106.19-47.98 5.51  +0.90               G5III              +0.045-0.012      -002V?
  5          BD+57 2865    123 21085          61  V640 Cas 000101.8+575245000616.0+582612117.03-03.92 5.96  +0.67 +0.20         G5V                +0.263+0.030 +.047-012V          0.8   1.4      *
  6          CD-4914337    142214963      W                000108.4-493751000619.0-490430321.61-66.38 5.70  +0.52 +0.05         G1IV               +0.565-0.038 +.050+003SB         5.7   5.4      *
(以下略)
↑

カタログの整形 †

カタログには星空を描くために必要な情報以外のものも記載されていますので、必要な情報だけを取り出しましょう。 星空を描くために必要な情報は以下です。

  • 位置 (赤経、赤緯)
  • 明るさ
  • 色

カタログは一行につき一天体のデータになっています。 星空のプロットには以下のデータを使用します。

Bytes説明
76-77赤経 時
78-79赤径 分
80-83赤径 秒
84赤緯 角度の符号
85-86赤緯 度
87-88赤緯 分
89-90赤緯 秒
103-107明るさ (等級)
110-114色 (等級)

上記のデータを抜き出す方法はいくつかありますが、ここでは awk というツールを使います。 以下を crop.awk というファイル名で保存してください。

{
    RA_hour = substr($0,76,2);
    RA_min  = substr($0,78,2);
    RA_sec  = substr($0,80,4);
    DEC_deg = substr($0,84,3);
    DEC_arcmin = substr($0,87,2);
    DEC_arcsec = substr($0,89,2);
    Vmag    = substr($0,103,4);
    B_V     = substr($0,110,4);
   
    print RA_hour, RA_min, RA_sec, DEC_deg, DEC_arcmin, DEC_arcsec, Vmag, B_V
}

これをコマンドライン上から次のようにして実行します。

awk -f crop.awk catalog.dat > crop.dat

実行すると、crop.dat というファイルができているはずです。このファイルの中身は次のようになっています。

00 05 09.9 +45 13 45  6.7 +0.0
00 05 03.8 -00 30 11  6.2 +1.1
00 05 20.1 -05 42 27  4.6 +1.0
00 05 42.0 +13 23 46  5.5 +0.9
(以下略)

第一から第三列が赤径、第三から第六列が赤緯、第七列が明るさ(等級)、第八列が色(B-V)です。

なお上の作業をするのが面倒な人は、整形済みのデータを添付しますのでこちらを利用してください。

filecrop.dat

↑

全天のプロット †

まずは全天プロットしてみることにします。

set term postscript eps enhanced colour

# define convert functions
RA2POS( hour, min, sec) = hour + min/60.0 + sec/60.0/60.0
DEC2POS( deg, arcmin, arcsec) = deg + arcmin/60.0 + arcsec/60.0/60.0

set size square

set xlabel "RA"
set ylabel "DEC"

set xrange [24:0]
set yrange [-90:90]

set output "plot01.eps"

plot "crop.dat"\
     using (RA2POS($1,$2,$3)):(DEC2POS($4,$5,$6))\
     notitle\
     with points pt 13 ps 0.3

以上を plot01.gp というファイルに保存してください。コマンドライン上から

gnuplot plot01.gp

とすると plot01.eps というファイルができているはずです。

plot01.png

これは天球面を無理矢理に正方形に直してプロットしたものです。 地球は球面だけど世界地図を長方形に書くのと同じです。

天の川がぼんやりと見えているのがわかるでしょうか :)

↑

オリオン座 †

一部をクローズアップして星座を探してみましょう。 例として有名な星座であるオリオン座をプロットしてみましょう。

オリオン座の位置は理科年表(H.17)によると赤径 5h20m、赤緯 3°です。 そのあたりをクローズアップしてみます。

set term postscript eps enhanced colour

# define convert functions
RA2POS( hour, min, sec) = hour + min/60.0 + sec/60.0/60.0
DEC2POS( deg, arcmin, arcsec) = deg + arcmin/60.0 + arcsec/60.0/60.0

set size square

set xlabel "RA"
set ylabel "DEC"

RA = RA2POS(5,20,0)
DEC = DEC2POS(3,0,0)
delta = 27
dRA = delta / 360.0 * 24
dDEC = delta

set xrange [RA+dRA:RA-dRA]
set yrange [DEC-dDEC:DEC+dDEC]

set xtics 1
set ytics 3

set output "plot02.eps"

plot "crop.dat"\
     using (RA2POS($1,$2,$3)):(DEC2POS($4,$5,$6))\
     notitle\
     with points pt 13 ps 0.7

以上を plot02.gp というファイルに保存してください。コマンドライン上から

gnuplot plot02.gp

とすると、plot02.eps ができているはずです。

plot02.png

このプロットでは異なる明るさの星を同じ点の大きさでプロットしているので、この中からオリオン座を探すのは困難です。 :(
明るさ毎にプロットしてやる必要がありそうです。

↑

明るさ毎にプロット †

明るさ毎にプロットするには、データファイルを明るさ毎にわけるか、gnuplot に何らかの方法でデータをフィルタリングさせる必要があります。ここでは後者の方法を試みます。

set term postscript eps enhanced colour

# define convert functions
RA2POS( hour, min, sec) = hour + min/60.0 + sec/60.0/60.0
DEC2POS( deg, arcmin, arcsec) = deg + arcmin/60.0 + arcsec/60.0/60.0

set size square

set xlabel "RA"
set ylabel "DEC"

RA = RA2POS(5,20,0)
DEC = DEC2POS(3,0,0)
delta = 27
dRA = 24 * delta / 360.0
dDEC = delta

set xrange [RA+dRA:RA-dRA]
set yrange [DEC-dDEC:DEC+dDEC]

set xtics 1
set ytics 3

set output "plot03.eps"

plot "crop.dat"\
       using (($7 < 1.0) ? RA2POS($1,$2,$3) : 1/0):(DEC2POS($4,$5,$6))\
       notitle\
       with points lt 1 pt 13 ps 1.2,\
     ""\
       using ((1.0 <= $7 && $7 < 2.0) ? RA2POS($1,$2,$3) : 1/0):(DEC2POS($4,$5,$6))\
       notitle\
       with points lt 1 pt 13 ps 1.0,\
     ""\
       using ((2.0 <= $7 && $7 < 3.0) ? RA2POS($1,$2,$3) : 1/0):(DEC2POS($4,$5,$6))\
       notitle\
       with points lt 1 pt 13 ps 0.8,\
     ""\
       using ((3.0 <= $7 && $7 < 4.0) ? RA2POS($1,$2,$3) : 1/0):(DEC2POS($4,$5,$6))\
       notitle\
       with points lt 1 pt 13 ps 0.6,\
     ""\
       using ((5.0 <= $7) ? RA2POS($1,$2,$3) : 1/0):(DEC2POS($4,$5,$6))\
       notitle\
       with points lt 1 pt 13 ps 0.2

以上をファイル plot03.gp に保存してください。コマンドライン上から

gnuplot plot03.gp

としてやると、plot03.eps ができているはずです。

plot03.png

どうですか？オリオン座は見つけられましたか？ :)
ちなみに牡牛座やプレアデス星団(すばる)なども見えています。

他の星座も見てみたい方は以下のデータ(理科年表H.17)を参考にプロットしてみてください。

星座赤径赤緯形状
オリオン座5h20m+3°鼓形
カシオペア座1h00m+60°W字形
こぐま座15h40m+78°しっぽに北斗七星
しし座10h30m+15°逆？形
はくちょう座20h30m+43°十字架形
南十字12h20m-60°十字架型
↑

全天プロット (再) †

星の明るさに応じて点の大きさを変える技法を用いて、全天プロットをやり直してみましょう。 描画範囲を全天に変更するだけです。

set term postscript eps enhanced colour

# define convert functions
RA2POS( hour, min, sec) = hour + min/60.0 + sec/60.0/60.0
DEC2POS( deg, arcmin, arcsec) = deg + arcmin/60.0 + arcsec/60.0/60.0

set size square

set xlabel "RA"
set ylabel "DEC"

set xrange [0:24]
set yrange [-90:90]

set xtics 4
set ytics 30

set output "plot06.eps"

plot "crop.dat"\
       using (($7 < 1.0) ? RA2POS($1,$2,$3) : 1/0):(DEC2POS($4,$5,$6))\
       notitle\
       with points lt 1 pt 13 ps 1.2,\
     ""\
       using ((1.0 <= $7 && $7 < 2.0) ? RA2POS($1,$2,$3) : 1/0):(DEC2POS($4,$5,$6))\
       notitle\
       with points lt 1 pt 13 ps 1.0,\
     ""\
       using ((2.0 <= $7 && $7 < 3.0) ? RA2POS($1,$2,$3) : 1/0):(DEC2POS($4,$5,$6))\
       notitle\
       with points lt 1 pt 13 ps 0.8,\
     ""\
       using ((3.0 <= $7 && $7 < 4.0) ? RA2POS($1,$2,$3) : 1/0):(DEC2POS($4,$5,$6))\
       notitle\
       with points lt 1 pt 13 ps 0.6,\
     ""\
       using ((4.0 <= $7) ? RA2POS($1,$2,$3) : 1/0):(DEC2POS($4,$5,$6))\
       notitle\
       with points lt 1 pt 13 ps 0.1

以上をファイルに保存して、プロットを実行して下さい。 以下のようになるでしょう。

plot06.png

見やすくなりましたね☆ ちなみにこの図から星座を見つけ出すことのできる人もいるでしょう。 北斗七星なんかは分かりやすい・・かな？

↑

天球儀 (全天の3Dプロット) †

星を球面座標上に描いてみることにしましょう。

set term postscript eps enhanced colour

# define convert functions
RA2RAD( hour, min, sec) = 2*pi/24.0*(hour + min/60.0 + sec/60.0/60.0)
DEC2RAD( deg, arcmin, arcsec) = 2*pi/360.0*(deg + arcmin/60.0 + arcsec/60.0/60.0)

set size square

unset xtics
unset ytics
unset ztics

set mapping spherical
set parametric
set isosamples 20,20
set ticslevel 0

set hidden3d 

# To see the Orion
set view 75, 150, 1, 1

set urange [-pi:pi]
set vrange [-pi:pi]

set output "plot04.eps"
splot sin(u)*cos(v),sin(u)*sin(v),cos(u)\
        notitle\
        with lines lt 3 lw 2,\
      "crop.dat"\
        using (($7 < 1.0) ? RA2RAD($1,$2,$3) : 1/0):(DEC2RAD($4,$5,$6)):(1)\
        notitle\
        with points lt 1 pt 13 ps 1.0,\
      ""\
        using ((1.0 <= $7 && $7 < 2.0) ? RA2RAD($1,$2,$3) : 1/0):(DEC2RAD($4,$5,$6)):(1)\
        notitle\
        with points lt 1 pt 13 ps 0.8,\
      ""\
        using ((2.0 <= $7 && $7 < 3.0) ? RA2RAD($1,$2,$3) : 1/0):(DEC2RAD($4,$5,$6)):(1)\
        notitle\
        with points lt 1 pt 13 ps 0.6,\
      ""\
        using ((3.0 <= $7 && $7 < 4.0) ? RA2RAD($1,$2,$3) : 1/0):(DEC2RAD($4,$5,$6)):(1)\
        notitle\
        with points lt 1 pt 13 ps 0.4,\
      ""\
        using (4.0<=$7 ? RA2RAD($1,$2,$3) : 1/0):(DEC2RAD($4,$5,$6)):(1)\
        notitle\
        with points lt 1 pt 13 ps 0.1

上記をファイル plot04.gp に保存して実行してください。 plot04.eps が生成されます。

plot04.png

なお、gnuplot 4.0 以降では画面に表示させた 3D グラフをマウスでドラッグすることにより視点を変更することができます。 次のスクリプトを実行すれば、クリクリ回転させて遊ぶことができるでしょう。 (改良中)

#!/usr/bin/gnuplot

# define convert functions
RA2RAD( hour, min, sec) = 2*pi/24.0*(hour + min/60.0 + sec/60.0/60.0)
DEC2RAD( deg, arcmin, arcsec) = 2*pi/360.0*(deg + arcmin/60.0 + arcsec/60.0/60.0)

set mapping spherical
set parametric
set isosamples 20,20
set ticslevel 0

set hidden3d

set urange [-pi:pi]
set vrange [-pi:pi]

splot sin(u)*cos(v),sin(u)*sin(v),cos(u)\
        notitle\
        with lines lt 3 lw 2,\
      "crop.dat"\
        using (($7 < 1.0) ? RA2RAD($1,$2,$3) : 1/0):(DEC2RAD($4,$5,$6)):(1)\
        notitle\
        with points lt 1 pt 13 ps 2.0,\
      ""\
        using ((1.0 <= $7 && $7 < 2.0) ? RA2RAD($1,$2,$3) : 1/0):(DEC2RAD($4,$5,$6)):(1)\
        notitle\
        with points lt 1 pt 13 ps 1.7,\
      ""\
        using ((2.0 <= $7 && $7 < 3.0) ? RA2RAD($1,$2,$3) : 1/0):(DEC2RAD($4,$5,$6)):(1)\
        notitle\
        with points lt 1 pt 13 ps 1.3,\
      ""\
        using ((3.0 <= $7 && $7 < 4.0) ? RA2RAD($1,$2,$3) : 1/0):(DEC2RAD($4,$5,$6)):(1)\
        notitle\
        with points lt 1 pt 13 ps 1.0
#      ""\
#        using (4.0<=$7 ? RA2RAD($1,$2,$3) : 1/0):(DEC2RAD($4,$5,$6)):(1)\
#        notitle\
#        with points lt 1 pt 13 ps 0.2



pause -1 "Hit the return key!"
添付ファイル: fileplot01.png 413件 [詳細] filecrop.dat 354件 [詳細] fileplot02.png 342件 [詳細] fileplot03.png 326件 [詳細] fileplot04.png 332件 [詳細] filedisplay01.png 314件 [詳細] fileplot06.png 327件 [詳細]
トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Modified by 物理のかぎプロジェクト PukiWiki 1.4.6 Copyright © 2001-2005 PukiWiki Developers Team. License is GPL.
Based on "PukiWiki" 1.3 by yu-ji Powered by PHP 5.3.29 HTML convert time to 0.025 sec.