3月12日は月面モザイク撮影の後、そのままμ-180Cの直焦点でシリウスBを撮影しました。阿南市科学センターが昨年末から「シリウスBチャレンジ」という観測キャンペーンをやっていて触発されたからです。
今年に入ってから公式 Twitter アカウント(@siriusb_wd)もできて、僕はそれで知りました。
シリウスBは2年前に撮っているのですが、今回はセンサーサイズの大きなカメラがあるので周囲の星も写るような写真も撮って、位置や離角を測定してみようと思いました。
機材は月面モザイク撮影の時の状態そのままで、μ-180C 直焦点に ASI294MM Pro (Bin 1)、フィルターなし、という構成です。
まず、周囲の星の位置を記録するために4秒露出で撮ったものがこれです。
縦の高さが 8K 相当(4320 pixel)になるように少しクロップしています。シリウスの光条がなかなかかっこいいです。インスタでもちょっとウケました。まあ、6本の細い光条は副鏡スパイダー由来の回折像だし、8本の太めの光条は方向からしておそらくCMOSセンサー由来のパターンなのでシリウスそのものとは何の関係もない模様ですが…
これを4倍に拡大(クロップ)します。
同じサイズでシリウスBが写るように1/15秒露出で撮って画像処理したものがこれです。
シリウスBに光条が半分かぶさっていますね。2020年に撮った写真ではもう少し離れていたのですが… μ-180C は鏡筒にアリガタが固定されているため、同じ架台に載せれば光条の方向は同じになります。シリウスBは地球から見て時計回りに周っているので来年あたりは光条と重なってしまいそうです。
さて、シリウスBの位置を測定しましょう、ということで、GW中から色々やっていたのですが…
まずは最初の写真を2倍に拡大(クロップ)したものを*1 astrometry.net でプレートソルブをしてみたのですが… ダメです。シリウスの位置のアノテーションが大きくズレてしまいます。
ローカルにインストールした astrometry.net でも試したのですが同様でした。クソ眩しいシリウスが邪魔しているのでしょうか？原因は後で判明するのですが…
これではダメだということで、astropy でなんとか計算する方法を考えます。最初に astrometry.net を使ったのは WCS (天球座標と画像上の座標との変換器)を作るためでした。WCS は astrometry.net がプレートソルブ時に生成する wcs.fits というファイルからロードできます。
調べてみると astropy には基準になる赤経/赤緯がわかっている星の画像上の位置を何点か指定して WCS を生成する方法があるようです。ならばと Stellarium で周囲の星の赤経/赤緯を調べようとしたのですが、シリウスの周りはほとんど星が描かれていません…
シリウスが明るすぎてデータがないのでしょうか？やや離れた位置にいくつか同定できる星があったのでそれを指定して WCS を生成してシリウスAとBのXY座標を天球座標に変換してみたのですが、シリウスAの位置が Wikipedia に書かれた座標とだいぶ違う位置になってしまいました。
まず考えたのは歪曲収差の影響です。でも Dall-Kirkham に歪曲収差ってあるのかな？ちょっと情報がみつけられなかったのですが、コマ収差で星像の重心位置もズレそうなので画面の端の方の星を指定するのは避けたいところ。
それなら真ん中のシリウスAの座標を指定したらよいのでは？ということで、やってみたのですが、シリウスAとBの離角がかなりズレてしまいました。そんなものなのかなー、と思いつつ、赤経/赤緯のグリッド線を描画してみるとなんと赤経線と赤緯線が直交してません。
なにこれ？と思って指定した基準星の座標の入力ミス等がないかチェックしたり、基準星の同定ミスがないかチェックしたりしてみたのですが、問題は見つからず。とにかく基準星を中心近くにみつけようと、Aladin Lite で探しました。
デフォルトの DSS2 等ではシリウスの周りは白飛びしてしまって星が同定できないのですが、2MASS に切り替えると比較的近くまで星が見えるので、これで SIMBAD に載っている星を選んで指定してみたところ、赤経線と赤緯線が直交するようになり、縦横1分角のクロスヘアを描画して縦横比を確認しても問題なさそうでしたが、シリウスAの位置が30秒角近くズレてしまい振り出しに戻りました…
シリウスは地球から近い(8.6光年)ので年周視差の影響？と思ったのですが年周視差は400ミリ秒角以下なので30秒角近いズレは説明がつきませんし、方向的にも合いません。
わけがわからなくてふて寝してたのですが、まさか固有運動？と思って調べたところ年間1.3秒角ぐらい動くようです。SIMBAD に載ってる座標って何時の時点の座標なんでしょう？観測時の座標だとしたらソースが2007年発表の論文なのでそれより前。それともepochがJ2000なので2000年時点の座標なんでしょうか？2000年以後の観測だと観測時点の座標から固有運動の分巻き戻した値？
よくわからないけど、Stellarium は固有運動考慮して位置を表示してるって言うから Stellarium の表示が正しい？と思ったら、固有運動の大きい恒星の座標が間違っているというバグがあるみたいですね…
なんか、1分角以内の誤差ならほとんどのユーザーは気にしないんだから自分で直せとか言われてて "low importance" ラベル貼られて放置されてます…
とりあえずシリウスの固有運動のせいならシリウスAとBの間の位置関係の測定には支障がないはず、と思うことにして、図を描いてみました。
「Sirius A (2007)」が現在 SIMBAD に載っているシリウスA座標で、ソースは VAN LEEUWEN F. Validation of the new Hipparcos reduction. (Astronomy and Astrophysics, volume 474, 653-664 (2007/11-1)) です。
最初に astrometry.net が出した画像を重ねるとアノテーションのシリウスAの位置はこの「Sirius A (2007)」に一致しました。どうやらアノテーションのズレは、astrometry.net が固有運動を考慮していないせいだったようです。
「Sirius A (Stellarium 0.22.1 (2022/3/12))」は Stellarium 0.22.1 で撮影日時の空を表示した時に表示された座標です。「Sirius A (2007)」よりは写真上のシリウスAに近いですが、それでもだいぶズレています。
「Sirius A (observed (2022/3/12))」は写真上のシリウスAの座標です。写真上の位置としては6本の光条の交点を使っています。正確には1点に交わらなかったので、中心部に出来た三角形の中心点を使っています。
「Sirius B (observed (2022/3/12))」は写真上のシリウスBの座標です。光条と少し重なって星像の中心点がわかりづらかったので、反対側に伸びた光条を反転したものを減算して光条のカブリをキャンセルした像を見て中心点を決めました。
シリウスAとBの離角を計算すると約11.6秒角(東に約10.3秒角、北に約5.5秒角)と出ました。シリウスBチャレンジのサイトに載っている値とほぼ一致しました。
観察しよう | シリウスBチャレンジ より。
シリウスBは2021年～2024年頃において、地球から見たとき主星（シリウスA）から最も離れて見える時期を迎えます。このとき離角は11.3秒角に達します。今期を逃せば、次回はまた約50年後となります。なおシリウスBの発見以来、同じような見頃となる時期は過去1870年代、1920年代、1970年代にありました。
観察しよう | シリウスBチャレンジ より。
シリウスAの位置の件は本当に気になりますが、とりあえずそれっぽい値がでてきたのでよしとします。また撮ってみたいのですが、上でも書いたように μ-180C ではスパイダーの光条の関係で2025年頃まで観測は難しそうです。その間シリウスAは固有運動で5秒角くらい動く計算。そのあたりも含めて測定できたら面白いかな…
ちなみにタイトルは「追放系」ラノベタイトルのパロディです。実は「追放系」は一度も読んだことないんですけど… すみません。「なろう」はダンジョン脱出を目指すパーティが全員方向音痴だったという謎な小説しか読んだことがありません。*2
データとスクリプト
最後にこの計算に使ったデータとスクリプトを貼っておきます。ある程度汎用的に作ってありますが動作結果は無保証です。
データ
references.json
座標系の決定に使う基準星のデータです。
[ { "name": "BD-16 1589", "x": 3826, "y": 3398, "ra": "06h44m58.68333s", "dec": "-16d47m50.6410s" }, { "name": "BD-16 1586", "x": 5149, "y": 2276, "ra": "06h44m39.08151s", "dec": "-16d43m43.7368s" }, { "name": "Gaia DR2 2947057475918406784", "x": 2457, "y": 888, "ra": "06h45m19.97874s", "dec": "-16d38m53.3484s" } ]
taegets.json
位置計算の対象となる星のデータです。"ra", "dec" (天球座標)が書いてある星は "x", "y" (写真上の位置)の計算対象、"x", "y" が書いてある星は "ra", "dec" の計算対象です。
座標系が正しく決定されているか確認するために天球座標のわかっているいくつかの星のデータを追加してあります。
[ { "name": "Sirius A (2007)", "ra" : "06h45m08.91728s", "dec" : "-16d42m58.0171s" }, { "name": "Sirius A (Stellarium 0.22.1 (2022/3/12))", "ra" : "06h45m08.45s", "dec" : "-16d43m39.2s" }, { "name": "Sirius A (observed (2022/3/12))", "x": 3240, "y": 2161 }, { "name": "Sirius B (observed (2022/3/12))", "x": 3193, "y": 2135 }, { "name": "TYC 5949-2705-1", "ra" : "06h45m39.85092s", "dec" : "-16d37m27.3042s" }, { "name": "TYC 5949-2765-1", "ra": "06h45m47.85991s", "dec": "-16d41m08.7441s" }, { "name": "CSS 254", "ra": "06h44m28.87639s", "dec": "-16d42m58.3959s" }, { "name": "TYC 5949-2747-1", "ra": "06h44m29.28203s", "dec": "-16d37m20.2212s" }, { "name": "TYC 5949-2643-1", "ra": "06h44m32.59814s", "dec": "-16d35m32.3589s" } ]
astro-image-masurement.py
references.json と targets.json と元画像を与えると計算結果(results.json)と星のマーカーとグリッドを記入したSVG画像を出力するスクリプトです。
#!/usr/bin/env python import json import math import numpy as np from astropy import units as u from astropy.coordinates import Angle from astropy.coordinates import SkyCoord from astropy.wcs import WCS from astropy.wcs.utils import fit_wcs_from_points import svgwrite from PIL import Image import base64 from argparse import ArgumentParser argparser = ArgumentParser(description='get RA/DEC of target position.') argparser.add_argument('ref_json', metavar='reference_stars.json', help="reference stars data file.") argparser.add_argument('target_json', metavar='target_stars.json', help="target stars data file.") # argparser.add_argument('image_width', help="width of image (pixels).") # argparser.add_argument('image_height', help="height of image (pixels).") argparser.add_argument('in_image', help="input jpg/png file.") argparser.add_argument('out_svg', help="output svg file.") args = argparser.parse_args() ref_stars = None target_stars = None with open(args.ref_json, 'r', encoding='utf-8') as f: ref_stars = json.load(f) with open(args.target_json, 'r', encoding='utf-8') as f: target_stars = json.load(f) def ra_dec_to_angle(str_or_deg): if (type(str_or_deg) is str): return Angle(str_or_deg) else: return Angle(str_or_deg, u.deg) ref_ra = [] ref_dec = [] ref_x = [] ref_y = [] proj_point = 'center' for star in ref_stars: ra = ra_dec_to_angle(star['ra']) dec = ra_dec_to_angle(star['dec']) ref_ra.append(ra.deg) ref_dec.append(dec.deg) ref_x.append(star['x']) ref_y.append(star['y']) if ('proj_point' in star) and star['proj_point']: proj_point = SkyCoord(ra=ra.deg, dec=dec.deg, frame="icrs", unit=u.deg) target_xy = [] xy_index = [] target_radec = [] radec_index = [] for i, star in enumerate(target_stars): if 'x' in star: xy_index.append(i) target_xy.append([star['x'], star['y']]) elif 'ra' in star: radec_index.append(i) target_radec.append((ra_dec_to_angle(star['ra']).deg, ra_dec_to_angle(star['dec']).deg)) stars = SkyCoord(ra=ref_ra, dec=ref_dec, frame="icrs", unit=u.deg) pixels_x = np.array(ref_x) pixels_y = np.array(ref_y) wcs = fit_wcs_from_points(xy=[pixels_x, pixels_y], world_coords=stars, proj_point=proj_point) radec_result = wcs.wcs_pix2world(target_xy, 0) # for xy in target_xy: # print(wcs.pixel_to_world(xy[0], xy[1])) for i, star_index in enumerate(xy_index): star = target_stars[star_index] star['ra'] = Angle(radec_result[i][0], u.deg).to_string(unit=u.hour) star['dec'] = Angle(radec_result[i][1], u.deg).to_string(unit=u.degree) xy_result = wcs.wcs_world2pix(target_radec, 0) for i, star_index in enumerate(radec_index): star = target_stars[star_index] star['x'] = xy_result[i][0] star['y'] = xy_result[i][1] output = json.dumps(target_stars, indent=2, ensure_ascii=False) print(output) image = Image.open(args.in_image) w = image.width h = image.height corner_points = wcs.wcs_pix2world([[0, 0], [w, 0], [0, h], [w, h]], 0) ra_min = Angle(min(list(map(lambda p: p[0], corner_points))), u.deg) ra_max = Angle(max(list(map(lambda p: p[0], corner_points))), u.deg) dec_min = Angle(min(list(map(lambda p: p[1], corner_points))), u.deg) dec_max = Angle(max(list(map(lambda p: p[1], corner_points))), u.deg) # print("RA: " + ra_min.to_string(unit=u.hour) # + " - " + ra_max.to_string(unit=u.hour)) # print("DEC: " + dec_min.to_string(unit=u.deg) # + " - " + dec_max.to_string(unit=u.deg)) ra1 = ra_min.hms ra2 = ra_max.hms dec1 = dec_min.dms dec2 = dec_max.dms ra_h_start = int(ra1.h) ra_h_last = int(ra2.h) ra_m_start = int(ra1.m) ra_m_last = int(ra2.m) ra_s_start = int(ra1.s) ra_s_last = int(ra2.s) ra_scale = [] ra_scale_deg = [] for ra_h in range(ra_h_start, ra_h_last + 1): if ra_h == ra_h_last: ra_m_stop = ra_m_last + 1 else: ra_m_stop = 60 for ra_m in range(ra_m_start, ra_m_stop): if (ra_h == ra_h_last) and (ra_m == ra_m_last): ra_s_stop = ra_s_last + 1 else: ra_s_stop = 60 for ra_s in range(ra_s_start, ra_s_stop): ra = (ra_h, ra_m, ra_s) ra_scale.append(ra) ra_scale_deg.append(Angle(ra, unit="hourangle").deg) ra_s_start = 0 ra_m_start = 0 dec_d_start = int(dec1.d) dec_d_last = int(dec2.d) dec_m_start = int(dec1.m) dec_m_last = int(dec2.m) dec_s_start = int(dec1.s) dec_s_last = int(dec2.s) dec_scale = [] dec_scale_deg = [] for dec_d in range(dec_d_start, dec_d_last + 1): if dec_d == dec_d_last: dec_m_stop = dec_m_last + 1 else: dec_m_stop = 60 if (dec_d >= 0) else 1 for dec_m in range(dec_m_start, dec_m_stop): if (dec_d == dec_d_last) and (dec_m == dec_m_last): dec_s_stop = dec_s_last + 1 else: dec_s_stop = 60 if (dec_d >= 0) else 1 for dec_s in list(range(dec_s_start, dec_s_stop)): dec = (dec_d, dec_m, dec_s) dec_scale.append(dec) dec_scale_deg.append(Angle(dec, unit=u.deg).deg) dec_s_start = 0 if (dec_d >= 0) else -59 dec_m_start = 0 if (dec_d >= 0) else -59 # draw output svg drw = svgwrite.Drawing(args.out_svg, size=(w, h)) # draw image with open(args.in_image, 'rb') as f: mime = '' if args.in_image.upper().endswith('.JPG'): mime = 'image/jpeg' elif args.in_image.upper().endswith('.PNG'): mime = 'image/png' image_data = 'data:' + mime + ';base64,' image_data += base64.standard_b64encode(f.read()).decode() drw.add(drw.image(image_data)) # draw grid ra_lines_h = [] ra_lines_m = [] ra_lines_s = [] for i,ra in enumerate(ra_scale): coords = [] for dec_deg in dec_scale_deg: coords.append((ra_scale_deg[i], dec_deg)) line = wcs.world_to_pixel(SkyCoord(coords, frame="icrs", unit=u.deg)) ra_h, ra_m, ra_s = ra # print((ra_h, ra_m, ra_s)) if (ra_m == 0) and (ra_s == 0): ra_lines_h.append(line) elif ra_s == 0: ra_lines_m.append(line) else: ra_lines_s.append(line) dec_lines_d = [] dec_lines_m = [] dec_lines_s = [] for i,dec in enumerate(dec_scale): coords = [] for ra_deg in ra_scale_deg: coords.append((ra_deg, dec_scale_deg[i])) line = wcs.world_to_pixel(SkyCoord(coords, frame="icrs", unit=u.deg)) dec_d, dec_m, dec_s = dec # print((dec_d, dec_m, dec_s)) if (dec_m == 0) and (dec_s == 0): dec_lines_d.append(line) elif dec_s == 0: dec_lines_m.append(line) else: dec_lines_s.append(line) def draw_lines(drw, lines, style): g = drw.g() for line in lines: points = [] for i, x in enumerate(line[0]): points.append((x, line[1][i])) svg_line = drw.polyline(points) svg_line.update({ 'style' : style }) g.add(svg_line) return g grid_g = drw.g() ra_g = drw.g() ra_g.add(draw_lines(drw, ra_lines_h, "stroke: #fcc; stroke-width: 2;")) ra_g.add(draw_lines(drw, ra_lines_m, "stroke: #ffc; stroke-width: 2;")) ra_g.add(draw_lines(drw, ra_lines_s, "stroke: #ccc; stroke-width: 1;")) grid_g.add(ra_g) dec_g = drw.g() dec_g.add(draw_lines(drw, dec_lines_d, "stroke: #fcc; stroke-width: 2;")) dec_g.add(draw_lines(drw, dec_lines_m, "stroke: #ffc; stroke-width: 2;")) dec_g.add(draw_lines(drw, dec_lines_s, "stroke: #ccc; stroke-width: 1;")) grid_g.add(dec_g) drw.add(grid_g) # draw center mark # center = SkyCoord.from_pixel(w/2, h/2, wcs) # min = Angle((0, 1, 0), unit=u.deg) # c0 = center.directional_offset_by(Angle(0, unit=u.deg), min).to_pixel(wcs) # c90 = center.directional_offset_by(Angle(90, unit=u.deg), min).to_pixel(wcs) # c180 = center.directional_offset_by(Angle(180, unit=u.deg), min).to_pixel(wcs) # c270 = center.directional_offset_by(Angle(270, unit=u.deg), min).to_pixel(wcs) def to_xy(pixel): x = (pixel[0]).item() y = (pixel[1]).item() return (x, y) # l1 = drw.line(start=to_xy(c0), end=to_xy(c180), stroke='red', stroke_width='1') # l2 = drw.line(start=to_xy(c90), end=to_xy(c270), stroke='red', stroke_width='1') # drw.add(l1) # drw.add(l2) def draw_markers(stars, color, font_size): for star in stars: x = star['x'] y = star['y'] l1 = drw.line(start=(x - 4, y), end=(x + 4, y), stroke=color, stroke_width='1') l2 = drw.line(start=(x, y - 4), end=(x, y + 4), stroke=color, stroke_width='1') drw.add(l1) drw.add(l2) t1 = drw.text(star['name'], x=[x+24], y=[y-(font_size * 2)], fill=color, font_size=font_size) t2 = drw.text(star['ra'], x=[x+24], y=[y-font_size], fill=color, font_size=font_size) t3 = drw.text(star['dec'], x=[x+24], y=[y], fill=color, font_size=font_size) g = drw.g() g.add(t1) g.add(t2) g.add(t3) drw.add(g) draw_markers(ref_stars, "#8f8", 48) draw_markers(target_stars, "#ff8", 24) drw.save(pretty=True)
results.json
出力結果です。
[ { "name": "Sirius A (2007)", "ra": "06h45m08.91728s", "dec": "-16d42m58.0171s", "x": 3171.2872323943716, "y": 2032.3461331144615 }, { "name": "Sirius A (Stellarium 0.22.1 (2022/3/12))", "ra": "06h45m08.45s", "dec": "-16d43m39.2s", "x": 3198.84000411331, "y": 2223.511470771049 }, { "name": "Sirius A (observed (2022/3/12))", "x": 3240, "y": 2161, "ra": "6h45m07.84766535s", "dec": "-16d43m25.56154398s" }, { "name": "Sirius B (observed (2022/3/12))", "x": 3193, "y": 2135, "ra": "6h45m08.56224723s", "dec": "-16d43m20.10826886s" }, { "name": "TYC 5949-2705-1", "ra": "06h45m39.85092s", "dec": "-16d37m27.3042s", "x": 1143.6019365086354, "y": 469.298876585322 }, { "name": "TYC 5949-2765-1", "ra": "06h45m47.85991s", "dec": "-16d41m08.7441s", "x": 593.4770269367737, "y": 1486.5362595217844 }, { "name": "CSS 254", "ra": "06h44m28.87639s", "dec": "-16d42m58.3959s", "x": 5830.79180786849, "y": 2077.298236424214 }, { "name": "TYC 5949-2747-1", "ra": "06h44m29.28203s", "dec": "-16d37m20.2212s", "x": 5833.743239682289, "y": 511.1390967611835 }, { "name": "TYC 5949-2643-1", "ra": "06h44m32.59814s", "dec": "-16d35m32.3589s", "x": 5622.878222448987, "y": 8.085561764319664 } ]
出力画像の方はこんな感じです。ただし、マーカーのテキストが重なる部分や、緑と黄色以外のマーカーの色は手動で(InkScapeで)修正してあります。
入力画像は4秒露出画像の上に1/15秒露出画像を重ねたものです。1/15秒露出画像に微かに写った周囲の星を使って手動で位置合わせしています。
separation.py
計算結果(results.json)と、二つの星の "name" を指定して離角を計算するスクリプトです。
#!/usr/bin/env python import sys import json from astropy import units as u from astropy.coordinates import Angle from astropy.coordinates import SkyCoord from argparse import ArgumentParser argparser = ArgumentParser(description='get RA/DEC of target position.') argparser.add_argument('stars_json', metavar='stars.json', help="stars data file.") argparser.add_argument('star_name_a', help="name of star A.") argparser.add_argument('star_name_b', help="name of star B.") args = argparser.parse_args() stars = None with open(args.stars_json, 'r', encoding='utf-8') as f: stars = json.load(f) star_a = None star_b = None for star in stars: if star['name'] == args.star_name_a: star_a = star elif star['name'] == args.star_name_b: star_b = star coord_a = SkyCoord(Angle(star_a['ra']).deg, Angle(star_a['dec']).deg, unit=u.deg) coord_b = SkyCoord(Angle(star_b['ra']).deg, Angle(star_b['dec']).deg, unit=u.deg) for_ra = SkyCoord(Angle(star_b['ra']).deg, Angle(star_a['dec']).deg, unit=u.deg) for_dec = SkyCoord(Angle(star_a['ra']).deg, Angle(star_b['dec']).deg, unit=u.deg) print(star_a['name'] + ": " + star_a['ra'] + " / " + star_a['dec']) print(star_b['name'] + ": " + star_b['ra'] + " / " + star_b['dec']) print("separation: " + coord_b.separation(coord_a).to_string(unit=u.deg) + " (" + coord_a.separation(for_ra).to_string(unit=u.deg) + " / " + coord_a.separation(for_dec).to_string(unit=u.deg) + ")")
results.json、"Sirius A (observed (2022/3/12))"、"Sirius B (observed (2022/3/12))" を指定し以下の出力を得ました。
Sirius A (observed (2022/3/12)): 6h45m07.84766535s / -16d43m25.56154398s Sirius B (observed (2022/3/12)): 6h45m08.56224723s / -16d43m20.10826886s separation: 0d00m11.62396968s (0d00m10.26536042s / 0d00m05.45327512s)
