天文坐标系转换

天文坐标系转换

• 黄道坐标系
• 赤道坐标系
• 银道坐标系
• 坐标系转换
• 实现坐标系转换代码

黄道坐标系

黄道坐标系，是以黄道作基准平面的天球坐标系统，多用作研究太阳系天体运动情况之用。黄道是地球上的人看太阳于一年内在恒星之间所走的视路径，即地球的公转轨道平面和天球相交的大圆。黄道和天赤道成23度26分的角，相交于春分点和秋分点。黄道即是太阳周年视运动轨道，通俗来说，由于地球上的人通常感觉不到地球的运动（公转），就像坐在行驶车辆中的人感觉的是周围的物体向后运动一样，地球人所看到的是太阳在恒星组成的星空背景上向后运动，每年转一圈，并将其称为太阳周年视运动，将太阳运行线路（即地球公转轨道在天球上的反映）称为黄道。

由于地球的自转轴没有垂直于轨道平面，所以赤道平面不与黄道平行，而有23°26’的夹角，这就是所知的黄赤交角。赤道平面和黄道平面与地球的交集所形成的大圆分别称为赤道和黄道，这两个平面的交叉点正好在一条地球直径线的两端点，就是著名的二分点(春分点与秋分点)。太阳从南向北经过的二分点称为春分点或是白羊座第一点，黄道经度，通常以字母λ标示，就以这一点为起点向东从0°到360°。

黄道纬度，通常以字母β标示，以黄道为测量的基础平面向北从0°到90°，向南从0°到-90°。春分点同样的也被定义为赤道坐标的原点，赤经的测量也是向东由0到24时，通常以字母 α \alpha α或R.A.表示；赤纬以字母δ或Del.表示，由赤道平面向北从0°到90°，向南从0°到-90°。简单的转动型式可以让 α \alpha α，δ和 λ,β互相转换。

赤道坐标系

地球赤道平面延伸后与天球相交的大圆，称为天赤道。天赤道的几何极称为天极。天赤道是赤道坐标系中的基圈，北天极P 是赤道坐标系的基本极。由于所取的主圈﹑主点以及随之而来的第二坐标的不同，赤道坐标系又有第一赤道坐标系和第二赤道坐标系之分。

第一赤道坐标系的主圈是子午圈，主点是天赤道与子午圈在地平圈之上的交点F ，天体的第二坐标是大圆弧FB =t 或球面角FP ，t 称为天体的时角。由主点F 开始按顺时针方向量度时角t ，从0°～360°，或从0～24。周日运动不会改变天体的赤纬，而仅仅使时角发生变化。

第二赤道坐标系的主点是春分点，它是黄道对赤道的升交点，过春分点的赤经圈就是该坐标系的主圈，春分点的时圈与天体时圈之间的球面角 P 或大圆弧B =α，是天体在第二赤道坐标系中的第二坐标，称为天体的赤经，赤经α是由春分点开始按逆时针方向量度的，从0°～360°，或从0～24。第一赤道坐标系是右旋坐标系，第二赤道坐标系为左旋坐标系。

天体的周日运动不影响春分点与天体之间的相对位置，因此也就不会改变天体的赤经和赤纬，而在不同的测站、不同的观测时间，天体的时角)却是变化的。所以，在各种星表中通常列出的都是天体在第二赤道坐标系中的坐标──赤经和赤纬，供全球各地的观测者使用。

银道坐标系

在有关恒星动力学和星系结构的某些理论工作中，常常采用一种球面坐标系──银道坐标系。银河系的主要部分是一个扁平的圆盘状结构，它的平均平面称为银道面。银道面是银道坐标系的基本平面，它与天球相交的大圆称为银道，也就是银道坐标系中的基圈。天球上与银道相平行的小圆称为银纬圈。

银道的几何极称为银极，又有南﹑北银极之分。作为银道坐标系的极是北银极L ，过两个银极所作的半个大圆称为银经圈，也就是银道坐标系中的副圈。所有的银经圈都与银道相垂直。

银道与天赤道在天球上相交于两点。由北银极向银道面看去，按逆时针方向从赤道以南向北通过赤道的那一个点，称为银道对天赤道的升交点﹔另一点就是降交点。1958年以前，采用银道升交点作为银道坐标系的主点，过该点的银经圈就是这一坐标系的主圈。

天体的银经圈与银道交于B 点，大圆弧B =b 就是天体在银道坐标系中的第一坐标，称为银纬。由银道起沿银经圈向南北银极分别量度银纬b ，从0°～±90°，南银纬取为负值。过升交点的银经圈与天体的银经圈所交的球面角L 或银道上的大圆弧B ，是天体在银道坐标系中的第二坐标，称为银经。1958年以前，银经由升交点起算，从0°～360°。量度方向是逆时针的，银道坐标系也是一种左旋坐标系。

