제가 천문에 관련된 프로그래밍을 조금 할줄 알다보니 태양의 위치 및 일출/일몰에 대한 문의가 자주 들어옵니다. 개인적인 소스는 있지만 높은 정확도를 위한 것이라 일반인들이 이해하기에는 수준이 높습니다.

출처 : http://very-bored.com/index.php?option=com_content&task=view&id=135&Itemid=29



대부분 문의유형을 살펴보면 태양전지판을 효율적으로 사용하기 위해 태양의 방향에 따라 입사각이 90도가 이뤄지도록 구동하는 장치를 개발을 해야한다던가 수족관의 전등을 일출/일몰 시간에 맞게 자동으로 On/Off시켜줘야 한다는 뭐 아주 간단하면서 정확도가 그리 중요하지 않은 수준의 구현을 요구하는 것들이 대부분 입니다.

일전에 관련되어 문의를 받았지만 요즘에는 개인적으로 너무 바쁘고 여유가 없어서 설명을 드리고 싶지만 드릴 수 있는 시간조차 없었습니다. 몇주가 답변을 못드리다가 잠깐 짬을 내어 답변을 드렸지요.

웹상에서 찾아보면 일몰/일출 시간 소스는 얼마든지 있습니다. 그 방법도 엄청 많이 나와있지요. 단지 어찌 찾아야하는 것인지 있어도 어찌 활용해야하는지 모르시는 분들이 대부분이라 생각합니다. 그래서 여기에 이 부분에 대해 설명한 사이트들을 소개하고 마지막으로 웹에 공개된 C++소스에 주석처리한 것을 공유하겠습니다.

태양위치 계산은 정확도에 따라서 계산방법은 천차만별입니다. 그러므로 완벽한 정답은 없습니다. 상황에 맞게 쓰시면 됩니다.
 
아래내용은 태양의 일출/일몰 계산 방법을 간략하게 소개하고 있습니다.
http://www.stargazing.net/kepler/sunrise.html
 
PHP 언어로 만들어진 일출/일몰 소스입니다.
http://www.phpschool.com/gnuboard4/bbs/board.php?bo_table=tipntech&wr_id=39302&sca=%C7%D4%BC%F6&page=29
 
C또는 C++로 만들어진 일출/일몰 소스입니다.
http://www.sci.fi/~benefon/stuff.html
아래 붙여놓은 소스는 여기에서 구현한 C++소스를 주석처리 해놓은 것 입니다. 참고하시고 main()함수부터 천천히 보세요. 이것으로 근사적 일출/일몰 시간을 계산할 수 있습니다.

구현하신다음 맞는지 확인하기 위해 아래 링크를 참고하세요.
http://www2009.kasi.re.kr/knowledge/solun_riset.aspx

// C++ program calculating the sunrise and sunset for
// the current date and a fixed location(latitude,longitude)
// Jarmo Lammi 1999 - 2000
// Last update January 6th, 2000

#include <iostream.h>
#include <math.h>
#include <time.h>

extern double pi;
double tpi = 2 * pi; 
double degs = 180.0/pi;
double rads = pi/180.0;

double L,g,daylen;
double SunDia = 0.53;     // Sunradius degrees(지구에서 본 태양의 반지름각. 단위:도)

double AirRefr = 34.0/60.0; // athmospheric refraction degrees (지구대기 굴절각)//

//   Get the days to J2000
//   h is UT in decimal hours
//   FNday only works between 1901 to 2099 - see Meeus chapter 7
// 이 함수는 J2000 1월 1일 12시 DT를 기준으로 날수(Day Number)를 구합니다. (Meeus의 책 챕터 7을 참고하세요)
// J2000에 대해서는 제 블로그 http://blog.jidolstar.com/487 를 참고하세요.
// 입력값은 년,월,일,시간(실수값)입니다. 
double FNday (int y, int m, int d, float h) {
int luku = - 7 * (y + (m + 9)/12)/4 + 275*m/9 + d;
// type casting necessary on PC DOS and TClite to avoid overflow
luku+= (long int)y*367;
return (double)luku - 730531.5 + h/24.0;
};

//   the function below returns an angle in the range
//   0 to 2*pi
// 
// 들어온 값을 항상 0~2pi 값안으로 normalize시킵니다. 
double FNrange (double x) {
    double b = x / tpi;
    double a = tpi * (b - (long)(b));
    if (a < 0) a = tpi + a;
    return a;
};

// Calculating the hourangle
// 시간각을 계산합니다. 
double f0(double lat, double declin) {
double fo,dfo;
// Correction: different sign at S HS
dfo = rads*(0.5*SunDia + AirRefr); if (lat < 0.0) dfo = -dfo;
fo = tan(declin + dfo) * tan(lat*rads);
if (fo>0.99999) fo=1.0; // to avoid overflow //
fo = asin(fo) + pi/2.0;
return fo;
};

// Calculating the hourangle for twilight times
// 박명시간에 대한 시간각을 계산합니다. 
double f1(double lat, double declin) {
double fi,df1;
// Correction: different sign at S HS
df1 = rads * 6.0; if (lat < 0.0) df1 = -df1;
fi = tan(declin + df1) * tan(lat*rads);
if (fi>0.99999) fi=1.0; // to avoid overflow //
fi = asin(fi) + pi/2.0;
return fi;
};


