Python - 赤道・黄道座標の変換！

[ プログラミング, 暦・カレンダー ] [ Python, カレンダー ]

以前、赤道直交座標と黄道直交座標を相互に変換したり、直交座標と極座標を相互に変換したりする RubyGems ライブラリを作成しました。

今回は、同様のことを Python で行ってみました。（但し、PyPI ライブラリではない）

0. 前提条件

  • Python 3.6.5 での作業を想定。
  • 天文学的な計算については疎いため、誤りがあるかもしれない。

1. アルゴリズムについて

過去記事を参照。

2. Python スクリプトの作成

  • ライブラリ内で NumPy の matrix を使用。

次のスクリプトは実行部分。

conv_coord.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

#! /usr/local/bin/python3.6
"""
赤道・黄道座標変換

  date          name            version
  2018.04.27    mk-mode         1.00 新規作成

Copyright(C) 2018 mk-mode.com All Rights Reserved.
"""
import math
import sys
import traceback
from lib import coord as lcd


class ConvCoord:
    # 黄道傾斜角(単位: rad)
    EPS = 23.43929 * math.pi / 180
    # 元の赤道直交座標
    # (ある日の地球重心から見た太陽重心の位置(単位: AU))
    POS = [ 0.99443659220700997281, -0.03816291768957833647, -0.01655177670960058384]

    def exec(self):
        try:
            self.rect_ec = lcd.rect_eq2ec(self.POS, self.EPS)
            self.rect_eq = lcd.rect_ec2eq(self.rect_ec, self.EPS)
            self.pol_eq, self.r = lcd.rect2pol(self.rect_eq)
            self.pol_ec = lcd.pol_eq2ec(self.pol_eq, self.EPS)
            self.pol_eq = lcd.pol_ec2eq(self.pol_ec, self.EPS)
            self.rect_eq_2 = lcd.pol2rect(self.pol_eq, self.r)
            self.__display()
        except Exception as e:
            raise

    def __display(self):
        """ 結果出力 """
        try:
            pass
            print((
                "元の赤道直交座標:\n  {}\n"
                "黄道直交座標に変換:\n  {}\n"
                "赤道直交座標に戻す:\n  {}\n"
                "赤道極座標に変換:\n  {} (R = {})\n"
                "黄道極座標に変換:\n  {}\n"
                "赤道極座標に戻す:\n  {}\n"
                "赤道直交座標に戻す:\n  {}"
            ).format(
                self.POS,
                self.rect_ec,
                self.rect_eq,
                self.pol_eq, self.r,
                self.pol_ec,
                self.pol_eq,
                self.rect_eq_2
            ))
        except Exception as e:
            raise


if __name__ == '__main__':
    try:
        obj = ConvCoord()
        obj.exec()
    except Exception as e:
        traceback.print_exc()
        sys.exit(1)

次の3つのスクリプトはライブラリ部分。（lib ディレクトリを作成後、その配下へ）

lib/coord.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101

"""
Modules for coordinates
"""
import math
import numpy as np
import sys
from lib import matrix        as lmtx
from lib import trigonometric as ltrg


def rect_eq2ec(rect, eps):
    """ 直交座標：赤道座標 -> 黄道座標.

    :param  list rect: 赤道直交座標
    :param  float eps: 黄道傾斜角 (Unit: rad)
    :return list     : 黄道直交座標
    """
    try:
        mtx_rect = np.matrix([rect]).transpose()
        r_mtx = lmtx.r_x(eps)
        return lmtx.rotate(r_mtx, mtx_rect).A1.tolist()
    except Exception as e:
        raise

def rect_ec2eq(rect, eps):
    """ 直交座標：黄道座標 -> 赤道座標.

    :param  list rect: 黄道直交座標
    :param  float eps: 黄道傾斜角 (Unit: rad)
    :return list     : 赤道直交座標
    """
    try:
        mtx_rect = np.matrix([rect]).transpose()
        r_mtx = lmtx.r_x(-eps)
        return lmtx.rotate(r_mtx, mtx_rect).A1.tolist()
    except Exception as e:
        raise

