Hip-HR_diagram

ヒッパルコス＆チコカタログを用いた天文教材作成備忘録　2008　Apr.04
（本稿の内容はもちろん無保証です．間違いなどがあればご連絡ください）

まずHipparcosデータ本体は国立天文台のftpサイト
ftp://dbc.nao.ac.jp/DBC/NASAADC/catalogs/1/1239
配下に
hip_main.dat.gz
として置いてある．これを展開すると

H| 1| |00 00 00.22|+01 05 20.4| 9.10| |H|000.00091185|+01.08901332| | 3.54| -5.20| -1.88| 1.32| 0.74| 1.39| 1.36| 0.81| 0.32|-0.07|-0.11|-0.24| 0.09|-0.01| 0.10|-0.01| 0.01| 0.34| 0| 0.74| 1| 9.643|0.020| 9.130|0.019| | 0.482|0.025|T|0.55|0.03|L| | 9.2043|0.0020|0.017| 87| | 9.17| 9.24| | | | | | | | 1| | | | | | | | | |S| | |224700|B+00 5077 | | |0.66|F5 |S

という恐ろしい項数のデータとなる（これで星１個分，449カラムで１行）．このデータ諸元は
ftp://dbc.nao.ac.jp/DBC/NASAADC/catalogs/1/1239/ReadMe
に詳細がある．

そこでまずここから，必要な項のみを抜き出す処理をする．

基本的に色付けされているのがUNIXコマンド
まず
tr ' ' '_' < hip_main.dat> h_m2.txt
で空白を_に変更（こうしないとawkのsubstrで文字数に応じた切り出しができない．FORTRANならread文，BASICなら mid$関数を用いるところ．Cは文字列操作が難しいのでお勧めできない．ここはawkの一人舞台！）

H|___________1|_|00_00_00.22|+01_05_20.4|_9.10|_|H|000.00091185|+01.08901332|_|___3.54|___-5.20|___-1.88|__1.32|__0.74|__1.39|__1.36|__0.81|_0.32|-0.07|-0.11|-0.24|_0.09|-0.01|_0.10|-0.01|_0.01|_0.34|__0|_0.74|_____1|_9.643|0.020|_9.130|0.019|_|_0.482|0.025|T|0.55|0.03|L|_|_9.2043|0.0020|0.017|_87|_|_9.17|_9.24|_______|_|_|_|__________|_|__|_1|_|_|_|__|___|_______|_____|_____|____|S|_|_|224700|B+00_5077_|__________|__________|0.66|F5__________|S_
H|___________2|_|00_00_00.91|-19_29_55.8|_9.27|_|G|000.00379737|-19.49883745|+|__21.90|__181.21|___-0.93|__1.28|__0.70|__3.10|__1.74|__0.92|_0.12|-0.14|-0.24|-0.29|_0.01|_0.21|-0.02|-0.19|-0.28|_0.14|__2|_1.45|_____2|10.519|0.033|_9.378|0.021|_|_0.999|0.002|G|1.04|0.00|I|_|_9.4017|0.0017|0.015|120|_|_9.37|_9.44|_______|C|_|_|__________|_|__|_1|O|_|_|__|___|_______|_____|_____|____|_|_|_|224690|B-20_6688_|__________|__________|1.04|K3V_________|4_

次に各パートを抽出．（急ぐ方は下記のバグバージョンを読み飛ばしてください）

ーーまずバグバージョン－－－
　　　　      HID　　　　      赤経　　　　　　赤緯　　　　　　　実視等級　　    年周視差　　　年周視差(error) B-V　B-V(error) 　     V-I
awk '{print substr($0,9,6),substr($0,18,11),substr($0,30,11),substr($0,42,5),substr($0,80,7),substr($0,120,6),substr($0,246,6),substr($0,253,5),substr($0,261,4)}' h_m2.txt > h_m3.txt

