次のような状況を考えることにします。観測場所は、ニューヨークにあるBear Mountain State Park。そこに自前の望遠鏡を持ち込んで、M33を観望するとします。前後にスケジュールが詰まっているため、観望できるのは午後１１時前後のみです。この時間にM33は観望できる高さまで上がっているでしょうか。このような問題は、例えば以下の様な方法で調べることができます。
まず、作業に必要なパッケージを読み込みます。
>>> import numpy as np >>> from astropy import units as u >>> from astropy.time import Time >>> from astropy.coordinates import SkyCoord, EarthLocation, AltAz
次に、M33の座標値をネットから取得し、観測場所の座標値を指定します。天体の座標値は、（SIMBADに登録されているものであれば、SkyCoord.from_nameで以下のように所得することができます。）
>>> m33 = SkyCoord.from_name('M33')
>>> bear_mountain = EarthLocation(lat=41.3*u.deg, lon=-74*u.deg, height=390*u.m)
次に、観測場所の時間を指定します。まず１行目で、世界標準時からのオフセットをしていし、２行目で、世界標準時からオフセットを差し引いて観測場所の時間を表すオブジェクトを作ります。
>>> utcoffset = -4*u.hour # ニューヨークなのでアメリカ東部標準時（UTC-4h）を設定します。
>>> time = Time('2012-7-12 23:00:00') - utcoffset
設定した時間と場所の情報を使って、ICRSからAltAzへの座標変換を行い、フォーマットを整えて出力します。
>>> m33altaz = m33.transform_to(AltAz(obstime=time,location=bear_mountain))
>>> "M33's Altitude = {0.alt:.2}".format(m33altaz)
"M33's Altitude = 0.13 deg"
ここまでの流れをまとめておきます。
>>> import numpy as np
>>> from astropy import units as u
>>> from astropy.time import Time
>>> from astropy.coordinates import SkyCoord, EarthLocation, AltAz
>>> m33 = SkyCoord.from_name('M33')
>>> bear_mountain = EarthLocation(lat=41.3*u.deg, lon=-74*u.deg, height=390*u.m)
>>> utcoffset = -4*u.hour # Eastern Daylight Time
>>> time = Time('2012-7-12 23:00:00') - utcoffset
>>> m33altaz = m33.transform_to(AltAz(obstime=time,location=bear_mountain))
>>> "M33's Altitude = {0.alt:.2}".format(m33altaz)
"M33's Altitude = 0.13 deg"
とりあえず、観測したい場所と時間における高度角を求めることができました。しかし、この高度は観測するには低すぎるようです。もっといろいろな時間で計算して、観測に適切な時間の範囲を調べたいところです。そのような確認をするには、１つの場合を計算するだけでなく、時間に対して高度がどのように変化するか図にした方が分かりやすいでしょう。また、高度角だけでなくエアマスも合わせて計算するとより便利かもしれません。ではそのような作業を、実際にやってみましょう。まず、計算原点となる時刻のtimeオブジェクトを作ります。先ほど作った国際標準時への補正 utoffsetを利用します。
>>> midnight = Time('2012-7-13 00:00:00') - utcoffset
次に、－２時間から＋７時間までを、１００分割した時間のnumpy配列を作ります。この配列は、astropy.unitの単位を後ろにつけているので、単なる数値の配列ではなく、時間の単位を持っているところがポイントです。
>>> delta_midnight = np.linspace(-2, 7, 100)*u.hour
transform_toは配列を処理することができるので、今しがた作ったnumpy配列を使って、大量の変換を一括して行います。
>>> m33altazs = m33.transform_to(AltAz(obstime=midnight+delta_midnight, location=bear_mountain))
出来上がった (Alt, Az) の数値ペアを１００個含んだ配列を使って、まず、時間に対するエアマスの変化をプロットしてみましょう。matplotlibを利用して以下の様な方法で図を作ることができます。
>>> import matplotlib.pyplot as plt
>>> plt.plot(delta_midnight, m33altazs.secz)
>>> plt.xlim(-2, 7)
>>> plt.ylim(1, 4)
>>> plt.xlabel('Hours from EDT Midnight')
>>> plt.ylabel('Airmass [Sec(z)]')
 |
>>> from astropy.coordinates import get_sun >>> delta_midnight = np.linspace(-12, 12, 1000)*u.hour >>> times = midnight + delta_midnight >>> altazframe = AltAz(obstime=times, location=bear_mountain) >>> sunaltazs = get_sun(times).transform_to(altazframe) >>> m33altazs = m33.transform_to(altazframe)
以下各業ごとの説明です。
１行目： get_sunを読み込みます。
２行目： 時間のNumpy配列を作ります。-12時間から+12時間まで１時間刻みの配列にします。
３行目： midnightという変数に保存されている時間を起点に、-12時間から+12時間までの時間の配列を作り、timesという変数に保存します。
４行目： 時間と観測場所の情報を含んだAltAzオブジェクトを作り、altazframeという名前で保存します。
５行目： timesに保存された時間における太陽の位置情報をAltAzに変換して、sunaltazsという変数に保存します。
６行目： 観測天体の位置情報を、AltAz座標に変換して、m33altazsという変数に保存します。
上記のように必要な数値を計算した後は、matplotlibで図を描画します。
>>> plt.plot(delta_midnight, sunaltazs.alt, color='y', label='Sun')
>>> plt.scatter(delta_midnight, m33altazs.alt, c=m33altazs.az, label='M33', lw=0, s=8)
>>> plt.fill_between(delta_midnight, 0, 90, sunaltazs.alt < -0*u.deg, color='0.5', zorder=0)
>>> plt.fill_between(delta_midnight, 0, 90, sunaltazs.alt < -18*u.deg, color='k', zorder=0)
>>> plt.colorbar().set_label('Azimuth [deg]')
>>> plt.legend(loc='upper left')
>>> plt.xlim(-12, 12)
>>> plt.xticks(np.arange(13)*2 -12)
>>> plt.ylim(0, 90)
>>> plt.xlabel('Hours from EDT Midnight')
>>> plt.ylabel('Altitude [deg]')
