ほしぞloveログ

天体観測始めました。

タグ:python

何日か前の記事で少しだけ書きましたが、分光器のSHG700を使って、太陽望遠鏡のHαエタロンの性能を表すFWHM(Full Width Half Maximum, 半値全幅)を実測してみました。これは太陽望遠鏡のフィルターがFabry-Perot cavityを利用したものだと知った2017年頃からやってみたかったことで、一時は中古の研究用の分光測定器を買うことを本気で考えていました。長年の夢の一つが叶ったことになります。


測定方法

今回SHG700で測定したものは、
  1. 太陽光の散乱光
  2. エタロンの透過光
の2種類です。ここからFWHMまで持っていきます。

1の散乱光は、SHG700を鏡筒から外して単体にして。部屋の中の(直射日光ではない)白い壁に向けます。白く明るい壁ですが、所詮背景光なので光量は大したことはなくて、露光時間を12.8秒でG3M678Mのゲインを400とし、さらにライブスタックで10枚重ねて、十分フラウンホーファー線の構造が見えるようにしました。

2のエタロンの透過光ですが、最初1と同様に太陽の散乱光を使って測定しようと思いました。でも光量が十分でなく、エタロンが共振しない暗いところは十分に見ることができません。太陽光を直接入れて測定するのがいいのですが、あいにくこの日は曇りです。というか、晴れないので痺れを切らしてこの測定を開始したので、太陽が出てないです。代わりに下の写真のようにLEDの小さなライトを使いました。
G1ghTC8aoAAOVvC

PCの画面にも出ていますが、うまくエタロンのComb (櫛形) 構造が見えるようになりました。

ただし測定は結構難しくて、ライトの光の絞り具合とか、ライトと分光器の間の距離だとか、ライトの位置や角度など、うまく合わせないとなかなか綺麗な線が出ません。とりあえず今回は机の上に適当に置いてやりましたが、できるなら光学定盤などを使い安定した測定にしたいです。特に、エタロンは入射光に角度依存性があって、いつかそれも含めて測定したいので、光の角度をきちんと調整できる機構が欲しくなります。


撮影画像

測定した画像は以下のようになります。

まず1の背景光です。
Capture_00001_WithDisplayStretch
太陽光のスペクトルが綺麗に出ているので、この画像から波長のキャリブレーションをすることができそうです。でもこれだけだと、Hα線は目立つのでまだしも、どの線がどの波長なのかよくわかりません。JSol'Exの「Spectrum browser」で見る参照スペクトルと比べてみても、なんか違うように見えます。

下の画像を見るとわかると思いますが、左が今回撮った散乱光、真ん中の細長いSharpCapの画像が以前撮った太陽を直接見たもの、右がJSol'Exの参照画面です。
Fraun_comp_cut

左と真ん中は同じような構造になっているので、まずは背景光がきちんと取れていると判断します。でも右の参照画面の線はかなり実測と違うことがわかります。なので下の画像のように、Hα線と目立って一致しているもう一本の線の波長を調べて、それを基準として他は波長が線形に変化していると仮定して、縦方向の波長を1次の直線でフィットすることにしました。

wavelength_select_cut


2のエタロンですが、本当は透過「率」を知りたいのですが、これは結構難しいとわかりました。まず、エタロンがある場合とない場合の画像を2枚撮影します。まずはエタロンありの画像をLEDライトの位置や角度を変えうまく撮れる状況を作ります。
Capture_00001

その撮影したままの状態をキープしながら、エタロンだけを動かして取り除きます。こうすることで同じ状況で基準光を撮影することができます。
Capture_00001

基準光は一見一定に見えますが、画像の上から下までで緩やかに暗くなっていくことがわかりました。エタロンの透過光のピーク位置もやはり同様に緩やかに暗くなっていくので、基準光で割ることにより、エタロンの透過光のピーク位置が平らに近くなります。

それでも、エタロンがないときにはSHG700の入射口径全体から光が入り、エタロンがあるときにはエタロン前後のレンズ径などに制限された光しか入らないので、透過率が低く出過ぎてしまいます。そのため透過「率」とすることは諦めて、ピーク位置を1とするようように規格化しました。

波長は画像の縦方向で変化しますが、スリットに長さがあるために画像の横方向にもフラウンホーファー線は広がっていて、しかも線が直線にはならずに曲線になっています。(どういった仕組みで曲線になるのか、どう調整したら直線にできるのかの方法は私はまだわかっていないので、こちらもいずれ解決したいです。おそらくスリット位置と回折格子の相対位置で決まるのではと推測しています。)しかも、エタロンの透過光の明るいところと暗いところの幅は横位置によって多少変わります。