def pol_eq2ec(pol, eps):
    """ 極座標：赤道座標 -> 黄道座標.

    :param  list  pol: 赤道極座標
    :param  float eps: 黄道傾斜角 (Unit: rad)
    :return list     : 黄道極座標 [λ, β]
    """
    try:
        alp, dlt = pol
        lmd = ltrg.comp_lambda(alp, dlt, eps)
        bet = ltrg.comp_beta(alp, dlt, eps)
        return [lmd, bet]
    except Exception as e:
        raise

def pol_ec2eq(pol, eps):
    """ 極座標：赤道座標 -> 黄道座標.

    :param  list  pol: 赤道極座標
    :param  float eps: 黄道傾斜角 (Unit: rad)
    :return list     : 黄道極座標 [α, δ]
    """
    try:
        lmd, bet = pol
        alp = ltrg.comp_alpha(lmd, bet, eps)
        dlt = ltrg.comp_delta(lmd, bet, eps)
        return [alp, dlt]
    except Exception as e:
        raise

def rect2pol(rect):
    """ 直交座標 -> 極座標

    :param  list rect: 直交座標
    :return list     : 極座標 [[λ, φ], 距離]
    """
    try:
        x, y, z = rect
        r = math.sqrt(x * x + y * y)
        lmd = math.atan2(y, x)
        phi = math.atan2(z, r)
        if lmd < 0:
            lmd %= math.pi * 2
        d = math.sqrt(x * x + y * y + z * z)
        return [[lmd, phi], d]
    except Exception as e:
        raise

def pol2rect(pol, r):
    """ 極座標 -> 直交座標

    :param  list pol: 極座標
    :return list    : 直交座標
    """
    try:
        lmd, phi = pol
        r_mtx = lmtx.r_y(phi)
        r_mtx = lmtx.r_z(-lmd, r_mtx)
        return lmtx.rotate(
            r_mtx, np.matrix([[r], [0.0], [0.0]])
        ).A1.tolist()
    except Exception as e:
        raise
lib/matrix.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83

"""
Modules for matrixes
"""
import numpy as np


R_UNIT = np.eye(3, dtype="float64")


def r_x(phi, r_src=R_UNIT):
    """ 回転行列生成(x軸中心)
          ( 1      0          0     )
          ( 0  +cos(phi)  +sin(phi) )
          ( 0  -sin(phi)  +cos(phi) )

    :param  np.matrix r_src: Rotation matrix
    :param  float       phi: Angle (Unit: rad)
    :return np.matrix r_dst: Rotated matrix
    """
    try:
        s, c = np.sin(phi), np.cos(phi)
        r_mx = np.matrix([
            [1,  0, 0],
            [0,  c, s],
            [0, -s, c]
        ], dtype="float64")
        return r_mx * r_src
    except Exception as e:
        raise

def r_y(theta, r_src=R_UNIT):
    """ 回転行列生成(y軸中心)
          ( +cos(theta)  0  -sin(theta) )
          (     0        1      0       )
          ( +sin(theta)  0  +cos(theta) )

    :param  np.matrix r_src: Rotation matrix
    :param  float     theta: Angle (Unit: rad)
    :return np.matrix r_dst: Rotated matrix
    """
    try:
        s, c = np.sin(theta), np.cos(theta)
        r_mx = np.matrix([
            [c, 0, -s],
            [0, 1,  0],
            [s, 0,  c]
        ], dtype="float64")
        return r_mx * r_src
    except Exception as e:
        raise

def r_z(psi, r_src=R_UNIT):
    """ 回転行列生成(z軸中心)
          ( +cos(psi)  +sin(psi)  0 )
          ( -sin(psi)  +cos(psi)  0 )
          (     0          0      1 )

    :param  np.matrix r_src: Rotation matrix
    :param  float       psi: Angle (Unit: rad)
    :return np.matrix r_dst: Rotated matrix
    """
    try:
        s, c = np.sin(psi), np.cos(psi)
        r_mx = np.matrix([
            [ c, s, 0],
            [-s, c, 0],
            [ 0, 0, 1]
        ], dtype="float64")
        return r_mx * r_src
    except Exception as e:
        raise

