物理のかぎしっぽ CO/gnuplotでプラネタリウム のバックアップソース(No.7)
[[CO]]

* gnuplotでプラネタリウム [#f8934c49]

星空を見るためのプラネタリウムソフトはたくさんありますが、ここでは自分で星空を描いてみたいと思います。

** カタログの入手 [#q70e6026]

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

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

これは Bright Star Catalogue 第5版です。
全天の肉眼で見える程度の星に関する情報が約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)です。

** 全天のプロット [#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 + 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 というファイルができているはずです。

&ref(plot01.png);
トップ   新規 一覧 単語検索 最終更新   ヘルプ   最終更新の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.014 sec.