今回は1の画像も2の画像も、真ん中あたりの斜めになっていない場所の10ライン程度の縦線を抜き出して、横方向に平均値をとりました。エタロンについては真ん中ら辺が明るい線が一番細いようなので、こちらも真ん中ら辺を選ぶのが一番良さそうです。

背景光のフラウンホーファー線を見ている限り、SHG700の回折格子を触りさえしなければ、撮影ごとの波長のズレのようなものはなさそうなこともわかりました。


波長のキャリブレーション

1の背景光画像のフラウンホーファー線では、波長がリニアに変化すると仮定して、上で決めた基準の2点Hαの6562.81Åと6643.63Åを元に1次の直線でフィットします。この時のあるところの数値と次の数値との差が、1ピクセルあたりに変化する波長となり、今回は0.089Å/pixelとなりました。しかしながら、SIMSPEC SHGで求めた0.091Å/pixelや普段撮影動画をJSol'Exで再構築した際にはこれまで0.091Å/pixelと出ていて、2%ほど結果が異なることがわかりました。

この違いの原因は2点だけを基準として波長が1次的に変化すると仮定したことかと思いますが、今のところはっきりとした原因は不明です。まあ今回は基準点のHα周りのFWHMを求めるのが最大の目的で、そこまで影響はないはずなので、とりあえずこのズレは無視することにします。


エタロンフィッティング

エタロンの透過光ですが、透過光を数値化したものを、基準光で割ったものをグラフにします。
etalon_ok

ここからHα周りのピークを抜き出して、フィッティングします。ピークの高さは右に行くに従って上がっていくようですが、基準光でのノーマライズがうまくいっていないのか、それともこうなるのが正しいのかよくわかりませんでした。Hα周りだけに絞ってしまえば、局所的にはほぼ同じ高さとしてしまっていいでしょう。

フィッティングはFitykというソフトでローレンツ関数やVoigtを使う例がいくつか示されているので、私も同様に試してみましたが、いくつか問題がありそうです。

下の画像は実際にFitykでフィットしてみたものです。
Voigt_cut

一つ目の問題は、これらの関数は基本的にピークの両側は0になることを想定していることです。ところが、エタロンの応答を表す関数は繰り返し構造になるため、ピークとピークの間の透過率が0になりません。ピークとピークの真ん中のちょうど反共振の位置では、エタロンの透過率は、同じ特性の鏡を2枚使うと仮定して、鏡の強度反射率Rと強度透過率Tを使って

(T/(1+R))^2

のような形に書けます。例えばここで、強度透過率T=0.3、強度反射率R=0.7とすると、ピークの真ん中でも(0.3/1.7)^2=0.0311と、3%ほど光を通してしまいます。

Fitykでは、別途定数を用いてフィッティングさせるような手法が取られているようですが、これだと個別の赤い2本の線のうち曲線の方を見てもらえばわかりますが、明らかに実測のピークより細い線でフィッティングされてしまっています。これは結果として、FWHMが小さく出過ぎてしまい、実際よりも性能がいいという間違った結果を出してしまいます。

今回の上の結果では、グラフ右にあるFWHMの数値を見ると、0.65ÅとPSTとしてはにわかに信じられないくらいのいい値が出てしまっています。例えばこのページでも同様の間違いをしていて、HeliostarのエタロンのFWHMがなんと0.3Åと、これも良すぎる値を出してしまっています。ピークの高さの半分のところの幅を見るだけでも、少なくとも0.4Åはあることがパッと見るだけでわかるので、明らかな間違いです。このグラフが出た時に何でこんな良すぎる値になるのかおかしいと思ったのですが、実際に自分でFitykを使ってみることでなぜこんな間違いに陥ったのかがよくわかりました。

二つ目の問題点は、ローレンツ関数やVoigt関数だと、一つのピークのみしかフィットすることができないことです。原理的に、エタロンの透過光の応答のような周期的なものを表すことはできません。このため、周期構造から求めることができる、FSR(Free Spectral Range)をきちんとフィッティングして求めることができません。

FSRはFinesse、FWHMとともにとても重要なパラメーターで、

Finesse = FSR / FWHM

