astropy を使ってドップラートラッキングをする

# 背景
### はじめに
- [天の川銀河](http://astro-dic.jp/milky-way-galaxy/)をただよう[恒星](http://astro-dic.jp/fixed-star/)や[星間ガス](http://astro-dic.jp/interstellar-matter/)は、銀河中心の周りを回転しています。
    - [:movie_camera:  かっこいいシミュレーション動画](http://www.jicfus.jp/jp/research/findings/image_movie/milkyway1/)

- 私たちがこういった構造を観測しようとする場合、例えば、[電波分光装置](http://astro-dic.jp/heterodyne-receiver/)を使って[星間ガス](http://astro-dic.jp/interstellar-matter/)を観測します。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/107797/b86ccb69-b2cf-b6f1-60e7-ee6a8e8d0406.jpeg" width=50%>   <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/107797/7fa68761-7535-89b3-349a-65b86badf2d1.png" width=30%>
(Credit: NASA/JPL-Caltech/ESO/R. Hurt, いらすとや)




### 天体の速度と、観測者の速度

- 観測で得られた[輝線](http://astro-dic.jp/emission-line/)は、周波数がズレており、ズレの分だけガスは速度を持っています ([ドップラー効果](https://ja.wikipedia.org/wiki/%E3%83%89%E3%83%83%E3%83%97%E3%83%A9%E3%83%BC%E5%8A%B9%E6%9E%9C))。

- ただし、[ドップラー効果](https://ja.wikipedia.org/wiki/%E3%83%89%E3%83%83%E3%83%97%E3%83%A9%E3%83%BC%E5%8A%B9%E6%9E%9C)は、天体の速度と、観測者の速度の 2 つが影響します。
    - 天体の速度は 100〜10 万年くらいは変わらないのに対し、
    - 観測者は、地球の自転 (1日周期)、公転 (1年周期) で速度が変わります。

- 天体の速度だけを記録したいので、地球の速度とは関係のない、普遍的な静止系を定義してやると便利です。

- 静止系として、電波天文学で用いられているのは、[局所静止基準 (LSR)](http://astro-dic.jp/local-standard-of-rest/)です。
    - [局所静止基準](http://astro-dic.jp/local-standard-of-rest/)は、太陽系の中心を原点とした静止系です。
    - ただし、[局所静止基準](http://astro-dic.jp/local-standard-of-rest/)の速度は、太陽系の速度を採用しませんでした。
        - なぜなら、局所静止基準には、太陽系の速度ではなく、天の川銀河の中で、太陽の位置での、銀河としての大局的な速度を採用したいからです。
        - 太陽系の速度は、この場所での銀河としての大局的な速度からは外れており、適しません。
    - [局所静止基準](http://astro-dic.jp/local-standard-of-rest/)の速度には、太陽系の位置での銀河としての大局的な速度 (太陽系周辺のいくつかの恒星の速度から求める) が使われます。

### 天体速度を取得する (ドップラートラッキング)

- したがって、「天体の速度を記録する」とは、局所静止基準に対する速度 $V_{LSR}$ を記録する、ということになります。
    - 観測者が地球上で観測した天体信号から、LSR に対する観測者の速度 $V_{observer}$ を引けば $V_{LSR}$ が得られます。
        -  $V_{observer}$ は、
            -  LSRに対する太陽速度 : $V_{sun}$
            -  太陽に対する地球の公転速度 : $V_{earth}$
            -  地球に対する観測者の自転速度 : $V_{location}$
            -  からなります。
    -  $V_{observer}$ を引く処理のことを、ドップラートラッキングと呼びます。

### 太陽速度が変わってしまう...！

- LSR に対する太陽速度 $V_{sun}$ は、観測精度が向上したり、LSR の定義の違いによって、変わってしまいます。
    - 例えば、最近の研究 [Schönrich et al. (2010)](https://ui.adsabs.harvard.edu/abs/2010MNRAS.403.1829S/abstract) によれば、$(U, V, W) = (11.1, 12.24, 7.25)$ ${\rm km\,s}^{-1}$ とされました。
        - $U$ : 太陽速度のうち、銀河中心方向への成分 (銀経 = $0^\circ$ の方向)
        - $V$ : 太陽速度のうち、銀河面に沿って、銀河中心へ垂直な成分 (銀経 = $90^\circ$ の方向)
        - $W$ : 太陽速度のうち、銀河面に垂直な成分

- これでは、採用する結果が異なれば、同じ天体の観測であっても $V_{LSR}$ が変わってしまって不便です。
- そこで、電波天文学における $V_{sun}$ は次のように定義することになっています。
    - 太陽は局所静止系に対して1900年分点で赤経$18{\rm h}$、赤緯$+30^{\circ}$の方向に$+20\,{\rm km\,s}^{-1}$で運動している ([天文学辞典](http://astro-dic.jp/local-standard-of-rest/)より引用)
    - この時、$(U, V, W) = (10.27, 15.32, 7.74)$ ${\rm km\,s}^{-1}$ です。


# astropy で計算してみよう！

それでは、[astropy](https://www.astropy.org/) を使って、実際に $V_{observer}$ を計算し、ドップラートラッキングをして見ましょう。

### 太陽速度 Vsun = (U, V, W) を計算する

```python
import astropy.coordinates
from astropy.units import km, s, deg
from numpy import cos, sin

# 太陽の運動方向を定義
dir_sun = astropy.coordinates.SkyCoord(
    ra = 18 * 15 * deg,                   # R.A. = 18 h
    dec = 30 * deg,                       # Dec. = 30 deg
    frame = 'fk4',                        
    equinox = astropy.time.Time('B1900'), # 1900 年分点を指定する
).galactic

# 太陽の速さを定義
v_sun = 20 * km / s

# 太陽速度を、銀河面に直交した座標系に成分分解し、U, V, W を計算
U = v_sun * cos(dir_sun.b) * cos(dir_sun.l)
V = v_sun * cos(dir_sun.b) * sin(dir_sun.l)
W = v_sun * sin(dir_sun.b)

# CartesianDifferential 型にしておく
v_bary = astropy.coordinates.CartesianDifferential(U, V, W)

# >>> v_bary
# <CartesianDifferential (d_x, d_y, d_z) in km / s
#     (10.27059164, 15.31741091, 7.73898381)>
```

- これで、$V_{sun}$ (変数名は v_bary) が計算できました。
- とても簡単ですね。

### Vobserver を計算する
```python
import time
import astropy.coordinates
import astropy.time
from astropy.units import deg, m
from numpy import cos, sin

# 観測者の地球上の座標を設定します。例は長野県野辺山です。
loc1p85 = astropy.coordinates.EarthLocation(
        lon = 138.472153 * deg,
        lat = 35.940874 * deg,
        height = 1386 * m,
)

# 観測時刻を定義します
tobs = astropy.time.Time(time.time(), format='unix')

# 観測方向を定義します。例は Ori-KL です。
target = astropy.coordinates.SkyCoord(
    ra = 83.809 * deg,
    dec = -5.372639 * deg,
    frame = 'icrs',
    obstime = tobs,     # 時刻を設定する
    location = loc1p85, # 観測者座標も設定する
)

# V_observer (LSR 系からみた観測者の速度) の計算
v_observer = astropy.coordinates.SkyCoord(
    # 観測者位置を、観測時刻での GCRS 座標系に変換する
    loc1p85.get_gcrs(target.obstime)
).transform_to(
    # それを、LSR 座標系に変換する
    astropy.coordinates.LSR(v_bary=v_bary) # 太陽速度を指定する
                                           # デフォルトでは、Schönrich et al. (2010) が使われてしまう
).velocity

# V_obs のうち、天体方向の成分を取り出す
v_correction = v_observer.d_x * cos(target.icrs.dec) * cos(target.icrs.ra) + \
               v_observer.d_y * cos(target.icrs.dec) * sin(target.icrs.ra) + \
               v_observer.d_z * sin(target.icrs.dec)

# tobs, ra, dec には配列も渡せます
```
- あらかじめ計算しておいた $V_{sun}$ (変数名 v_bary) を用いて、$V_{observer}$ を計算し、さらに、$V_{observer}$ のうち、観測天体の方向の成分を v_correction として計算しました。
- これで、観測された速度から v_correction を引けば $V_{LSR}$ を取得できます。
- コードの大半は座標の定義です。計算はとても簡単ですね。astropy 最高。

#### メモ
- WARNING: failed to download http://maia.usno.navy.mil/ser7/finals2000A.all が出る場合
    - [IERS](https://ja.wikipedia.org/wiki/%E5%9B%BD%E9%9A%9B%E5%9C%B0%E7%90%83%E5%9B%9E%E8%BB%A2%E3%83%BB%E5%9F%BA%E6%BA%96%E7%B3%BB%E4%BA%8B%E6%A5%AD) のテーブルの保存元が、[アメリカ海軍天文台](https://ja.wikipedia.org/wiki/%E3%82%A2%E3%83%A1%E3%83%AA%E3%82%AB%E6%B5%B7%E8%BB%8D%E5%A4%A9%E6%96%87%E5%8F%B0)から、NASAに変わった
        - 新しいダウンロード先リンク : ftp://cddis.gsfc.nasa.gov/pub/products/iers/finals2000A.all
    - 新しい astropy では対応されているので、アップデートしましょう