これで
#HID　　赤経　　　　　赤緯   実視等級年周視差 error B-V　error　V-I
_____1 00_00_00.22 +01_05_20.4 _9.10 ___3.54 __1.39 _0.482 0.025 0.55
_____2 00_00_00.91 -19_29_55.8 _9.27 __21.90 __3.10 _0.999 0.002 1.04
_____3 00_00_01.20 +38_51_33.4 _6.61 ___2.81 __0.63 -0.019 0.004 0.00
ができる．この_をはずすために
tr '_' ' ' < h_m3.txt> h_m4.txt
#HID　　赤経　　　　　赤緯   実視等級年周視差　SE 　B-V　　SE　V-I
     1 00 00 00.22 +01 05 20.4 9.10    3.54   1.39 0.482 0.025 0.55
     2 00 00 00.91 -19 29 55.8 9.27   21.90   3.10 0.999 0.002 1.04
     3 00 00 01.20 +38 51 33.4 6.61    2.81   0.63 -0.019 0.004 0.00
ができる．これでは赤経，赤緯が扱いにくいので，これを小数に直すがここで失敗．昔恥ずかしいバグを作ったことがあるので要注意．つまり赤緯の-は度の桁にしか及んでいないことに注意．
awk '{print $1, $2+$3/60+$4/3600, $5+$6/60+$7/3600,$8,$9,$10,$11,$12,$13}' h_m4.txt > h_m5.txtは
駄目！
ーーバグバージョンおわり－－－

各パートを抽出からやり直し．（符号を別項に抽出する）
#　　　　　　　HID　　　　　　赤経　　　　　　　赤緯（符号）　　赤緯　　　　　　実視等級　　　年周視差　　　　　SE 　　　　　　B-V　　　　　SE　　　　V-I
awk '{ ( if substr($0,249,1) != "_" ) print substr($0,9,6),substr($0,18,11),substr($0,30,1),substr($0,31,10),substr($0,42,5),substr($0,80,7),substr($0,120,6),substr($0,246,6),substr($0,253,5),substr($0,261,4)}' h_m2.txt > h_m3.txt
_____1 00_00_00.22 + 01_05_20.4 _9.10 ___3.54 __1.39 _0.482 0.025 0.55
_____2 00_00_00.91 - 19_29_55.8 _9.27 __21.90 __3.10 _0.999 0.002 1.04
_____3 00_00_01.20 + 38_51_33.4 _6.61 ___2.81 __0.63 -0.019 0.004 0.00
_____4 00_00_02.01 - 51_53_36.8 _8.06 ___7.75 __0.97 _0.370 0.009 0.43
_____5 00_00_02.39 - 40_35_28.4 _8.55 ___2.87 __1.11 _0.902 0.013 0.90
_____6 00_00_04.35 + 03_56_47.4 12.31 __18.80 __4.99 _1.336 0.020 1.55

(上のコマンドの( if substr($0,249,1) != "_" ) の部分がないと，ところどころにB-Vのデータがない行ができるのでそれをはじくための条件設定)

再び_を空白に戻す．
tr '_' ' ' < h_m3.txt> h_m4.txt
#HID　　赤経　　　　　赤緯   実視等級年周視差　SE 　B-V　　SE　V-I
     1 00 00 00.22 + 01 05 20.4 9.10    3.54   1.39 0.482 0.025 0.55
     2 00 00 00.91 - 19 29 55.8 9.27   21.90   3.10 0.999 0.002 1.04
     3 00 00 01.20 + 38 51 33.4 6.61    2.81   0.63 -0.019 0.004 0.00
     4 00 00 02.01 - 51 53 36.8 8.06    7.75   0.97 0.370 0.009 0.43
     5 00 00 02.39 - 40 35 28.4 8.55    2.87   1.11 0.902 0.013 0.90
     6 00 00 04.35 + 03 56 47.4 12.31   18.80   4.99 1.336 0.020 1.55