というとてもシンプルな関係があります。Finesseはπで割って2をかけると、エタロン内での光の折り返し回数をすぐに計算できる、非常に重要な指標となります。FSRはエタロンの2枚の鏡の間の距離と反比例関係にあるので、FSRがわかるとエタロン間の距離を直接求めることができます。このように、複数のピークを含めてフィッティングしてFSRを求めることはかなり意義があると言えます。

では、なぜこれまであまり周期的な関数でフィットされてこなかったのでしょうか?これは推測なんですが、単に関数が結構複雑になるためにあまり挑戦してこなかっただけなのかと思います。少なくともFItykのような既存のソフトでフィットするのはかなり大変になりそうです。

今回は周期的な関数を書き下して、自分でpythonでコードを書いて、いくつかのピークをまとめてフィッティングしてみました。結果は以下のようになります。
fit_result_ok

フィッティング曲線がきちんと周期的に出ること、ピークとピークの間が0にならないことがわかるかと思います。ただし、ピークとピークの間の暗い部分が実測とフィッティング曲線でずれてしまっています。これは鏡のロスを考えないで、R+T=1という理想的な鏡を考えてしまったことに由来します。ロスを考えるとさらに複雑になるので、今回は諦めました。それでもFWHMの推定は、ピークの高さをきちんと0を基準に考えているので現実により近い値になっているはずです。


パラメータなど

実際の計算手順としては、フィッティングパラメータとして使った鏡の反射率と透過率、鏡間の距離がまずわかります。鏡の反射率からFinesseが計算でき、鏡間の距離からFSRがけいさんできます。FinesseとFSRがわかると、FWHMがわかるというわけです。下に少しだけ式を書いておきました。

代表的なパラメータはグラフの中に書き込んでおきましたが、今回分かったエタロンの特性を表すパラメーターは以下の通りです。
  • 鏡の振幅反射率、振幅透過率: r, t
  • 鏡の強度反射率、強度透過率: R = r^2 = 0.70, T = t^2 = 0.30
  • キャビティーの鋭さを表すFinesse = π r/(1-R) = 8.75
  • エタロンを構成する鏡と鏡の間の距離 = 0.313 [mm]
  • 周期の幅を表すFSR (Free Spectrul Range) = 6.88 [Å]
  • エタロンの性能を表すFWHM = FSR/Finesse=0.787 [Å]
  • 光の折り返し回数: Finesse *2/π = 5.6 [回 (片道)]
目的のFWHMは0.787 [Å] と出ましたが、公称1 [Å] 以下という値と比べてもそこそこ信頼性のある数字になったのかと思います。FWHMだけでなく、他の重要なパラメータもかなりの精度でわかり、PSTエタロンの特性がかなり特定できたと言っていいかと思います。長年の疑問にやっと答えが出たと言えそうで、かなり嬉しいです。
逆に、今回の測定でまだわからないことは
  • 光の入射角度の依存性
  • Hαからのピークの中心波長のずれ (入射光の角度と、エタロン回転調整をいじっていないため)
  • 個々の鏡の反射率と透過率がどれくらい違うか (2枚の鏡の反射率と透過率を同じと仮定したため)
  • 鏡のロス
などになります。今後の課題としたいと思います。


まとめと今後

手に入れたSHG700で、手持ちのPSTエタロンの透過特性を、うまく測定することができました。角度依存性などの課題はまだ残されていますが、目的のFWHMが測定でき、これまでわからなかった鏡の反射率、ミラー間の距離やFSR、フィネスまで確定できたのはかなり満足感が高いです。

今後やりたいことが、エタロン以外にも太陽望遠鏡でに必須の、BFの測定とかERFの測定です。他にも、ナローバンドフィルターやワンショットナローバンドフィルターも、メーカーが謳っている半値幅が本当に出ているのか、実測してみたいと思っています。




一連の皆既月食記事の8本目。いい加減にそろそろ終わりにするつもりです。前回の記事はこちらから。



Hough変換

今回はHough変換を使ったサークル抽出で、月食中の月の位置合わせに挑戦です。FS-60CBで撮影した4時間分の広角の画像から月が画面中心に来るように位置変換をします。その後、タイムラプス映像にしてみます。

