* gnuplotでプラネタリウム [#f8934c49]
星空を見るためのプラネタリウムソフトはたくさんありますが、自分で星空を描けたら素敵ですね。
このページに書いてあることをマスターすれば下のような絵を作ることも可能です。(もっと凝った絵も作れます) :)
&ref(display01.png);
この記事は CO が書きました。
#contents
** カタログの入手 [#q70e6026]
以下から全天の星についての情報が書かれているカタログ(データファイル)を入手します。
カタログ本体は 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 *
(以下略)
** カタログの整形 [#a139d1ed]
カタログには星空を描くために必要な情報以外のものも記載されていますので、必要な情報だけを取り出しましょう。
星空を描くために必要な情報は以下です。
- 位置 (赤経、赤緯)
- 明るさ
- 色
カタログは一行につき一天体のデータになっています。
星空のプロットには以下のデータを使用します。
,Bytes,説明
,76-77,赤経 時
,78-79,赤径 分
,80-83,赤径 秒
,84,赤緯 角度の符号
,85-86,赤緯 度
,87-88,赤緯 分
,89-90,赤緯 秒
,103-107,明るさ (等級)
,110-114,色 (等級)
// 76- 77 I2 h RAh ?Hours RA, equinox J2000, epoch 2000.0 (1)
// 78- 79 I2 min RAm ?Minutes RA, equinox J2000, epoch 2000.0 (1)
// 80- 83 F4.1 s RAs ?Seconds RA, equinox J2000, epoch 2000.0 (1)
// 84 A1 --- DE- ?Sign Dec, equinox J2000, epoch 2000.0 (1)
// 85- 86 I2 deg DEd ?Degrees Dec, equinox J2000, epoch 2000.0 (1)
// 87- 88 I2 arcmin DEm ?Minutes Dec, equinox J2000, epoch 2000.0 (1)
// 89- 90 I2 arcsec DEs ?Seconds Dec, equinox J2000, epoch 2000.0 (1)
// 103-107 F5.2 mag Vmag ?Visual magnitude (1)
// 110-114 F5.2 mag B-V ? B-V color in the UBV system
上記のデータを抜き出す方法はいくつかありますが、ここでは 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)です。
なお上の作業をするのが面倒な人は、整形済みのデータを添付しますのでこちらを利用してください。
&ref(crop.dat);
** 全天のプロット [#gcc50a74]
まずは全天プロットしてみることにします。
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 + (deg < 0 ? -arcmin : arcmin)/60.0 + (deg < 0 ? -arcsec : 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 というファイルに保存してください。コマンドライン上から
&ref(plot01.gp);
gnuplot plot01.gp
とすると plot01.eps というファイルができているはずです。
&ref(plot01.png);
これは天球面を無理矢理に正方形に直してプロットしたものです。
地球は球面だけど世界地図を長方形に書くのと同じです。
天の川がぼんやりと見えているのがわかるでしょうか :)
** オリオン座 [#c36eaa58]
一部をクローズアップして星座を探してみましょう。
例として有名な星座であるオリオン座をプロットしてみましょう。
オリオン座の位置は理科年表(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 + (deg < 0 ? -arcmin : arcmin)/60.0 + (deg < 0 ? -arcsec : 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 というファイルに保存してください。コマンドライン上から
&ref(plot02.gp);
gnuplot plot02.gp
とすると、plot02.eps ができているはずです。
&ref(plot02.png);
このプロットでは異なる明るさの星を同じ点の大きさでプロットしているので、この中からオリオン座を探すのは困難です。 :(~
明るさ毎にプロットしてやる必要がありそうです。
*** 明るさ毎にプロット [#fc55f1e4]
明るさ毎にプロットするには、データファイルを明るさ毎にわけるか、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 + (deg < 0 ? -arcmin : arcmin)/60.0 + (deg < 0 ? -arcsec : 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 に保存してください。
&ref(plot03.gp);
コマンドライン上から
gnuplot plot03.gp
としてやると、plot03.eps ができているはずです。
&ref(plot03.png);
どうですか？オリオン座は見つけられましたか？ :)~
ちなみに牡牛座やプレアデス星団(すばる)なども見えています。
他の星座も見てみたい方は以下のデータ(理科年表H.17)を参考にプロットしてみてください。
,星座,赤径,赤緯,形状
,オリオン座,5h20m,+3°,鼓形
,カシオペア座,1h00m,+60°,W字形
,こぐま座,15h40m,+78°,しっぽに北斗七星
,しし座,10h30m,+15°,逆？形
,はくちょう座,20h30m,+43°,十字架形
,南十字,12h20m,-60°,十字架型
** 全天プロット (再) [#m9edc39d]
星の明るさに応じて点の大きさを変える技法を用いて、全天プロットをやり直してみましょう。
描画範囲を全天に変更するだけです。
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 + (deg < 0 ? -arcmin : arcmin)/60.0 + (deg < 0 ? -arcsec : 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 "plot04.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
以上をファイルに保存して、プロットを実行して下さい。
&ref(plot04.gp);
以下のようになるでしょう。
&ref(plot04.png);
見やすくなりましたね☆
ちなみにこの図から星座を見つけ出すことのできる人もいるでしょう。
北斗七星なんかは分かりやすい・・かな？
** 天球儀 (全天の3Dプロット) [#vd6acd7d]
星を球面座標上に描いてみることにしましょう。
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 + (deg < 0 ? -arcmin : arcmin)/60.0 + (deg < 0 ? -arcsec : 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 "plot05.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
上記をファイル plot05.gp に保存して実行してください。
&ref(plot05.gp);
plot05.eps が生成されます。
&ref(plot05.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 + (deg < 0 ? -arcmin : arcmin)/60.0 + (deg < 0 ? -arcsec : 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!"