//   Find the ecliptic longitude of the Sun
//  태양의 황경을 구합니다.
double FNsun (double d) {

//   mean longitude of the Sun (태양의 평균 황경)

L = FNrange(280.461 * rads + .9856474 * rads * d);

//   mean anomaly of the Sun (태양의 평균근점이각)

g = FNrange(357.528 * rads + .9856003 * rads * d);

//   Ecliptic longitude of the Sun (태양의 황경계산)

return FNrange(L + 1.915 * rads * sin(g) + .02 * rads * sin(2 * g));
};

// Display decimal hours in hours and minutes
// 이것은 실수(Real Number)로 표현된 시간을 시:분으로 표기 하기 위한 겁니다. 
void showhrmn(double dhr) {
int hr,mn;
hr=(int) dhr;
mn = (dhr - (double) hr)*60;
if (hr < 10) cout << '0';
cout << hr << ':';
if (mn < 10) cout << '0';
cout << mn;
};


// 메인함수입니다. 
int main(void){
double y,m,day,h,latit,longit;

time_t sekunnit;
struct tm *p;

//  get the date and time from the user
// read system date and extract the year
// 사용자의 시스템에서부터 날짜와 시간을 얻습니다.
/** First get time **/
time(&sekunnit);

/** Next get localtime **/

 p=localtime(&sekunnit);

 y = p->tm_year; //년 
 // this is Y2K compliant method
 y+= 1900;
 m = p->tm_mon + 1; //월 

 day = p->tm_mday; //일

 h = 12; //낮 12시를 기준으로 계산합니다.

//관측장소의 경도(단위 degree),위도(단위 degree),타임존(단위시간) 값을 입력받습니다.
double tzone=2.0;
cout << "Input latitude, longitude and timezone and month\n";
cin >> latit;
cin >> longit;
cin >> tzone;

// testing
// m=6; day=10;

double d = FNday(y, m, day, h);

//   Use FNsun to find the ecliptic longitude of the
//   Sun
// 태양의 황경 
double lambda = FNsun(d);

//   Obliquity of the ecliptic (황도기울기 계산)

double obliq = 23.439 * rads - .0000004 * rads * d;

//   Find the RA and DEC of the Sun (태양의 적경,적위 계산)

double alpha = atan2(cos(obliq) * sin(lambda), cos(lambda)); //태양의 적경
double delta = asin(sin(obliq) * sin(lambda)); //태양의 적위 

// Find the Equation of Time
// in minutes
// Correction suggested by David Smith
// 균시차 계산 
double LL = L - alpha;
if (L < pi) LL += tpi;
double equation = 1440.0 * (1.0 - LL / tpi);
double ha = f0(latit,delta);
double hb = f1(latit,delta);
double twx = hb - ha;	// length of twilight in radians (라디안 단위 박명길이)
twx = 12.0*twx/pi;		// length of twilight in hours (시간단위 박명길이)
cout << "ha=" << ha << "  hb=" << hb << endl; 
// Conversion of angle to hours and minutes (하루길이를 시간단위로 계산) //
daylen = degs*ha/7.5;
     if (daylen<0.0001) {daylen = 0.0;}
// arctic winter     //

double riset = 12.0 - 12.0 * ha/pi + tzone - longit/15.0 + equation/60.0; //뜨는시간
double settm = 12.0 + 12.0 * ha/pi + tzone - longit/15.0 + equation/60.0; //지는시간 
double noont = riset + 12.0 * ha/pi; //정오 시간 

//정오시에 태양의 고도(maximum altitude)
double altmax = 90.0 + delta * degs - latit; 
// Correction for S HS suggested by David Smith
// to express altitude as degrees from the N horizon
if (latit < delta * degs) altmax = 180.0 - altmax;

double twam = riset - twx;	// morning twilight begin 시민박명(아침)
double twpm = settm + twx;	// evening twilight end 시민박명(저녁)

if (riset > 24.0) riset-= 24.0;
if (settm > 24.0) settm-= 24.0;

cout << "\n Sunrise and set\n";
cout << "===============\n";

cout.setf(ios::fixed);
cout.precision(0);
cout << "  year  : " << y << '\n'; //년
cout << "  month : " << m << '\n'; //월
cout << "  day   : " << day << "\n\n"; //일
cout << "Days until Y2K :  " << d << '\n'; //J2000을 기준으로 센 날수 
cout.precision(2);
cout << "Latitude :  " << latit << ", longitude:  " << longit << '\n'; //관측자의 경도와 위도 출력 
cout << "Timezone :  " << tzone << "\n\n"; //타임존 출력 (한국은 UT+9h)
cout << "Declination   : " << delta * degs << '\n'; //태양의 적위
cout << "Daylength     : "; showhrmn(daylen); cout << " hours \n"; //하루길이 

//시민박명(아침) 
cout << "Civil twilight: ";
showhrmn(twam); cout << '\n';
//일출시간 
cout << "Sunrise       : "; //
showhrmn(riset); cout << '\n';

 //정오때 태양 고도(Amendment by D. Smith)
cout << "Sun altitude at noontime "; 
showhrmn(noont); cout << " = " << altmax << " degrees" 
    << (latit>=0.0 ? " S" : " N") << endl;

//일몰시간 
cout << "Sunset        : ";
showhrmn(settm); cout << '\n';
//시민박명(저녁)
cout << "Civil twilight: ";
showhrmn(twpm); cout << '\n';

return 0;
}

어떤가요? 조금만 공을 들이면 완벽히 이해는 못해도 구현은 얼마든지 가능합니다.

글쓴이 : 지돌스타 (http://blog.jidolstar.com/705)

CDS는 프랑스의 스트라스부르 관측소 데이터 센터(Centre de Données astronomiques de Strasbourg)에서 제공하는 천체목록 데이타 서비스이다. 여기서는 VizieR(카탈로그), Simbad(천체대상), Aladin(Sky atlas)를 기본해서 다양한 천체목록 서비스를 하고 있다. 이들에 대한 가이드문서를 참고하면 이용하는데 도움이 된다.

CDS는 학술적 내용의 접근에 용이하도록 되어 있기 때문에 일반인들보다는 천문학자들에게 매우 유용하다. 본인은 이곳을 통해 다양한 천체정보를 종종 얻어오곤 한다.

CDS와 같은 서비스가 있는 이유는 천문학에서 다루는 천체들을 정리한 목록과 이름이 정말 많기 때문이다. 같은 대상이라 하더라도 그것들을 정리한 목록의 수가 많기 때문이다. 가령, 안드로메다은하만 하더라도 M31, NGC 224, Andromeda Galaxy, UGC 454, PGC 2557, LEDA 2557(참고:Wikipedia) 등 정말 많다. 이렇게 엄청 많은 목록이 존재할 수 밖에 없었던 이유는 천체목록이 생성된 시기에 따른 기술적 진보가 있었기 때문이다. 여전히 지금도 수많은 천체대상이 많이 발견되고 있기 때문에 이러한 목록은 없어지지 않고 항상 현존하게 된다.

CDS에서 제공하고 있는 목록만 하더라도 엄청 많다. (참고 : Catalogues and files available at CDS) 은하, 성운/성단, 별의 단순 분류에 따른 목록도 있지만 학술적 목적에 의해 만들어진 목록도 꽤 된다는 것을 알 수 있다.

CDS에서는 많은 개발자들의 참여를 유도할 수 있는 서비스가 준비되어 있다.

이 개발자 코너를 통해 CDS에 제공하는 데이터를 다양한 형태로 가공해 다른 서비스를 만들 수 있다.

본인도 천체대상 검색이 필요했다. 사용자가 찾기 원하는 천체의 이름을 입력하면 그에 대한 정보를 얻어오는 기능이다. 그러나 엄청많은 이름의 천체목록을 다 검색하도록 만드는 것은 왠지 불가능해보인다. (언제 데이터베이스를 만들고 그들 의 연관 알고리즘을 언제 만들어!) CDS 개발자 코너에서 제공하고 있는 CDS XML Web Services는 이것을 가능하게 했다.

CDS XML Web Services : http://cds.u-strasbg.fr/cdsws.gml

이 서비스는 XML 형태로 SOAP(Simpe Object Access Protocal)을 이용해 서비스가 진행되며 Simbad, NED, VizieR에서 제공하는 천체목록을 검색해서 결과를 반환해준다.

위 서비스에서 Name Resolver 서비스가 있다.

Name Resolver는 천체이름을 입력하면 다른 천체이름과 천체에 대한 몇가지 이름정보를 반환해주는 역할을 한다.

현재 접근 주소는 4개가 있다. 그러므로 어느 한 주소로 접근하지 못하더라도 다른 주소로 접근하면 되겠다.


    위 주소로 접근해서 사용하는 예제도 CDS Name Resolver 서비스 문서에 매우 잘 나와 있으며 JavaPerl 예제도 있다.
     
    아래는 PHP 예제이다.

    <?php
    $wsdl_url = "http://cdsws.u-strasbg.fr/axis/services/Sesame?wsdl";
    $client = new SoapClient($wsdl_url);
    $tmp = $client->Sesame('M 31', 'x', true, "N");
    var_dump($tmp);
    ?>
    


    위 소스가 제대로 돌아가려면 php.ini 에 Soap가 동작할 수 있을 수 있게 모듈을 열어주어야 한다.  위와 같이 할 경우 M31 대상을 XML형태로 NED에서 찾는 것을 의미한다. XML 형태로 $temp에 저장되므로 별도로 XML 파서를 이용해야할 것이다. 한가지 예이긴 하지만 문서를 잘보면 사용방법을 익힐 수 있다. 읽어오는 방식도 XML뿐 아니라 HTML, String형태로도 반환하며 각종 옵션도 있으므로 꼼꼼히 챙겨보자. 단 deprecated로 표시된 것은 다음에 없어질 것이니 사용하지 말자. 위 형태 말고 다음처럼 쓸수도 있다.

    sesame("m31","x",true,"N"), sesame("m31","xi",false,"A"), sesame("m31","Hpi",true,"NS"), ...


    위 결과는 다음과 같다.


    <?xml version="1.0" encoding="UTF-8"?>
    <Sesame xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
     xsi:noNamespaceSchemaLocation="http://vizier.u-strasbg.fr/xml/sesame_4.xsd">
    <Target option="N">
      <name>M 31</name>
      <!-- Q5133035 #1 -->
      <Resolver name="N=NED">
        <INFO>from cache</INFO>
        <otype>G</otype>
        <jpos>00:42:44.32 +41:16:08.5</jpos>
        <jradeg>010.6846833</jradeg>
        <jdedeg>+41.2690361</jdedeg>
        <refPos>1992ApJ...390L...9C</refPos>
        <errRAmas>375</errRAmas><errDEmas>375</errDEmas>
        <Vel><v>-300.09229</v><e>3.89730</e><r>1991RC3.9.C...0000d</r></Vel>
        <MType>SA(s)b</MType>
        <oname>MESSIER 031 </oname>
        <alias>NGC 0224</alias>
        <alias>Andromeda Galaxy</alias>
        <alias>UGC 00454</alias>
        <alias>CGCG 535-017</alias>
        <alias>CGCG 0040.0+4100</alias>
        <alias>MCG +07-02-016</alias>
        <alias>GIN 801</alias>
        <alias>B3 0040+409</alias>
        <alias>2MASX J00424433+4116074</alias>
        <alias>IRAS  00400+4059</alias>
        <alias>IRAS F00400+4059</alias>
        <alias>KTG 01C</alias>
        <alias>HOLM 017A</alias>
        <alias>PGC 002557</alias>
        <alias>UZC J004244.3+411608</alias>
        <alias>87GB 004002.2+405940</alias>
        <alias>87GB[BWE91] 0040+4059</alias>
        <alias>6C B004001.6+410004</alias>
        <alias>MY 0040+409A</alias>
        <alias>CXOM31 J004244.3+411608</alias>
        <alias>RX J0042.6+4115</alias>
        <alias>1RXS J004241.8+411535</alias>
        <alias>EXSS 0039.9+4059</alias>
        <alias>1H 0039+408</alias>
        <alias>1ES 0039+409</alias>
        <alias>XSS J00425+4102</alias>
        <alias>LGG 011:[G93] 001</alias>
        <alias>[PFJ93] 44</alias>
        <alias>[MHH96] J004241+411531</alias>
        <alias>[VCV2001] J004244.3+411610</alias>
        <alias>MESSIER 031:[KGP2002] r1-010</alias>
        <alias>MESSIER 031:[PFH2005] 321</alias>
        <alias>MESSIER 031:[VG2007] 001</alias>
        <alias>0039+408</alias>
        <alias>0040+4059</alias>
        <alias>LEDA 002557</alias>
        <nrefs>3044</nrefs>
      </Resolver>
    </Target>
    </Sesame>
    <!--- ====Done (2009-Dec-29,08:30:21z)==== -->
    


    결과에서 alias(별칭) 부분이 바로 검색한 M31의 다른 이름들이다. 정말 많다. ^^;
    간단히 위치정보(jpos, jradeg, jdedeg)도 알 수 있다.

    이를 이용해 다양한 정보를 제공하는 또 다른 천문학 관련 애플리케이션을 만들 수 있게 된다.
    혹시 지나가다가 이것을 이용해서 만든 프로그램이 있다면 알려주었으면 한다. ^^

    글쓴이 : 지돌스타(http://blog.jidolstar.com/641)

    오늘날 우리가 일상생활에서 사용하는 시간체계는 UTC(Coordinated Universal Time; 세계협정시)이다. 이는 원자시계를 근간으로 하지만 더불어 지구의 운동까지 고려한 UT(Universal Time;세계시)와 0.9초 차이가 안나도록 윤초(leap year)를 매년 6월 말 및 12월 말에 넣어주는 방식을 택하고 있다. 즉, UTC는 원자시계 + 지구운동을 같이 고려하고 있는 것이다. 원자시계를 기반으로하는 시간체계는 TAI(International Atomic Time, 국제원자시)이다. 이 시간은 2009년 이래로 UTC와 +64초 차이가 나게 되었다. TAI는 매우 정교한 시간을 만들어주기 때문에 느려지거나 빨라지지 않는다. 하지만 UTC는 UT와의 차이를 좁혀주므로 지구의 운동에 민감하게 된다.  UTC와 UT는 거의 구별을 안하지만 엄밀히 말하면 다르다. UT는 지구자전에 의해 연속적으로 변하는 값이지만 UTC는 UT와 0.9초 벌어지지 않기 위해 불연속적으로 변한다.

     

    천문계산을 하는데 있어서 일상에서 사용하는 시간과 과학을 목적으로하는 시간체계는 차이가 있다. 일상생활에 사용하는 UT나 UTC는 과학계산에 적합하지 않다. 왜냐하면 균일한 시간이 아니기 때문이다. TAI를 사용하기 이전에 시간의 균일하지 않은 문제를 해결하기 위해 ET(Ephemeris Time,역표시), TDT(Terrestiral Dynamical Time), DT(Dynamical Time)등을 만들게 되었다.

     

    식(eclipse, 일식/월식등)을 예견하기 위해 태양, 달의 위치는 Terrestiral Dynamical Time(DT)기반으로 계산된다. 왜냐하면 DT는 균일한 시간체계이기 때문이다. 우리가 평소에 쓰는 것은 UTC이고 이것은 UT와 차이없도록 만들어진다. 이는 지구의 운동에 관련된 시간체계이기 때문에 균일한 시간이 아니다. 그래서 일상생활의 시간에 어떤 천문 현상을 예측하기 위해 항상 UT와 DT의 차이값을 알아야만 한다. 이 차이값을 delta-T(ΔT)라고 하며 다음과 같은 관계를 가진다.

     

    ΔT = DT – UT = 32.184 + (TAI-UTC) - (UT1-UTC)

     

    이 식을 보면 알겠지만 ΔT는 DT와 UT의 차이면서 TAI, UTC, UT와 관계가 있다는 것을 알 수 있다. 즉 DT를 몰라도 UT, TAI, UTC만 알면 충분히 계산해 낼 수 있는 시간이다. UT의 경우 지구의 운동을 측정하면 되고 TAI는 원자시계이므로 이미 알고 있는 값이다. 또한 UTC도 윤초를 더한 날이 언제인지 알고 있으므로 충분히 알아낼 수 있는 값이다. 하지만 이런 관계는 현대시간에서만 이런 계산이 가능하다. 일반적으로 ΔT를 계산하는데는 과거, 현재, 미래까지 통용할 수 있는 식을 원한다. 그럼 먼저 과거의 ΔT는 어떻게 구할까?

     

    과거의 ΔT는 역사적 기록으로부터 추론된다. 일찍이 우리나라를 비롯한 유럽, 중동, 중국등의 나라에서는 수십건의 식(eclipse)을 관측하고 기록해왔다. 역사적 기록은 낮은 정확성에도 불구하고 과거 ΔT를 계산하는데 중요한 자료로 사용된다. 1609년부터는 망원경을 이용한 천체관측을 시작하게 되어 달의 별을 가림현상(달의 엄폐)등을 관측하여 전보다 높은 정확도를 가진 관측자료를 만들게 되어 ΔT의 오차를  줄일 수 있게 되었다.

    현대에는 지구 운동과 독립적인 원자시계 및 퀘이사 전파측정으로 거의 완벽한  ΔT를 계산할 수 있게 되었다. 하지만 지구의 운동은  어떻게 변할지 예상하기 힘들다. 1965년부터 1980년까지 ΔT는 평균적으로 1년에 0.99초씩 증가했다. 또 1985년부터 2000년까지 ΔT는 평균적으로 1년에  0.63초씩 증가했다. 2000년부터 2005년까지는 단지 0.18초 증가했다. 이렇게 불규칙하게 변하는 ΔT이기 때문에 미래의 ΔT의 정확한 예측은 이론적으로는 힘들다. 대신 달에 의한 조수차로 지구의 자전 주기가 길어지는 것은 예상이 가능하다. 예측 계산에 따르면 2010년에는 ΔT는 +67초, 2050년에는 +93초, 2100년에는 203초, 그리고 2200년에는 442초가 될 것이라 예상할 수 있다.

     

    ΔT를 계산하는 방법은 매우 다양하다. 본인은 이 ΔT를 계산하기 위해 NASA(미국항공우주국)에서 제공하는 “다항식을 이용한 ΔT 계산”을 이용해서 B.C. 2000년 부터 A.D. 3000년까지 ΔT를 계산해보았다.

     

    다음 프로그램에서 버튼을 누르면 ΔT를 계산해서 출력해준다. 누르고 반응이 없더라도 기다리면 출력된다.

    여기서 사용한 식은 Morrison과 Setphenson[2004]가 연구결과로 만든 것으로 달의 평균운동으로 인해 결정된 조석항(the secular tidal term in the mean motion of the moon)을 -26 arcsec/cy^2로 가정하고 계산한 것이다. 하지만 ELP-2000/82 lunar ephemeris에서 사용하는 달의 식변화는 이 값을 -25.858 arcsec/cy^2((Chapront, Chapront-Touzé, and Francou,2002)으로 계산했다. 그러므로 필요한 경우 최종 계산된 결과값인 ΔT에 아래와 같은 작은 보정값 c를 구해 더해줄 필요가 있다.  참고로 cy는Julian century를 뜻한다.

     

    c = - 0.91072 * (ndot + 26.0 ) * t^2
    t = (year-1955)/100
    ndot=-25.858

     

    위 보정값 c를 아래처럼 ΔT에 더한다.

     

    보정된 ΔT = ΔT + c

     

    위 ndot는 다른 값을 넣을 수 있다. 필요에 따라서 보정값을 조정하면 되겠다.

     

    앞서 설명했지만 ΔT는 정확하게 구할 수 있는 값이 아니다. 연구하는 사람마다 다르게 나오며 시간이 지날 수록 지구의 운동을 측정을 통해 계속 보정하여 다른 식으로 바뀌어질 수 있다. 그때는 이 식을 사용하는 사람이 알아서 기존에 사용하던 식을 대체할 필요가 있을 것이다.

     

    해당 년도에 대해 역사적으로 추론 또는 근래 계산법에 의한 ΔT의 오차범위에 대한 정보는 “Historycal values of DELTA T(ΔT)”를 참고한다.

     

    참고글

     

     글쓴이 : 지돌스타(http://blog.jidolstar.com/486)

    세차운동(歲差運動, Precession)은 일반적으로 강체에 돌림힘(Torque)이 작용할 때, 회전하는 물체의 축이 어떤 부동축의 둘레를 회전하는 현상을 말한다. 천문학에서 말하는 세차운동은 지구의 자전축이 불변인 황도면의 축의 둘레를 2만 5800년 주기로 회전하는 운동을 말한다. 이 현상은 지구에 작용하는 달, 태양의 중력과 지구가 타원체라는 것에서부터 기인한다.

     

    천문 계산을 하는데 있어서 세차를 고려해야 과거,미래의 천체위치를 제대로 예측할 수 있다. 왜냐하면 세차운동으로 인해 지구 자전축이 변한다는 것은 천체 위치의 기준점인 춘분점(Vernal Equinox)의 위치도 변한다는 것을 의미하기 때문이다.

     

    여기서는 적도좌표계(Equatorial Coordinate System)황도좌표계(Ecliptic Coordinate System)에서 세차운동 계산 방법을 학습하는 것을 주 목표로 한다. 세차운동의 공식 유도 및 사용되는 각종 다항식이 어떻게 만들어졌는지 알아내는 과정은 생략한다.

     

    이 글의 참고 서적은 Jean Meeus의 Astronomical Algorithms 2nd임을 밝혀둔다. 하지만 이 책에서 소개한 대로만 하지 않으며 더 효과적으로 계산하는 방법도 소개하겠다.

     

    이 글은 세차운동 계산시 정밀한 방법(Rigorous Method)을 이용한다. 낮은 정밀도를 가지는 방법은 책을 참고하길 바란다.

     

    Julian Century

    1 Julian Century는 100 Julian Year와 같다.  Julian Year는 86000 SI 초를 기반으로하는 365.25일을 1로 측정하는 단위이다. Julian Century는 Julian Year의 100배 이므로 36525일을 1로 측정하는 단위이다.  보통 J2000.0 = JDE 2451545.0 = 2000년 1월 1일 12h DT를 0으로 시작하게 된다. 여기서 JDE는 Julian Ephemeris Date, DT는 Dynamic Time을 의미한다.

     

    그러므로 Julian Date(JD)로부터 Julian Century는 다음과 같이 계산한다.

     

    T = (JD0 – J2000) / 36525

     

    세차운동의 경우 어떤 시작원기(the initial epoch)부터 마지막 원기(the final epoch)까지의 세차운동값을 계산해야한다. 그러므로 시작원기와 마지막원기의 Julian Century차이는 다음과 같다.

     

    t = (JD – JD0) / 36525

    T, t에서 JD0는 시작원기, JD는 마지막 원기를 의미한다.

     

     

    적도좌표계에서 세차운동적용 (일반 공식 이용)

    Meeus의 책에 보면 다음과 같은 적도좌표계에서 세차운동을 계산하기 위한 수치적 표현공식이 적혀있다.

     

    ζ = (2306”.2181 + 1”.39656T – 0”.000139 T2 ) t
          + (0”.30188 – 0”.000344T)t2 + 0”.017998 t3

    z = ζ + (0”.79280+0”.000411*T)t2 + 0”.000205t3

    ϑ =(2004”.3109 – 0”.85330T + 0”.000217T2)t
           - ( 0”.42665+0”.000217T )t2 + 0.041833t3

     

    이 식은 IAU(International Astronomical Union, 1976)에서 지정한 값으로 관측을 통해 얻어지는 데이타로 만들어지는 식이다. 참고로 IAU에서는 IAU2006 세차모델을 사용하도록 권하고 있고 미해군천문대에서는 2009년부터 천문역서에 적용하고 있다.(참고 : IAU2006 세차모델)

     

    시작원기의 적경,적위를 (α0, δ0)라고 하고 세차운동이 적용된 마지막원기의 적경,적위를 (α, δ)로 한다면 위 수치적 표현공식을 이용해 다음과 같이 (α, δ)를 계산할 수 있다. 이 식을 유도하는 과정은 구면천문학에 관련된 책을 참고하기 바란다.

     

    A = cos(δ0)  sin( α0 + ζ )

    B = cos( ϑ ) cos( δ0 ) cos( α0 + ζ ) – sin( ϑ ) sin( δ0 )

    C = sin( ϑ ) cos ( δ0 ) cos( α0 + ζ ) + cos( ϑ ) sin( δ0 )

    tan(  α – z ) = A/B

    sin( δ ) = C

     

    이 식만으로도 정확한 계산을 할 수 있다. 그러나 컴퓨터 프로그램을 이용하여 계산을 한다면 이 식은 매우 비효율적인 식이 된다. 가령, 수천, 수만개의 별이 있다고 하자. 그 별들은 어떤 시작원기를 가지는 각각 다른 α0, δ0값을 가지고 있을 것이다. 이들 α0, δ0을 세차운동이 적용된 새로운 α, δ를 계산하기 위해 일일히 cos, sin을 계산한다. 컴퓨터 프로그램을 한 사람은 알겠지만 cos, sin은 컴퓨터의 많은 자원을 소비하는 급수형태의 계산식이다. cos, sin을 너무 많이 쓰면 계산의 효율이 급격하게 떨어져서 1분이면 도출할 수 있는 계산을 1시간~2시간 이상식 소비할 수도 있게 된다.

     

    컴퓨터 프로그램을 이용한 계산을 하려면 되도록이면 cos, sin과 같이 매우 비싼 자원을 소비하는 함수 사용을 지양하도록 해야한다. 그래야 빠른 계산을 할 수 있기 때문이다.

     

     

    적도좌표계에서 세차운동적용 ( 변환행렬이용 )

     

    앞서 세차운동이 적용된 새로운 α, δ를 계산하는 식은 한개의 천체 위치를 계산할 때마다 α0, δ0이 포함된 cos, sin을 계산해야한다. 그래서 새로운 α0, δ0가 주어져도 cos, sin을 최소한으로 이용하여 계산하는 방법을 여기서 소개한다. 그것은 변환행렬(Transformation Matrix)를 이용하는 방법이다.

     

    변환행렬을 이용하면 다음과 같은 계산이 가능해진다. 세차운동을 적용하는 3x3 변환행렬을 P라고 하자. 그리고 α0, δ0를 직교좌표계로 바꾼 백터(vector)를 r0=[x0,y0,z0]라고 하고 세차운동이 적용된 α, δ를 직교좌표계를 바꾼 백터를 r=[x,y,z]라고 한다. 그럼 다음과 같이 계산하면 되겠다.

     

     

    변환행렬 P는 세차운동이 적용되어 앞서 설명한 공식과 다르게 α0, δ0이 분리가 되었다. 행렬 P만 구하면 cos, sin과 같이 큰 자원을 소비하는 함수를 계산할 때마다 이용할 필요가 없어지며 대신, 단순하게 행렬과 백터 곱으로 수천,수만개의 별을 빠른 시간안에 세차운동을 적용할 수 있게 된다.

     

    그럼 P는 어떻게 구할까? 그것은 회전변환행렬을 이용하면 된다. 앞서 구한 ζ, z, ϑ 값은 모두 세차운동에 의해 계산된 각도이다. 이 각도를 이용해 시작원기때의 적도좌표를 3번 회전변환하여 마지막원기때의 적도좌표를 구하는 것이다.

     

    회전변환행렬에서 X축 회전 행렬을 Rx, Y축 회전 행렬을 Ry, Z축 회전행렬을 Rz라고 하자. 이들 행렬은 모두 3x3 행렬이다.

     

     



     

     

    다른 참고를 보면 행렬안에 sin의 부호가 바뀐것을 볼 수 있는데 회전의 대상을 어떻게 정하느냐에 따라서 다른 결과이다. 천문계산에서 회전행렬은 위 정의대로 하는 것이 좋다.

     

    적도좌표계에서 세차운동을 적용한 변환행렬 P는 다음과 같이 주어진다.

     

     

     

    위 변환행렬은 ζ, z, ϑ 이 적용된 행렬곱이다. 이 행렬곱의 결과는 다음과 같다.

     

     

    컴퓨터 프로그램으로 계산할 때는 이렇게 변환행렬을 만들어 계산하면 되겠다.

     

     

    황도좌표계에서 세차운동적용 ( 변환행렬이용 )

     

    황도좌표계에서도 적도좌표계와 동일한 방식이 적용된다.

     

    Π = 174˚.876383889 + 3289”.4789 + 0”.60622 T2 
          - (869”.8089 + 0”.50491T)t + 0”.03536 t2

    π = (5029”.0966 + 2”.22226T – 0”.000042T2)t
         + (1”.11113 – 0”.000042T)t2 – 0”.000006t3

    η = (47”.0029 – 0”.06603T – 0”.000598T2)t
         + (-0”.03302 + 0”.000598T)t2 + 0.000060t3

     

    시작원기 황경(ecliptic longitude),황위(ecliptic latitude)가 λ0,β0 이고 마지막원기 λ,β이라고 하자. 이때 각각의 직교좌표값을 r0, r이라고 한다면 변환행렬 P는 다음과 같다.

     


     

    결과적으로 P는 다음과 같이 결정된다.

     

     

     

    응용하기

    지금까지 세차운동이 적용된 2가지 변환행렬을 알아보았다. 이 행렬은 필요에 따라서 선택해서 사용한다. 한가지 예를 들어보겠다. 적도좌표계에서 세차운동을 적용한 변환행렬을 Pequ, 황도좌표계에서 세차운동을 적용은 변환행렬을 Pecl이라고 하자. 황도좌표계에서 적도좌표계로 변환하는 행렬을 EclToEqu라고 가정한다. 만약 어떤 행성의 황경과 황위를 계산했다고 하면 다음과 같은 2가지 방법으로 새로운 적경, 적위를 계산할 수 있을 것이다.

     

    r = EclToEqu Pecl  r0

    r = Pequ Ecl2Equ  r0

    여기서 r0는 황경,황위를 담은 직각좌표계에서의 벡터값, r은 적경,적위를 담은 직각좌표계에서의 벡터값

     

    세차운동 계산기(Precession Calculator)

    이 프로그램은 지금까지 설명한 내용을 토대로 만든 프로그램이다.  

     

     

    참고


    글쓴이 : 지돌스타(http://blog.jidolstar.com/487)

    여기서 다루는 문제는 천문 관련 계산할 때 반드시 숙지하고 있어야할 각도 변환에 관련된 것이다. 매우 기초적인 내용이다. 꼭 천문 계산이 아니더라도 수학계산을 위해 이러한 내용은 잘 알고 있어야 한다.

     

    혹시 다음 값을 어떻게 읽는가?

     

    • 2h 44m 11.986s
    • 49°13'42.48"

     

    아마도 대부분 사람들이 다음처럼 읽을 것이다.

     

    • 2시 44분 11.986초
    • 49도 13분 42.48초

     

    2h 44m 11.986s 는 시(hour)이며, 49°13'42.48"는 도(degree)이다.

     

    잘 알고 있는 사실이지만 24시간은 360도에 대응한다. 그리고 1시간은 15°에 대응한다. 또한 시간에서의 1분은 도에서 15분과 대응한다. 문제는 여기서 발생한다. 한국어로 번역할때 시간의 1분과 도의 1분은 명확히 다른 단위라는 것이다. 영문으로는 시간의 분은 minute를 쓰고 도에서의 분은 arcminute를 사용한다.

     

    이러한 시간과 도의 차이를 잘 모르고 용어를 섞어쓰면 나중에 분명히 문제가 발생한다. 그래서 도와 시간의 표현에는 영어를 쓰는 것이 안전하다. 다음처럼 쓰자.

     

    • 2h = 2 hours 또는 2h
    • 44 m = 2 minutes 또는 2m
    • 11.986s = 11.986 seconds 또는 2 s
    • 49° = 49 degrees 또는 49 deg
    • 13' = 13 arcminutes 또는 13 arcmin
    • 42.48" = 42.48 arcseconds 또는 42.48 arcsec

     

    아래값을 보자.

     

    • 1h 21m 23.2s
    • 1°21'23.2"

     

    위 값은 둘다 1.3564 값을 가진다.(둘다 분에는 60을 곱하고 초에는 3600을 곱해서 시나 도에 더하면 된다.) 하지만 단위가 하나는 hours이고 하나는 degrees이다. 1h 21m 23.2s 를 각도로 표시하기 위해서는 15가 곱해져야 한다. 결과적으로 1.3564 h * 15°/h = 20.346°가 된다.

     

    앞서 설명한 도와 시간은 사람이 보기에 편한 단위이다. 누구나 360도와 24시에 대해서 복잡하다는 느낌을 받는 사람은 없을 것이다. 하지만 수학적 계산을 할 때는 이러한 값을 사용할 수 없다. 수학적 계산을 위해서는 각도와 시간은 항상 라디안(radian)값으로 변환해야한다. 지름 1의 원의 둘레값으로 π=3.141592653… 이다. 이 값이 radian이다. 즉 지름과 원의 둘레값의 단위를 동일하게 봄으로써 수학적 계산을 가능하게 한다. 360도와 24시는 2π=6.2831853…와 대응한다.

     

    각도를 표현하는 단위는 지금까지 살펴본 결과 도(degree), 시(hour), radian 이라는 것을 알았다. 천문학에서 이것 외에는 사용하지 않는다. 가령, 적경(Right Asension),적위(Declination)은 각각 시와 도를 사용한다. 앞서 설명했듯이 시와 도는 수학계산에는 적합하지 않으므로 radian으로 변환해야한다. 가령, cos, sin과 같은 수학함수를 이용할 때 도, 시 단위로는 계산할 수 없으므로 반드시 radian으로 변환한 다음에 사용해야한다.

     

    천문 프로그래밍을 한다면 각도 단위 변환을 하기 위해 다음과 같은 방법으로 변환값을 미리 만들어 사용한다.

     

    /**
     * Pi
     */
    public static const PI:Number = 3.141592653589793238462643383279502884197;
    
    /**
     * Pi * 2
     */
    public static const TPI:Number	= 6.2831853071795864769252867665590;
    
    /**
     * Half of PI
     */
    public static const HPI:Number = 1.5707963267948966192313216916398;
    
    /**
     * Math.sqrt(2);
     */
    public static const SQRT2:Number = 1.5707963267948966192313216916395;
    
    /**
     * Radian -> Degree
     */
    public static const R2D:Number = 57.295779513082320876798154814105;
    
    /**
     * Degree -> Radian
     */
    public static const D2R:Number = 0.017453292519943295769236907684886;
    
    /**
     * Arcsecond -> Radian
     */
    public static const S2R:Number = 4.8481368110953599359e-6;
    
    /**
     * Radian -> Arcsecond
     */
    public static const R2S:Number = 206264.80624709635515647335733078; //3600 * 180/PI
    
    /**
     * Arcminute -> Radian
     */
    public static const M2R:Number = 2.9088820866572159615394846141477e-4; // 1/60 * PI /180
    
    /**
     * Radian -> Arcminute
     */
    public static const R2M:Number = 3437.7467707849392526078892888463; //60 * 180/PI
    
    /**
     * Arcminute -> degree
     */
    public static const M2D:Number = 1/60;
    
    /**
     * degree -> Arcminute
     */
    public static const D2M:Number = 60;
    
    /**
     * Radian -> Hours
     */
    public static const R2H:Number = 3.8197186342054880584532103209403;
    
    /**
     * Hours -> Radian
     */
    public static const H2R:Number = 0.26179938779914943653855361527329;
    
    /**
     * Degree -> Hour
     */
    public static const D2H:Number = 1/15;
    
    /**
     * Hour -> Degree
     */
    public static const H2D:Number = 15;

     

     

    만약 시(hour)에서 도(degree)로 변환하려면 14h * H2D = 210 deg 가 된다. 이러한 방법을 이용해여 수학적 계산을 마무리 짓고 나온 결과값이 각도인 경우 그에 맞게 시(hour) 또는 도(degree)로 변경해주어 사람이 보기 편리하도록 만들어주면 되는 것이다.

     

     

    + Recent posts