Hough変換は画像の中から特徴的な形を抽出するアルゴリズムで、その中に円を抽出する関数があります。Hough変換は色々な環境で使えますが、今回はOpenCVで用意されている関数をPythonで使うことにしました。検索するとサンプルプログラムなどたくさん出てくるのと、今回は対処療法で組んでいったのの積み重ねなので、あまりに汚いコードで人様に見せれるようなものではないです。なので、どうパラメータを取ったかだけの説明にとどめることにします。

基本的にはあるフォルダにある、月が映っているたくさんのファイルを順次読み込み、カラー画像をグレースケールに変換して、HoughCircles関数を呼び出すだけです。Hough関数は以下のようにしました。

HoughCircles(gray, cv2.HOUGH_GRADIENT, dp=1, minDist=1500, param1=30, param2=2, minRadius=350, maxRadius=400)

パラメータがいくつかありますが、少しだけコツを書いておきます。
  • HOUGH_GRADIENTとHOUGH_GRADIENT_ALTは両方試しましたが、ALT付きの方が誤検出が多かったのでALT無しの方にしました。
  • dpは1以下も試しましたが、位置精度に違いはあまりなかったです。大きな値にすると精度が悪くなりました。
  • 円を複数検出するわけではないので、minDistは1000とかの相当大きな値にしておきます。
  • param1は最終的に50程度にしました。大きすぎたり小さすぎたり値では検出できなかったりしますが、位置精度にはあまり影響ないようです。
  • param2は位置精度に影響があるようです。大きすぎると精度が悪くなりますが、5以下くらいだと精度はこれ以上変わらないようです。
  • minRadiusとmaxRadiusは、月の大きさは一定なのでその半径を挟むような値を取ると効率が良いでしょう。

とりあえず結果を示します。


結局出た精度はこれくらいです。もうちょっとビシーっと止まってくれればいいのですが、まだちょこちょこズレています。

実は精度を出そうとして色々試してみてます。Hough変換のパラメータをいじるのはもちろんなのですが、大きなものは
  • グレー化後に、あらためて画像の2値化
  • グレー化後に、あらためて輝度やコントラストをいじる
  • Adobe Premiereのブレ補正
  • 一旦Hough変換でラフに位置合わせして、少し大きめに切り抜いて、再度Hough変換
くらいでしょうか。でも結局どの方法も精度を劇的に上げることはできませんでした。Hough変換で目で見てもおかしいズレがあるので、そこだけでもずれが直ったらと思ったのですが、やはりダメなものはダメで、どうも苦手な画像があるようです。

今回はFS-60CBで広角で撮ったものから抜き出しているので、解像度が良くないこと。さらに動画の中から最初の1枚だけを抜き出しているので、大気揺らぎなどで多少のブレがあります。もしかしたらスタックをすると平均化されるので、細かいブレは少なくなるのかもしれません。スタックしてから、Hough変換ではなく、特徴点を抜き出して合わせるとかいう処理をした方が精度が出るかも知れませんが、ちょっと力尽きたので、今回はここまでにします。

ついでに同じ位置合わせの手法で、いくつか画像を抜き出して天王星の潜入画像をつくってみました。位置合わせ後、比較明合成しただけです。これくらいの精度ではそこそこ合って見えるんですけど、やはり長時間のタイムラプスだときついかもです。

uraus_in

まとめ

今回の皆既月食に関する記事はこれで一旦終わりにします。

撮影後、1ヶ月以上にわたって楽しむことができてかなり満足しました。その一方、そろそろ月も飽きてきました。ただ、まだ未処理ファイルが大量にあるので、気が向いたらもう少し追加で書くかもしれません。

先の記事でも書きましたが、今回は月食に関してはかなり満足して撮影できました。大きな課題は地球本影を位置補正することなく撮影することですが、次回月食ではこれ一本に絞ることにするかもしれません。

ここ数日、BlurXTerminatorがすごいことになっています。次は少しこちらの方に時間を費やそうかと思います。












まだまだ自分の中では皆既月食マイブーム状態です。今回は地球の影が止まるような位置をどう計算するかの一連の記事の、少し脱線するような記事になります。


「ほんのり光芒」さんとのコラボ? 

前回の記事でコメントをいただいた「ほんのり光芒」のみゃおさんがブログ記事の方で、地球本影の固定に関してかなり細かい検討をされています。


上記ページからリンクを辿れますが、かなり以前から考えられていたようで、私のようなにわか月食撮影とは歴史が全然違います。天リフさんのピックアップにも取り上げられていて、ほしぞloveログとのちょっとしたコラボのような様相を呈しています。




視野ズレの計算プログラム

前回までで、地球の自転による観測者の位置変化で視野に大きなずれができて、地球の影の形がひしゃげることを、多少定量的に見積もってみました。


大まかな見積もりと、実際にどれくらいずれるかは、少なくともオーダーレベルでは一致しているようなので、今回はもう少し精度良く計算できないかを考えてみます。といっても、自分だけで考えるのはそろそろ限界で、Webで少し検索してみました。すると老猫こてつさんという方が各種軌道計算などされて、その中で求めていた月の視野位置のずれそのものをpythonで計算してくれていることがわかりました。


プログラムのソースを公開してくれているので、そのまま計算することができ、とても助かります。ただし、それら計算式をどう求めているのかの記述はないので、そのプログラムの出典を調べるためにブログの過去記事を読んでいくと、どうやら中野主一さんの昔のBASICのコードをpythonに置き換えてくれているようです。


参考書籍

中野主一さんといえば、最近では小説「星になりたかった君と」で「長野秀一」という名前で出てきた重要な役割を担う老人のモデルとなった方。おそらく昔の天文少年にとっては憧れのアマチュア天文家で、多くの天文雑誌で連載をもち、アマチュアながら元国立天文台長の古在由秀先生から計算を依頼されるなど、軌道計算の大家です。私はその当時マイコン少年でしたが、大人になって天文を始めてから、昔使っていたマイコンで実用的な天体計算をしていたことを知って、とても感動したことを覚えています。実は中野主一さん、星を初めて少しした2018年に一度お会いして、お話しさせていただいたことがあります。私はかなり緊張していたのですが、気さくにお話ししていただきました。講演も聞かせてもらったのですが、当時どう計算を進めていたか、プロから依頼された当時の様子などのお話で、ユーモアたっぷりの講演で今でも心に残っています。

そんなわけで早速、サイトに載っていた参考文献を注文しました。長沢先生の「天体の位置計算」は以前から持っていたもので、まだ普通に販売されているようですが、他の3冊は流石に古本でかろうじて見つかるくらいでした。
72AF0666-ACA0-4136-8D0B-6D5DB7F44CF5

その中で一番古い「マイコン宇宙講座」は昭和55年の初出のもので、月の計算そのものの章があり、式の解説もあったので理解しやすかったです。

「マイコンが解く天体の謎」は昭和57年出版で、使われている言語はなんとF-BASICですよ。内容はプラネタリウムのようなものを実現することが中心ですが、実は私、中1の時にFM-NEW7を中古で買って、しゃぶり尽くした口です。その数年前に出た本のようなので、恐らくFM-7が出た時で、実際にはFM-8で組まれた時代のプログラムですね。

一番新しい「天体の軌道計算」は1992年なので、前の2さつからはかなり経っていて、プログラムも複雑になり、さらに精度を求めているような内容になっています。

とりあえずは、月の視野ズレの計算方法が載っている「マイコン宇宙講座」をもとに、老猫こてつさんのpythonコードを使わせていただいて計算を進めようと思います。

ただしこれはまだ、「月」の視野ズレを追いかけるプログラムで、一番求めたい月食中の「地球の本影」を追うものではありません。でもこれらの計算の延長上に、それもそう遠くないところに地球本影を求めることができるのではと期待しています。


赤道儀の制御

あともう一つ、仮に地球本影の視野ズレも含んだ位置を計算できたとして、それをどう赤道儀に伝えるかですが、彗星を追うメカトーフ法というのが応用できるかもしれません。ただし、地球の本影を追うというようなものは見つからなかったので、実際にどういう方法で赤道儀に伝えるのか、どういうデータ形式なのか、地球本影に応用できるかなど、もう少し調べる必要があります。

ガイドソフトのPHD2にもメカトーフに相当するような機能があるらしいのですが、どうも1次の傾きでガイド信号に補正信号を加えていくようなものらしいです。視野ズレのような複雑な動きはできないかもしれませんが、1次補正だけでも近似的にそこそこ地球の影の形はうまく出るのではと思います。もしくは撮影途中で何度か係数を変えるかとかでしょうか。


今後

視野ズレの計算方法や赤道儀の制御など、まだ直接ではないですが答えにつながりそうな幾つかの見通しは出てきました。次回の日本での地球本影が見える月食は2023年10月23日の部分月食だそうです。1年近くあるので、じっくり準備したいと思います。


 
 
 
 
 
 
 
 

このページのトップヘ