次に赤経，赤緯を小数に直す．
awk '{print $1, $2+$3/60+$4/3600, $5($6+$7/60+$8/3600),$9,$10,$11,$12,$13,$14}' h_m4.txt > h_m5.txt
やっとこれで
#HID　赤経　　赤緯実視等級年周視差 SE B-V SE V-I
1 6.11111e-05 +1.089 9.10 3.54 1.39 0.482 0.025 0.55
2 0.000252778 -19.4988 9.27 21.90 3.10 0.999 0.002 1.04
3 0.000333333 +38.8593 6.61 2.81 0.63 -0.019 0.004 0.00
4 0.000558333 -51.8936 8.06 7.75 0.97 0.370 0.009 0.43
5 0.000663889 -40.5912 8.55 2.87 1.11 0.902 0.013 0.90
6 0.00120833 +3.9465 12.31 18.80 4.99 1.336 0.020 1.55
という描画用の最終データ h_m5.txtが完成．

ここまでの最終データ処理スクリプトは
#!/bin/bash
#　Hipparcosデータから恒星を抽出するスクリプト
#　by Y.Okamoto 27/March 2008
#HID　赤経　　赤緯実視等級年周視差 SE B-V SE V-I
#----------------------------------------------------------------------------------------------
# 0.000663889 -40.5912 8.55 2.87 1.11 0.902 0.013 0.90
tr ' ' '_' < hip_main.dat> h_m2.txt
awk '{ ( if substr($0,249,1) != "_" ) print substr($0,9,6),substr($0,18,11),substr($0,30,1),substr($0,31,10),substr($0,42,5),substr($0,80,7),substr($0,120,6),substr($0,246,6),substr($0,253,5),substr($0,261,4)}' h_m2.txt > h_m3.txt
tr '_' ' ' < h_m3.txt> h_m4.txt
awk '{print $1, $2+$3/60+$4/3600, $5($6+$7/60+$8/3600),$9,$10,$11,$12,$13,$14}' h_m4.txt > h_m5.txt

これを用いてもよいが，12万行近いデータで精度の悪い測定の星も含まれるのでこれを精度で間引く．年周視差（Parallax）とB-Vについて，その誤差（stand errorと表記）の割合を求め，求める誤差未満のもののみを選び出す．
awk '{ if ( ($6/$5 < 0.10) && ($8/$7 < 0.25) ) print $1,$2,$3,$4,$5,$6,$7,$8,$9 }' h_m5.txt > h_m6.txt
でdistanceで10%，色で25%の星が選択できる．これをh_m6.txtとする．星の総数２４１３４個

次にdistanceで5%，色で15%の星の選択は
awk '{ if ( ($6/$5 < 0.05) && ($8/$7 < 0.15) ) print $1,$2,$3,$4,$5,$6,$7,$8,$9 }' h_m5.txt > h_m7.txt
これをh_m7.txtとする．星の総数９２６３個．しかしこれは少し星が少なすぎて却下．

最終的な白黒HR図を描画する基本シェルスクリプトは以下のとおり（このスクリプトには誤差選択の上記部分を含んでいる）
#!/bin/bash
#　GMTでHR図を描くスクリプト
#　by Y.Okamoto 27/March 2008
# Hipparcos data
# #HID RA[deg] Dec[deg]V[mag]Parallax SE B-V SE   V-I
#----------------------------------------------------------------------------------------------
# 2 0.000252778 -19.4988 9.27 21.90 3.10 0.999 0.002 1.04
# $1 $2           $3     $4   $5    $6   $7    $8    $9
export PATH=$PATH:/usr/lib/gmt/bin
awk '{ if ( ($6/$5 < 0.10) && ($8/$7 < 0.25) ) print $7,$4-5.0*0.434294482*log(1000/$5)+5.0 }' h_m5.txt | psxy -Jx5.0c/-0.7c -R-0.5/2.0/-6/15 -Sc0.02c -Ba1.0f1.0g1.0> HR2.ps

これで作られる図は下記のとおり．あとはpsxyの各項を必要に応じて適当に直せばよい．