坐标系转换

• λ \lambda λ和 β \beta β代表黄经和黄纬
• α \alpha α 和 δ \delta δ 代表赤经和赤纬
• ε = 2 3 ∘ 2 6 ′ 20.51 2 ′ ′ \varepsilon=23^{\circ}26'20.512'' ε=23∘26′20.512′′ 为黄赤交角，也就是地轴的倾角

转换公式 - λ \lambda λ, β ⟶ \beta \longrightarrow β⟶ α , δ \alpha, \delta α,δ

sin ⁡ δ = sin ⁡ ε sin ⁡ λ cos ⁡ β + cos ⁡ ε sin ⁡ β cos ⁡ α cos ⁡ δ = cos ⁡ λ cos ⁡ β sin ⁡ α cos ⁡ δ = cos ⁡ ε sin ⁡ λ cos ⁡ β − sin ⁡ ε sin ⁡ β \sin{\delta} = \sin{\varepsilon}\sin{\lambda}\cos{\beta}+\cos{\varepsilon}\sin{\beta} \\ \cos{\alpha}\cos{\delta}=\cos{\lambda}\cos{\beta} \\ \sin{\alpha}\cos{\delta} = \cos{\varepsilon}\sin{\lambda}\cos{\beta}-\sin{\varepsilon}\sin{\beta} sinδ=sinεsinλcosβ+cosεsinβcosαcosδ=cosλcosβsinαcosδ=cosεsinλcosβ−sinεsinβ

转换公式 - α , δ ⟶ \alpha, \delta \longrightarrow α,δ⟶ λ \lambda λ, β \beta β

sin ⁡ β = cos ⁡ ε sin ⁡ δ − sin ⁡ ε sin ⁡ α cos ⁡ δ cos ⁡ λ = cos ⁡ α cos ⁡ δ / cos ⁡ β sin ⁡ λ cos ⁡ β = sin ⁡ ε sin ⁡ δ + cos ⁡ ε cos ⁡ δ sin ⁡ α \sin{\beta} = \cos{\varepsilon}\sin{\delta} -\sin{\varepsilon}\sin{\alpha}\cos{\delta}\\ \cos{\lambda}=\cos{\alpha}\cos{\delta}/\cos{\beta} \\ \sin{\lambda}\cos{\beta} = \sin{\varepsilon}\sin{\delta}+\cos{\varepsilon}\cos{\delta}\sin{\alpha} sinβ=cosεsinδ−sinεsinαcosδcosλ=cosαcosδ/cosβsinλcosβ=sinεsinδ+cosεcosδsinα

实现坐标系转换代码

所需模块 math numpy
函数功能：

1. 对于dec 单位由deg转换为dms(度分秒)
2. 对于ra 单位由deg转换为dms(时分秒)
3. 赤道坐标系转为黄道坐标系
4. 黄道坐标系转为赤道坐标系
5. 赤道坐标系转为银道坐标系

import math
import numpy as np

def deg2HMS(ra='', dec='', rou=True):
    '''convert deg to ra's HMS or dec's DHS'''
    RA, DEC, rs, ds = '', '', '', ''
    if dec:
        if str(dec)[0] == '-':
            ds, dec = '-', abs(dec)
        deg = int(dec)
        decM = abs(int((dec-deg)*60))
        
        if rou:
            decS = round((abs((dec-deg)*60)-decM)*60,1)
        else:
            decS = (abs((dec-deg)*60)-decM)*60
        
        if deg ==0:
            deg ="00"
        elif deg <10:
            deg = "0%s"%deg
        
        if decM ==0:
            decM ="00"
        elif decM <10:
            decM = "0%s"%decM
            
        if decS ==0:
            decS ="00.0"
        elif decS <10:
            decS = "0%s"%decS
            
        DEC = '{0}{1}:{2}:{3}'.format(ds, deg, decM, decS)
  
    if ra:
        if str(ra)[0] == '-':
            rs, ra = '-', abs(ra)
        raH = int(ra/15)
        raM = int(((ra/15)-raH)*60)
        if rou:
            raS = round(((((ra/15)-raH)*60)-raM)*60,1)
        else:
            raS = ((((ra/15)-raH)*60)-raM)*60
            
        if raH ==0:
            raH = "00"
        elif raH <10:
            raH = "0%s"%raH
            
        if raM ==0:
            raM = "00"
        elif raM <10:
            raM = "0%s"%raM
            
        if raS ==0:
            raS = "00.0"
        elif raS <10:
            raS = "0%s"%raS
        
        RA = '{0}{1}:{2}:{3}'.format(rs, raH, raM, raS)
  
    if ra and dec:
        return (RA, DEC)
    else:
        return RA or DEC

