파이썬으로 천체력을 갖고 놀아보자

영화 마션을 본 분들은 JPL 이란 집단을 알텐데요,
제트추진 연구소에서는 천체력을 FTP 로 release 하고 있습니다.

삼체문제는 어차피 유리화된 해답을 얻는게 불가능하기 때문에
적분을 통해 상호작용들을 풀어나가야하기에 이런 저런 계산법들이 동원되고
실측값과 맞춰보며 자꾸만 개선되고 있는 것이죠.



https://www.projectpluto.com/jpl_eph.htm

위와 같은 곳에서 개념을 잡으시면 되고,
우선, 지난번에 올린 달력 코드를 조금 개선합니다.

from    datetime import datetime
import  math
#-------------------------------+----------------------------------------------------------------

class   cal:
    celes       = '甲乙丙丁戊己庚辛壬癸包焦'
    topos       = '寅卯辰巳午未申酉戌亥子丑'
    dayws       = '火水木金土日月'

    Auto        = 'A'
    Julian      = 'J'
    Gregorian   = 'G'

    secs_1m     = 60
    secs_1h     = secs_1m * 60
    secs_1d     = secs_1h * 24
    days_5m     = 153
    days_j1y    = 365.25
    days_j4y    = int( days_j1y * 4 )
    days_g1y    = 365.2425
    days_g400y  = int( days_g1y * 400 )
    years_jd2ad = 4712
    years_j400f = int( math.ceil( years_jd2ad / 400 ) * 400 )
    jd_base     = ( years_j400f - years_jd2ad ) * days_j1y - 60
    jd_base_ad  = years_jd2ad * days_j1y + 60
    j2g         = 38

    def __init__( self, jd = 0. ):
        self.jd = jd

    def __str__( self ):
        return  f'JD {self.jd:f}'

    @staticmethod
    def _time( h, n, s ):
        return  ( h * cal.secs_1h + n * cal.secs_1m + s ) / cal.secs_1d - 0.5

    @staticmethod
    def g( gy = -4712, gm = 1, d = 1, h = 12, n = 0, s = 0 ):
        y = gy - ( gm < 3 ) + cal.years_j400f
        return  cal( cal.j( gy, gm, d, h, n, s ).jd - y // 100 + y // 400 + cal.j2g )

    @staticmethod
    def j( jy = -4712, jm = 1, d = 1, h = 12, n = 0, s = 0 ):
        y = jy - ( jm < 3 )
        m = ( jm + 9 ) % 12
        return  cal( d - 1 + ( cal.days_5m * m + 2 ) // 5 +
                     y * cal.days_j4y // 4 + cal.jd_base_ad + cal._time( h, n, s ) )

    @staticmethod
    def a( y = -4712, m = 1, d = 1, h = 12, n = 0, s = 0 ):
        return  cal.j( y, m, d, h, n, s ) if y * 10000 + m * 100 + d < 15821015 else \
                cal.g( y, m, d, h, n, s )

    @staticmethod
    def _ymd( c, b = 0 ):
        y = ( 4 * c + 3 ) // cal.days_j4y
        e = c - ( cal.days_j4y * y ) // 4
        m = ( 5 * e + 2 ) // cal.days_5m

        D = e - ( cal.days_5m * m + 2 ) // 5 + 1
        M = ( m + 2 ) % 12 + 1
        Y = b + y - cal.years_j400f + ( M < 3 )
        return  Y, M, D

    @staticmethod
    def hns( jd ):
        ss  = int( ( jd + 0.5 ) * cal.secs_1d ) % cal.secs_1d
        return  ss // cal.secs_1h, ss // cal.secs_1m % 60, ss % 60

    def gdt( self ):
        a = int( self.jd + 0.5 + cal.jd_base ) - cal.j2g
        b = ( 4 * a + 3 ) // cal.days_g400y
        c = a - ( b * cal.days_g400y ) // 4
        return  self._ymd( c, b * 100 ), self.hns( self.jd )

    def jdt( self ):
        return  self._ymd( int( self.jd + 0.5 + cal.jd_base ) ), self.hns( self.jd )

    def __sub__( self, other ):
        return  self.jd - other.jd

    def __add__( self, other ):
        return  cal( self.jd + other )

    def weekday( self ):
        return  cal.dayws[ int( self.jd - 0.5 ) % 7 ]

    def kanji( self ):
        return  cal.celes[ int( self.jd - 0.5 ) % 10 ] + \
                cal.topos[ int( self.jd - 0.5 ) % 12 ]

    def str( self, mode = Auto, show_zero_heading = True ):
        if mode == cal.Auto:
            mode = 'J' if self.jd < cal.g( 1582, 10, 15, 0 ).jd else 'G'
        d, t = ( self.jdt() ) if mode == 'J' else ( self.gdt() )
        l = '02d' if show_zero_heading else 'd'
        mode += '.AD ' if d[0] > 0 else '.BC '
        mode += str( d[0] if d[0] > 0 else ( 1 - d[0] ) ) + '.'
        mode += f'{d[1]:{l}}.{d[2]:{l}}.{t[0]:{l}}:{t[1]:{l}}:{t[2]:{l}}'
        return  mode

    @staticmethod
    def leap( y, mode = Auto ):
        if mode == cal.Auto:
            mode = cal.Julian if y < 1583 else cal.Gregorian
        return  ( y % 4 == 0 ) if mode == cal.Julian else \
            ( ( y % 4 == 0 ) and ( y % 100 != 0 or y % 400 == 0 ) )

    @staticmethod
    def days_in( y, m, mode = Auto ):
        M = ( m + 9 ) % 12
        M1 = ( M + 1 )
        DS = ( cal.days_5m * M1 + 2 ) // 5 - ( cal.days_5m * M + 2 ) // 5
        return  DS + ( cal.leap( y, mode ) - 2 ) * ( m == 2 )

    @staticmethod
    def endar( y, m, mode = Auto, more = False ):
        if mode == cal.Auto:
            mode = cal.Julian if y * 100 + m < 158210 else cal.Gregorian
        h = cal.j( y, m, 1 ) if mode == cal.Julian else cal.g( y, m, 1 )
        b = -( int( h.jd + 1 ) % 7 )
        e = cal.days_in( y, m, mode )

        lines = [ '' ]
        if more:
            lines.append( 'Julian' if mode == cal.Julian else 'Gregorian' )
        if y <= 0:
            head = 'BC '
            y = -y + 1
        else: head = 'AD '
        lines.append( f'{head}{str(y):14}.{m:2d}' )
        lines.append( '日 月 火 水 木 金 土' if more else 'SU MO TU WE TH FR SA' )

        while b < e:
            msg = ''
            kmsg = [ '', '' ]
            for _ in range( 7 ):
                kanji = cal( h.jd + b ).kanji()
                kmsg[ 0 ] += kanji[ 0 ] + ' '
                kmsg[ 1 ] += kanji[ 1 ] + ' '
                b += 1
                msg += f'{b:2d} ' if 0 < b <= e else '   '
            lines.append( msg )
            if more:
                lines.append( kmsg[ 0 ] )
                lines.append( kmsg[ 1 ] )
                lines.append( '-' * 20 )
        return  lines

    @staticmethod
    def datetime( dt ):
        return  cal.a( dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second )

    @staticmethod
    def now():
        return  cal.datetime( datetime.today() )

#-------------------------------+----------------------------------------------------------------

파이써닉하게도 JPL 의 천문력을 wrapping 해 놓은 게 있습니다.

pip install jplephem

JPL 에서 제공해주는 천문력은 텍스트 인코딩된것과 바이너리 인코딩된것이 있고,
바이너리 데이타쪽이 용량이 작고 처리가 빠르니 저는 bsp 파일을 사용하겠습니다.

ftp://ssd.jpl.nasa.gov/pub/eph/planets/bsp/

여기서 받으실 수 있는데, 이미 다운로드 코드를 데모에 삽입해 뒀으니 그냥 실행만 시키셔도 됩니다.
( 전송속도가 느린 곳에선 인내심이 필요할지도요 )

천문력 사용은 이걸로 충분한데, 잠시 테스트 해보니 win7 에서는 mmap 용량에 제한이 있더군요.
그러니 de431t 같은 기원전 13000년대를 지원하는 용량 큰 녀석은
맥이나 리눅스 win10 에서 돌리시는게 정신건강에 좋을 것 같습니다.
예제에는 작은걸 올려드립니다.

어떤 녀석인지 대충 보여드리기 위해서 예제를 만들어봤는데
OpenGL 로 대충 그리다보니 이거 파이썬 패키지가 좀 지저분하더군요.
맥에서는 큰 번거로움 없이 구동되었는데 windows 에서 호환성 문제가 좀 터졌더랬습니다.
암튼 대충만 체크하고 올려드리니 관심있으시면 돌려보시면 되겠습니다.

아래와 같은 패키지 삽질들이 필요할 수 있습니다.

conda install PyOpenGL
pip install PyOpenGL PyOpenGL_accelerate
conda install -c conda-forge freeglut

import  numpy as np
from    jplephem.spk import SPK
from    ftplib import FTP as ftp
from    platform import system

from    OpenGL.GL import *
from    OpenGL.GLUT import *
from    OpenGL.GLU import *

#-------------------------------+----------------------------------------------------------------

path                            = '' if system() == 'Windows' else '../../'
bsp_name                        = 'de421.bsp' # 'de431t.bsp'
dst_bsp                         = path + bsp_name
if not os.path.exists( dst_bsp ):
    server = ftp( 'ssd.jpl.nasa.gov' )
    server.login( '', '' )
    server.cwd( '/pub/eph/planets/bsp/' )
    print( 'downloading', bsp_name )
    with open( dst_bsp, 'wb' ) as f:
        server.retrbinary( 'RETR ' + bsp_name, f.write )
    print( 'done' )

kernel                          = SPK.open( dst_bsp )

print( kernel )
print( '' )

#-------------------------------+----------------------------------------------------------------

def textout( x, y, text, color, font = GLUT_BITMAP_8_BY_13 ):
    glColor3f( *color )
    glWindowPos2f( x, y )
    for ch in text:
        glutBitmapCharacter( font, ctypes.c_int( ord( ch ) ) )

#-------------------------------+----------------------------------------------------------------

jd              = kernel[ 0, 10 ].start_jd
scale           = 100000000
label_visible   = False

#-------------------------------+----------------------------------------------------------------

class   BODY:
    def __init__( self, name, color, size, segs, gain = 1 ):
        self.obj    = gluNewQuadric()
        self.name   = name
        self.color  = color
        self.size   = size
        self.segs   = segs
        self.gain   = gain

    def interscale( self, p, d = 1 ):
        global  scale
        return  np.divide( p, scale / d )

    def position( self, jd ):
        if type( self.segs[ 0 ] ) == list:
            pos = BODY( 'parent', 0, 0, self.segs[ 0 ] ).position( jd )
        else: pos = self.interscale( self.segs[ 0 ].compute( jd ) )
        if len( self.segs ) == 2:
            pos = np.subtract( pos, self.interscale( self.segs[ 1 ].compute( jd ), self.gain ) )
        return  pos

    def draw( self, jd ):
        global label_visible
        glPushMatrix()
        glTranslatef( *self.position( jd ) )
        glColor3f( *self.color )
        gluSphere( self.obj, self.size, 21, 21 )
        if label_visible:
            x, y, _ = gluProject( self.size + 0.05, -0.07, 0 )
            textout( x, y, self.name, self.color )
        glPopMatrix()

bodies  = [
    BODY( 'Sun',     ( 1.0, 0.5, 0.2 ), 0.08,  [ kernel[ 0, 10 ] ] ),
    BODY( 'Mercury', ( 0.2, 0.2, 0.8 ), 0.02,  [ kernel[ 0,  1 ], kernel[ 1, 199 ] ] ),
    BODY( 'Venus',   ( 0.7, 0.6, 0.5 ), 0.025, [ kernel[ 0,  2 ], kernel[ 2, 299 ] ] ),
    BODY( 'Earth',   ( 0.0, 0.5, 1.0 ), 0.025, [ kernel[ 0,  3 ], kernel[ 3, 399 ] ] ),
    BODY( 'Moon',    ( 0.7, 0.7, 0.7 ), 0.01,  [ [ kernel[ 0,  3 ], kernel[ 3, 399 ] ], kernel[ 3, 301 ] ], 20 ),
    BODY( 'Mars',    ( 1.0, 0.4, 0.0 ), 0.02,  [ kernel[ 0,  4 ] ] ),
    BODY( 'Jupiter', ( 0.8, 0.9, 0.6 ), 0.06,  [ kernel[ 0,  5 ] ] ),
    BODY( 'Saturn',  ( 1.0, 0.7, 0.6 ), 0.04,  [ kernel[ 0,  6 ] ] ),
    BODY( 'Uranus',  ( 0.1, 0.0, 0.5 ), 0.05,  [ kernel[ 0,  7 ] ] ),
    BODY( 'Neptune', ( 0.0, 0.2, 1.0 ), 0.05,  [ kernel[ 0,  8 ] ] ),
    BODY( 'Pluto',   ( 0.3, 0.5, 0.3 ), 0.05,  [ kernel[ 0,  9 ] ] )
]

#-------------------------------+----------------------------------------------------------------

def init_light():
    glShadeModel( GL_SMOOTH )
    glLightfv( GL_LIGHT0, GL_SPECULAR, GLfloat_3( 1., 1., 1. ) )
    glLightfv( GL_LIGHT0, GL_DIFFUSE,  GLfloat_3( .8, .8, .8 ) )
    glLightfv( GL_LIGHT0, GL_AMBIENT,  GLfloat_3( .3, .3, .3 ) )
    glLightfv( GL_LIGHT0, GL_POSITION, GLfloat_3( 0., 0., 0. ) )

def display():
    global jd
    glClear( GL_COLOR_BUFFER_BIT )
    glColor3f( 1.0, 1.0, 1.0 )
    glLoadIdentity()
    gluLookAt( 0.0, 5.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 )
    glScalef( 1.0, 1.0, 1.0 )
    textout( 10, 10, cal( jd ).str(), ( 0, 1, 0 ), GLUT_BITMAP_9_BY_15 )

    glEnable( GL_LIGHTING )
    glEnable( GL_LIGHT0 )

    for body in bodies: body.draw( jd )

    glFlush()
    jd += 1

def reshape( w, h ):
    glViewport( 0, 0, w, h )
    glMatrixMode( GL_PROJECTION )
    glLoadIdentity()
    glFrustum( -1.0, 1.0, -1.0, 1.0, 1.5, 20.0 )
    glMatrixMode( GL_MODELVIEW )

def loop( value ):
    glutPostRedisplay()
    glutTimerFunc( 10, loop, 1 )

def keyboard( key, x, y ):
    global label_visible, scale
    if key == b'\t':
       label_visible = not label_visible
    elif key == b'q':
        sys.exit()
    elif key == b'-':
        scale *= 2
    elif key == b'=':
        scale /= 2
    else: print( key, x, y )

#-------------------------------+----------------------------------------------------------------

if __name__ == '__main__' :
    glutInit()
    glutInitWindowSize( 800, 800 )
    glutInitWindowPosition( 0, 0 )
    glutCreateWindow( b'Solar System' )

    glutInitDisplayMode( GLUT_SINGLE | GLUT_RGB | GLUT_DEPTH )
    glutDisplayFunc( display )
    glutReshapeFunc( reshape )
    glutKeyboardFunc( keyboard )

    glClearColor( 0.0, 0.0, 0.0, 0.0 )
    glEnable( GL_COLOR_MATERIAL )
    init_light()

    glutTimerFunc( 10, loop, 1 )
    glutMainLoop()

#-------------------------------+----------------------------------------------------------------

ftp 저장되는 경로는 path 에서 수정하시고,
조명 세팅은 제대로 테스트 안해봤으니 쓰시는 분들이 알아서 맞춰 보시길.

그럼~

6 Likes

이런게 어디 필요해? 하는 분들이 계시겠지만,
어차피 시간과 공간 문제를 조금 깊게 들어가면 절입일시가 중요하고, 윤초 개념도 필요한 것이거든요.
우리나라 천문연도 기상청도 JPL 의 데이타를 받아 절입일시를 뽑고 각계각층에 정보를 제공하고 있습니다.

1 Like

달은 전역 스케일로 보면 지구에 들러붙어 돌게 되기 땜에 별도의 증폭값을 적용하고 있습니다~
크기는 내맘대로 대충 적은거예요. 실제 크기를 적용하면 안보임.

1 Like

혹여 windows 에서 float128 이 없다고 라벨 출력이 실패하는 경우엔, 그냥 exception 을 raise 하는 코드만 주석 처리하시면 될겁니다.

1 Like