Instantly share code, notes, and snippets.

# komasaru/conv_coord.py

Last active May 3, 2018 08:52
Show Gist options
• Save komasaru/b8fca575ddab23c56d94a9155446f7ea to your computer and use it in GitHub Desktop.
Python script to convert coordinates.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #! /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)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 """ 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 """ 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters