はやぶさ2の地球スイングバイの見え方をFortranで計算していた話
はやぶさ2の地球スイングバイについて
小惑星探査機「はやぶさ2」の地球スイングバイ実施結果について
「はやぶさ2」は、平成27(2015)年12月3日(木)の夕方から夜にかけて地球スイングバイを実施し、19時08分(日本時間)に地球に最接近、ハワイ諸島付近の太平洋上空約3,090kmを通過しました。
はやぶさ2は,近傍の小惑星「リュウグウ」ヘ向け,2014年12月3日に種子島宇宙センターからH-IIAロケット26号機で打ち上げられたものでした.
打ち上げ後,しばらくは地球軌道と同じような公転軌道を回っていましたが,打ち上げから約1年後,地球の引力を利用した地球スイングバイが実施されています.
これは,地球の引力により,はやぶさ2にさらなる加速をつけ,軌道も同時に変更するといったものです.
イメージ的には,磁石があるとき,その磁石に対し,パチンコ玉などの鉄球をある角度,速度でうまく入射してやると,磁石の引力によりパチンコ玉の軌道が変化し,同時にパチンコ玉が加速すると言ったような感じです.これは磁石で遊んだことがある人ならわかるような気がします.
イメージ的にはこんな感じ.
このしくみと全く同じようなことを,宇宙,地球という非常に大きなスケールで行うのが地球スイングバイです.
この利点は,うまく入射軌道を調節してやれば,燃料をあまり消費せずに大きな加速度が得られ,速度が上がるということです.
これは長い旅をするはやぶさ2にとっては非常に重要なことになります.
はやぶさ2の地球スイングバイを見る
実は,このときのはやぶさ2の地球スイングバイは,ちょうどハワイから日本の上空,ロシアにかけて地球に近づき,太平洋上空ではやぶさ2が地球に最接近するものでした.
ということは,そうです.日本からスイングバイが見れる可能性があったのです.
そして,実際そのとき,複数の天文台ではやぶさ2が見れたという報告が上がっておりました.
私も,はやぶさ2スイングバイを見れるという情報を聞き,トライをしてみることを決めました.
後から考えると,はやぶさ2は大型望遠鏡を持つ天文台のみでしか観測できなかったという無謀な挑戦だったのですが,
実はその当時は,はやぶさが何等級の光として見えるか定かではなかったのです.
太陽電池パドルの細かな向きによっても何等で見えるかが全く変わってくるためです.これははやぶさの姿勢をきちんと計算しても誤差要因によってかき消されてしまったり,日本のどこで見るかで全く違うことになり得るため,明るさは"見てみないとわからない"という状況でした.この希望にかけて私も観測を試みました.ミラーレスカメラを持ち,露光モードを長秒用の撮影にし,ずっとシャッターを切っていましたが,できた写真にははやぶさ2らしき軌跡は確認できませんでした.
はやぶさの軌道からどう見えるのか計算したくなった
観測に失敗した私は,実際自分の位置からはやぶさがどう見えるのか気になりました. 地球の緯度,経度からはやぶさ2スイングバイの見え方を計算してくれたサイトがありました.
2015年12月3日の「はやぶさ2」地球スイングバイの観測について(参考情報 その1)
このサイトではこのような軌道が計算されて公表されています.
(上記サイトより引用,ステラナビゲータ10/(株)アストロアーツ)
私はこれを参考に観測に挑戦していましたが,実はJAXAからはやぶさ2のスイングバイに関する軌道データが公表されています.
トピックス 「はやぶさ2」の地球スイングバイについての情報を公開します!
ここの軌道データは3種類
データ1(UTF-8):hy2_trj_EC_10min.txt
・黄道座標系
・11月3日から1月3日まで、10分間隔
・地球およびRyuguのデータ付き
データ2(UTF-8):hy2_trj_EQ_1min.txt
・赤道座標系
・12月1日から12月5日まで、1分間隔
・「はやぶさ2」のみ
データ2(UTF-8):hy2_trj_EQ_60min.txt
・赤道座標系
・11月3日から1月3日まで、60分間隔
・「はやぶさ2」のみ
この座標データから自分で軌道を計算して見ることにしました.
軌道を計算しよう
軌道を計算するのは,簡単といえば簡単です. ただの座標変換だからです. JAXAからは,ある座標系におけるはやぶさ2の位置と速度がxyz方向で示されています. したがってこの位置を座標変換し,自分のいる緯度,及び経度を参考に地球上での見え方にマッピングしてあげればよいのです. 私は今回,データ2を用いることにしました.
このデータは,赤道座標系,これは地球中心J2000赤道面基準慣性座標系と呼ばれています.
これは下記の富山大学の下記webサイトによると.
衛星航法システムにおける座標系
J2000.0(力学時2000年 1月1.5日,JD2451545.0)における 平均赤道面,平均春分点の方向を基準(平均春分点の方向がx軸) とする慣性座標系であり、衛星運動方程式を積分して衛星の軌道 を生成する座標系である.
とされています.つまり2000年 1月1.5日の天文の状況を考えてそれを基準にしますということです.
日々刻々と春分点や赤道面は地球の歳差運動などにより変化しうるので,これを日付も含めてfixしているということです.
この座標系から,地球の地上の座標に変換するだけです.
これをfortranで書いていましたので,今更ですが公開してみます.
PROGRAM haya2_swingby implicit none double precision coordinates(3),ecef_hy(3) integer year,month,day,hour, min,sec,i double precision, parameter :: longitude=140.872183d0, latitude=38.268737d0,height=0.0d0 double precision, parameter :: a=6378137.0000d0, b=6356752.314250d0, omega_e=-7292115.467d-11 double precision, parameter :: epoch_j2000=1449129600-946728000 double precision, parameter :: Pi=3.141592653589793238 double precision N_o,h_o,x_o,y_o,z_o double precision x,y,z double precision ro,azi,ele open(20, file="./hy2_trj_EQ_1min_mod.txt", status="old") DO i=0,959 read(20, '(I4,5(1X,I2),3F20.1)'), year, month, day, hour, min, sec, coordinates(1), coordinates(2), coordinates(3) ecef_hy(1) = coordinates(1) * 1000 * cos(omega_e*(epoch_j2000+i*60)) + coordinates(2) & * 1000 * sin(omega_e*(epoch_j2000+i*60)) ecef_hy(2) = coordinates(1) * 1000 * (-1)*sin(omega_e*(epoch_j2000+i*60)) + coordinates(2) & * 1000 * cos(omega_e*(epoch_j2000+i*60)) ecef_hy(3) = coordinates(3) * 1000 N_o=a**2/dsqrt(a**2*dcos(latitude*Pi/180)**2+b**2*dsin(latitude*Pi/180)**2) x_o=(N_o+height)*dcos(latitude*Pi/180)*dcos(longitude*Pi/180) y_o=(N_o+height)*dcos(latitude*Pi/180)*dsin(longitude*Pi/180) z_o=(b**2/a**2*N_o+height)*dsin(latitude*Pi/180) !print *, year, month, day, hour, min, sec, coordinates(1), coordinates(2), coordinates(3) x=ecef_hy(1)-x_o y=ecef_hy(2)-y_o z=ecef_hy(3)-z_o !print *,x,y,z ro = sqrt(x**2+y**2+z**2) azi = atan2(z,sqrt(x**2+y**2)) * 180 / Pi ele = (atan2(y,x)) * 180 / Pi print '(5I4,2F10.2,F15.1)',year,month,day,hour,min,azi, ele, ro/1000 write(10,*) ele,azi END DO print *,epoch_j2000,4*3600+42*60+5 OPEN(13,file="gnupdummy.plt") write(13,*) "reset" write(13,*) "set xrange [-180:180]" write(13,*) "set yrange [-90:90]" write(13,*) "plot 'fort.10' with points ps 0.5 pt 1"!with lines lw 1" CLOSE(13) call system("gnuplot gnupdummy.plt") END PROGRAM haya2_swingby
緯度: 38度16分7.454秒(38.268737) 経度: 140度52分19.86秒(140.872183)
は私の住んでいる宮城県の宮城県庁の座標です.これを使ってスイングバイ付近のデータで計算してみると
このようなデータが得られます.
横軸は-180度から180度の水平角度,縦軸は-90度から90度の仰角を表しています.
仰角がマイナスということは水平線の下にいるということです.
先ほどのwebサイトのデータを全く同じということにはなりませんでした.
これは,計算の方法で座標系の近似などを行ってしまったためだと思っています.
きちんと厳密に計算しているわけではないのです.
厳密にやりなおしてみたいと思いつつ今までやっていません.
もしかすると,完全に計算を間違えているのかもしれません.わかりません.
もし知見がある方がおられましたら是非ご指摘いただけると幸いです.
とはいえ,間違っているかもしれませんが,だいたいの傾向は同じものとなっています.
一応自分で参考程度の軌道情報を計算することができました.
当時,他の複数のwebサイトで検証をしていたのですが,この軌道に似ている形のものもあればそうでないものもあり,
計算する人が使う手法によっても軌道の形は変わってしまうようです.
天文,宇宙の世界は細かい世界だと実感した瞬間でした.
References
2017年になりました
新年のごあいさつ
2017年となりました.新年明けましておめでとうございます.
2017年1月1日のうるう秒挿入のせいで時差ボケしています
更新する頻度が少ない本ブログですが,ちまちまと来てくださる方も増え,嬉しく思っています.
2016年で印象的だった記事は,やはりMP640の修理の記事です.
Canon MP640プリンタ給紙機構の修理
Canon MP640のカムをDMM.makeに出品してみた話
DMM.makeに出品したCanon MP640のカムを用いて修理してくれた方がいた
ふと「修理しました!」という報告程度に上げてみた記事から,「私も修理したい」という声が上がり,実際に修理ができてしまったのはこちらとしてもびっくりで嬉しかった.
「壊れたから捨てる」ではなく,「壊れたけどなんとかして使えないだろうか」という考えをしてしまうのが私なのですが,これは私がものづくりをやっていたことから上がってくる考えなのでしょうか.壊れたものをできる限り直して使うというのが,ものづくりをする人間として,そのものを作ってくれた人に対する感謝ではないかと心の片隅で思っていたりします.
とはいえ,不要になったものは捨てざるを得ないわけですが...できる限り,です.
要するに簡単に同じような機能の新しいもので"交換"したくないという思いな気がします.新しい機能とかがついて進化したものは欲しくなって買ってしまったりします.私は.
なんというか,機能や性能に不満はないけど不具合が生じてしまって使えなくなり,同じようなものを買いたくない感じ.
機能とか性能が劣るとか,「こっちのほうがいいな」って思ってしまうようだとやっぱり買いたくなっちゃう.
私はNHKの「仕事の流儀」という番組が好きでよく見ています.あのプロフェッショナルのアレです.
その中でこんなプロフェッショナルがいまして
時計に命、意地の指先 - 松浦敬一
この人は動かなくなった時計のうち,他のお店が修理を断るような時計でも,直してしまうようなプロだそうです.
そしてこの人は,古くからあるものをなんでも修理して直して使っているそうです.この番組の中では掃除機や電灯がそうでした.
さすが修理職人といった感じ.それで私は,壊れたものを修理して使うところに共感を覚えたりしたわけでした.
小さい頃から修理とか分解(元に戻せなくなることもしばしば)が好きだった私なのですが,修理とか分解をすると本当に中身がよくわかったりするものです.
そこから学ぶことは山のようにあり,ネジの使い方,止め方,止める位置から熱設計,振動設計(耐久性),ケーブリングなどなど,一流の企業が作ったモノであれば,一流のテクニックがモノに詰まっているはずです.そのモノを見るということは,一流の技術を見てることに他ならない気がします.分解するだけで一流の技術が見れるなんて素晴らしくないですか?
壊れて捨てることが決定しても,分解して勉強してから捨てるようにしています.それが私なりの供養です.そこから技術を学ぶことが私のポリシーです.
新年始まって早々,「分解はいいぞ」という記事になってしまった.
まだまだ書きたいことがあるのでちょくちょく書いていきたい.
まとまりがありませんが,「分解はいいぞ」ということを言って,新年のご挨拶にさせていただきます.
本年もよろしくお願いいたします.
WSJT系の備忘録
いろいろリファレンスを読んだので,「へー」と思ったことを書いてみる.
すでに知っていた部分は書かない.
WSJT-X
TUNEボタン: TXしてJT65等の単一トーンを発生させる.クリック音やグリッチ音が含まれないか,綺麗な信号かを確認する.
PCのミキサを調整して,TX POWERが10%程度低下する程度に出力レベル調整すると良い.アンテナチューナ等の調整にも使える.
送信に関して,受信と送信にはもちろんフィルタがあり,受信に関してはバンド幅を設定できるものが多い.
できるだけ幅広い周波数を受信てできるようにしておく.最低4kHz程度あるとよい.
ただ,送信に関してはフィルタされているものが多い(2700Hz以上あたりで).TXフィルタの帯域をいじれればよいが,そうでないものなどには
Splitモードを用いると良い.WSJT-X v1.1からはVFOAとVFOBを用い,VFOBの周波数を送信する場所に合わせて自動的に調整するとことで,
TXのオーディオ信号が1kHzから2kHzの中に入るようになっている.CATとSplit TXをONにしておくことにより使うことができる
WSJT v10.0
Tol: frequency tolerance, デコード対象の周波数範囲
DT: Diff. Timeかな?
DF: Diff. Freq.かな.
AFC: Auto Freq. Controlとかでしょうか.デコードアルゴリズムを走らせる時に,周波数のドップラーシフト等に追従するようにしてくれる.
Clip: 受信した信号をデコードする前にクリップする.通常は0で,スタティックノイズなどがある場合はこの値を増加させるとクリップしてくれる.
Exclude: Ave. Messageのから,最近のメッセージをExcludeする.DFやDTが著しくおかしいデータが入ってしまった時に前までのデータが汚染することを防ぐ.
Freeze: Enableにし,赤のスパイク波をクリックする(周波数を決める)ことで,DFを考慮したTol[Hz]幅の周波数のみをデコード対象とする.
知見みたいなものとか.
メイン画面のスペクトル表示の赤いラインは周波数偏差(SYNC TONE)を,青いのはDTが何秒かを示し,緑は謎.
緑は,平坦だったらといあえずOKらしい.バックグラウンドのノイズを示すらしい.変にうねってたりするとノイズがリズミカルではないということらしい.
一度同期が取れたら,赤い波形(スパイク波)をクリックしてFreezeし,Tolを小さくしておく(100Hzもしくはそれ以下)
そうするとその時点のDF値をもとに解析をしてくれる.
JT65のショートハンドメッセージのテキスト欄に「?」が表示される原因として,赤のSYNC信号を選択してFreezeされていないか,Tolが100以下に設定されていない.
または,メッセージテキストが全て解読されていないなどの原因がある.
WSJTは平均化処理によるデコード機能がすごいらしい.
メイン画面で見えてなくても,複数回同じものを受信したものを相加平均をとることでS/Nをあげて,デコードする.
これのためにExcludeボタンとかがあるのか.なるほど.
ディープサーチ機能もある.
これはCALL3.TXTに記載されているコールサインとの比較をして,デコードを補助する的なもの.
ノーマルとアグレッシブとアグレッシブ+平均化データを用いるというモードがある.
ディープサーチも失敗することがあるので,本当にあっていそうなデータか疑って確かめる目が必要らしい.
送信開始から10秒以内であれば,送信メッセージを変更しても受信が成功することが多いらしい.
EME, JT65備忘録
ただの備忘録です.
referenceにあるページから引っ張ってきたやつを一部引用させていただいています.
JT65関連
JT65は流星痕反射通信用に開発されたデジタル通信技術を、ノーベル物理学賞受賞のK1JT/Dr.Joe TaylorがEME(月面反射通信)などの微弱信号通信を目的として改良・開発した狭帯域通信モードであるWSJTのうちの一つ.
JT65には3つモードがある.
JT65AはHF/50MHz用、JT65Bは144MHzと432MHz用mJT65Cは1296MHz用となっているらしい.
JT65プロトコルは65-tone(65個のトーン)を使い,1270.5HzのサブキャリアをAFSK変調する(65FSK).この変調方式はon-offキーイングよりも非常に効率が良く,PSKよりも周波数変動に対する許容度がはるかに高い特長がある.
ソフトウェアはとりあえず現在はWSJT-Xを使っていればよさそう.
JT65A(HF)運用周波数
JT65Aの国際的な中心運用周波数は以下.モードはUSB.
1.8MHz帯 : 1838kHz
3.5MHz帯 : 3576kHz
7MHz帯 : 7076kHz(但し、日本ではバンドプランによって受 信のみ、代わりに7039kHz)
10MHz帯 : 10139kHz,10147kHz
14MHz帯 : 14076kHz
18MHz帯 : 18102kHz,18106kHz
21MHz帯 : 21076kHz
24MHz帯 : 24092kHz
28MHz帯 : 28076kHz
50MHz帯 : 50300~51000kHz
EME関連.
QRPpな局は,グランドエフェクトをうまく利用するとよい.
これは大地反射により利得が大きく見える.最大で+6dB程度利得増加が見込めるらしい.
月の出と月の入の際に,月が地面に近づくときがチャンス.
月の仰角が10度前後までならグランドエフェクトを利用するために仰角は水平にしておくと良いらしい.
グランドエフェクトは仰角20度くらいまでなら使えるらしい.
月のコンディションに関しては,DGRDという評価項目がある.
これは月の背景ノイズ,月との距離などから算出されるらしい.
算出された結果によるコンディションがカレンダーになっている.
EA6VQ EME calendar
DGRD が -5.0dB から 0dB の間がターゲットになる.欧州の QRP 局は -2.5dB を切ったあたりからアクティブになる傾向があるらしい.
JT65Bのデコード限界値はS/N -32dBほどらしいので,相手局の出力やアンテナを考慮して,受信レベルから相手が受信すると想定されれるレベルを計算して-32dBを下回っていないか確認すると良いかも.
偏波面問題
EMEでは偏波面が変わってしまうことがある.
変わる要因は主に2つあり,Geometric RotationとFaraday Rotationである.
Geometric Rotation が |0-20| または |70-90| くらいの範囲で, Faraday Rotation がおとなしい時にだけお互いの偏波面は保存される.
ファラデー回転は電磁気学の法則であり記録媒体であるMOディスクとかに使われている技術ですよね.今はもう使われなくなってきましたが.
Geometric Rotationは月とお互いの位置によるものらしい.WSJTにおいては,Dpolという表示があり,これがGeometric Rotationを示しているらしい?
これが0の時は受信と送信は同じ偏波を使えば良いが,45の時は送信と受信で偏波面を変えてやるとよいとのこと.
EME関連チャットツール等
R言語を使ってアンテナの指向性パターンをグラフ化する
電磁界と戯れる場合,アンテナと戯れることになる.
アンテナを議論する時に,主な評価項目として,VSWRと利得,ビームパターン(指向性)が議論されることであろうと思う.
VSWRと利得は数値として得ることができるから良いものの,ビームパターンは測定には電波暗室が必要であるし,個人でやる上ではビームパターンを測定することすら難しいであろう.
仮に測定できた場合や,シミュレーションで角度と利得の数値が得られたとしても,グラフ化するソフトがあまりないのだ.
今回はこのグラフ化に関して書いていきます.
フリーで指向性パターンを描画できる簡単なソフトとして,plot32 (Plots32の詳細情報 : Vector ソフトを探す!)や,AP Chart(APchart Homepage)というソフトがある.
しかしながらこれらのソフトはイマイチ勝手が効かなかったりすることも多い.
そこで今回は統計解析言語であるR言語のグラフプロット機能を利用してみることにする.
R言語による指向性パターンの描画には protrixパッケージを用いて,レーダーチャートを利用して描画を行う.
レーダーチャートとはその名の通りレーダーのような図であり,一般に複数の項目のバランスなどを議論する際によく用いられます.
ゲームとかでよくあるこんなやつですね.
レーダーチャートをベースとして指向性パターンを描いていきたいと思います.
必要な情報は,アンテナを回したときの角度と,そのときの利得です.
利得でなくても,受信電界強度などでもよいです.せっかくR言語を使っているのですから,R言語にフリスの伝達公式のとおりに計算して貰えばいい話です.
まずはplotrixをロードします.
# パッケージの読み込み library(plotrix)
これでplotrixが使えるようになりました.
パッケージがそもそもインストールされていない場合は,
install.packages("plotrix")
でパッケージをインストールします.
アンテナの回転角と利得のデータは,CSVやタブ区切りなどで保存されている場合が多いかと思います.
回転角とそれに対応する利得のデータをロードします.
# ファイルのロード(必要に応じてread.csv()などを用いればいいと思います) hplane <- read.table("ファイルまでのパス") #余計なカラムがある場合は,必要なカラムのみにしておくと見やすいと思います. #たとえば1列目と9列目のみを抽出 hplane <- hplane[,c(1, 9)]
ここからはフリスの伝達公式の計算です.もともと利得がわかっている場合は計算する必要はありません.
電波暗室などで測定したとして,送信電力100mW(20dBm), 送信アンテナのゲインが真数で10, アンテナ間の距離が3m,波長が0.1223642686 m(2.45GHz)と仮定してみます.
hplaneの2列目に受信電界強度[mW]が入っているとすると,下記によって一気に計算してくれます.
send_power <- 100 sender_gain <- 10 distance <- 3 wavelength <- 0.1223642686 calced_hplane <- 10*log10((10^(hplane[,2]/10))/send_power/sender_gain*(4*pi*distance/wavelength)^2)
calced_hplaneには利得[dB]が入っていることになります.
さらに,角度軸の準備をします.
何度のところに線を引いておきたいかによってこの値をいじります.
45度ずつ,上が0度ならばこんな感じです.
azimath <- c(180, -135, -90, -45, 0, 45, 90, 135)
これでデータの準備はできました.
描画をします.たとえばこんな感じに引数を渡します.
radial.plot(calcd_hplane, labels=azimath, rp.type="p",start=pi/2,clockwise=TRUE, radial.lim = c(-40,20),radial.label = "", lwd=2)
凡例はたとえば以下のように,
legend("topright", legend = c("H-Plane", "E-Plane"), lty = c(1,3), lwd = 3)
これで電磁界解析ソフトを利用したシミュレーション結果をR言語とplotrixを用いて描画してみた結果が以下のようになります.
わりと綺麗に描画できていると思います.
引数はたくさんあるので,リファレンスを参照するとよいでしょう.
たとえばですが,
radial.plot( [描画するデータ], labels=[ラベルデータ], rp.type="p", #r:radial lines, p:polygon start=-pi/2, #描画スタート位置のオフセット,デフォルトは3時の方向 clockwise=TRUE, #時計回りに描画 radial.lim = c(-20,15), #ここでいうゲイン軸の範囲,つまり円に垂直な軸の範囲 lty = c(1,2), #line type, 直線,破線,いろいろあります. lwd = 2, #line width, 線の太さ main="Directivity" #主題です. )
等,いろいろあります.
plotコマンドに使い慣れていれば,それをサポートしているらしいので,似たような感覚で使えると思います.
R言語は統計解析の計算のみならず,可視化にも優れています.
こういう風に,アンテナの指向性パターンも綺麗に描画できました.めでたし.
# References
http://ww36.tiki.ne.jp/~natsu/
レーダーチャートを描く - langstat blog
グラフィックス | Rで各種グラフの描き方,折れ線,散布図,ヒストグラム,棒グラフ,円グラフ,ボックスプロット
https://cran.r-project.org/web/packages/plotrix/plotrix.pdf
ディスコーンアンテナを作ってみた話
題のとおりです.
まずはモノを見ていただきましょう.
100均の用品をふんだんに利用して安価に仕上げています.
エレメントには100均の園芸用洋ラン線を,形状の保持には100均のまな板に穴をあけて使っています.
同軸の芯線側(ディスク側)はこんな感じで接続,固定しています.
被覆側(コーン側)はこんな感じ.
すべて圧着端子で仕上げています.あとはボルトで接続.
これであとは防水関連の処理をしてあげれば屋外でも使えます.
ブチルテープ(自己融着テープ)や,ホットボンド等で絶縁をします.
割と形状もよく,思ったより良い出来映えでした.
送信用に使うにはなかなか難しいかもしれませんが,広帯域の受信用アンテナとして活用できると思います.
安価に仕上げているので,もしよろしければ作ってみてはいかがでしょうか.