#def HMS2deg(ra='', dec='', rou=True):

def convertEquatorial(el,eb):
    '''convert Ecliptic(el eb) to Equatorial raj decj'''
    #el Ecliptic longitude (degrees)
    #eb Ecliptic latitude (degrees)
    
    M_PI = math.pi
    sin = math.sin
    cos = math.cos
    asin = math.asin
    tan = math.tan
    atan2 = math.atan2
    
    ce = 0.91748213149438 # Cos epsilon
    se = 0.39777699580108 # Sine epsilon
    dr = M_PI/180.0
    
    sdec = sin(eb*dr)*ce + cos(eb*dr)*se*sin(el*dr)
    dec = asin(sdec) # in radians
    cos_ra = (cos(el*dr)*cos(eb*dr)/cos(dec))
    sin_ra = (sin(dec)*ce - sin(eb*dr))/(cos(dec)*se)
    ra = atan2(sin_ra,cos_ra)  # in radians
    
    #Get RA into range 0 to 360
    if ra < 0.0:
        ra = 2*M_PI+ra
    elif ra > 2*M_PI:
        ra = ra - 2*M_PI

    raj = math.degrees(ra) # convert to degrees
    decj = math.degrees(dec) #convert to degree

    ra = deg2HMS(ra=raj)
    dec = deg2HMS(dec=decj)

    return ra,dec

ra,dec = convertEquatorial(9.07039380,6.3091085)
print(ra,dec)

def convertEcliptic(raj,decj):
    '''convert raj deci to elog elat'''
    #rai in rad
    #decj in rad

    M_PI = math.pi
    sin = math.sin
    cos = math.cos
    asin = math.asin
    tan = math.tan
    atan2 = math.atan2
    
    
    deg2rad = M_PI/180.0
    epsilon = 23.439292*deg2rad
    #/*  epsilon = 23.441884*deg2rad;*/

    sinb = sin(decj)*cos(epsilon)-cos(decj)*sin(epsilon)*sin(raj)
    beta = asin(sinb)
    y = sin(raj)*cos(epsilon)+tan(decj)*sin(epsilon)
    x = cos(raj)
  
    lambdap = atan2(y,x)
    if lambdap<0:
        lambdaa = lambdap+2*M_PI
    else:
        lambdaa = lambdap

    elong = round(math.degrees(lambdaa),3) # convert to degree
    elat  = round(math.degrees(beta),3)

    return elong,elat

def convertGalactic(raj,decj):
    '''convert raj decj to gl gb'''
    # raj in rad
    # decj in rad
    M_PI = math.pi
    sin = math.sin
    cos = math.cos
    asin = math.asin
    tan = math.tan
    atan2 = math.atan2


    deg2rad = M_PI/180.0
    gpoleRAJ = 192.85*deg2rad
    gpoleDECJ = 27.116*deg2rad
    rot = [0.0]*9
    rot = np.array(rot).reshape(3,3)

    #/* Note: Galactic coordinates are defined from B1950 system - e.g. must transform from J2000.0 equatorial coordinates to IAU 1958 Galactic coords */
  
    #/* Convert to rectangular coordinates */
    rx = cos(raj)*cos(decj)
    ry = sin(raj)*cos(decj)
    rz = sin(decj)
    
    #/* Now rotate the coordinate axes to correct for the effects of precession */
    #/* These values contain the conversion between J2000 and B1950 and from B1950 to Galactic */
    rot[0][0] = -0.054875539726
    rot[0][1] = -0.873437108010
    rot[0][2] = -0.483834985808
    rot[1][0] =  0.494109453312
    rot[1][1] = -0.444829589425
    rot[1][2] =  0.746982251810
    rot[2][0] = -0.867666135858
    rot[2][1] = -0.198076386122
    rot[2][2] =  0.455983795705
    
    rx2 = rot[0][0]*rx + rot[0][1]*ry + rot[0][2]*rz
    ry2 = rot[1][0]*rx + rot[1][1]*ry + rot[1][2]*rz
    rz2 = rot[2][0]*rx + rot[2][1]*ry + rot[2][2]*rz

    #/* Convert the rectangular coordinates back to spherical coordinates */
    
    gb = asin(rz2)
    gl = atan2(ry2,rx2)
    if gl < 0:
        gl += 2.0*M_PI
    gl = round(math.degrees(gl),3)
    gb = round(math.degrees(gb),3)
        
    return gl,gb

el,eb = convertEcliptic(1.2,0.5)
gl,gb = convertGalactic(1.2,0.5)
print(el,eb)
print(gl,gb)

注，astropy模块提供了由银道坐标系转赤道坐标系功能