def rotate(r, pos):
    """ 座標回転

    :param  np.matrix r    : 回転行列
    :param  np.matrix pos  : 回転前直交座標
    :return np.matrix pos_r: 回転後直交座標
    """
    try:
        return r * pos
    except Exception as e:
        raise
lib/trigonometric.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

"""
Modules for trigonometrics
"""
import math
import numpy as np


def comp_lambda(alp, dlt, eps):
    """ λ の計算
        * λ = arctan((sinδ sinε + cosδ sinα cosε ) / cosδ cosα)

    :param  float alp: α (Unit: rad)
    :param  float dlt: δ (Unit: rad)
    :param  float eps: ε (Unit: rad)
    :return float lmd: λ (Unit: rad)
    """
    try:
        a = math.sin(dlt) * math.sin(eps) \
          + math.cos(dlt) * math.sin(alp) * math.cos(eps)
        b = math.cos(dlt) * math.cos(alp)
        lmd = math.atan2(a, b)
        if lmd < 0:
            lmd %= math.pi * 2
        return lmd
    except Exception as e:
        raise

def comp_beta(alp, dlt, eps):
    """ β の計算
        * β = arcsisn((sinδ cosε - cosδ sinα sinε)

    :param  float alp: α (unit: rad)
    :param  float dlt: δ (unit: rad)
    :param  float eps: ε (unit: rad)
    :return float bet: β (unit: rad)
    """
    try:
        a = math.sin(dlt) * math.cos(eps) \
          - math.cos(dlt) * math.sin(alp) * math.sin(eps)
        return math.asin(a)
    except Exception as e:
        raise

def comp_alpha(lmd, bet, eps):
    """ α の計算
        * α = arctan((-sinβ sinε + cosβ sinλ cosε ) / cosβ cosλ)

    :param  float lmd: λ (unit: rad)
    :param  float bet: β (unit: rad)
    :param  float eps: ε (unit: rad)
    :return float alp: α (unit: rad)
    """
    try:
        a = -math.sin(bet) * math.sin(eps) \
          +  math.cos(bet) * math.sin(lmd) * math.cos(eps)
        b =  math.cos(bet) * math.cos(lmd)
        alp = math.atan2(a, b)
        if a < 0:
            alp %= math.pi * 2
        return alp
    except Exception as e:
        raise

def comp_delta(lmd, bet, eps):
    """ δ の計算
        * δ = arcsisn((sinβ cosε + cosβ sinλ sinε)

    :param  float lmd: λ (unit: rad)
    :param  float bet: β (unit: rad)
    :param  float eps: ε (unit: rad)
    :return float dlt: δ (unit: rad)
    """
    try:
        a = math.sin(bet) * math.cos(eps) \
          + math.cos(bet) * math.sin(lmd) * math.sin(eps)
        return math.asin(a)
    except Exception as e:
        raise

3. Python スクリプトの実行

  • 黄道傾斜角εや元の赤道直交座標(3x1行列)は実行スクリプト内で指定する。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

$ ./conv_coord.py
元の赤道直交座標:
  [0.99443659220701, -0.038162917689578336, -0.016551776709600584]
黄道直交座標に変換:
  [0.99443659220701, -0.04159771108146677, -5.622172494390565e-06]
赤道直交座標に戻す:
  [0.99443659220701, -0.038162917689578336, -0.016551776709600584]
赤道極座標に変換:
  [6.24482770879939, -0.01663059980037209] (R = 0.9953062370542631)
黄道極座標に変換:
  [6.241379248856413, -5.648686087788284e-06]
赤道極座標に戻す:
  [6.24482770879939, -0.01663059980037209]
赤道直交座標に戻す:
  [0.99443659220701, -0.03816291768957855, -0.01655177670960066]

これらの座標変換は、天体の各種計算でよく使用します。

以上。