C++ - 平均黄道傾斜角の計算！
Updated:
平均黄道傾斜角（地球自転軸の傾き、地球公転面と赤道のなす角）の計算を C++ で行いました。
23.4度等と簡単に表すことが多いですが、実際は時々刻々と変化しております。
天文や暦等を正確に計算する際に必要になってきます。
過去には Ruby, Python で行っています。
0. 前提条件
- Debian GNU/Linux 10.6 (64bit) での作業を想定。
- GCC 9.2.0 (G++ 9.2.0) (C++17) でのコンパイルを想定。
1. 平均黄道傾斜角について
- 計算には、「暦象年表の改訂について（国立天文台）（PDF 1.7MB）」で紹介されている計算式を使用する。
2. C++ ソースコードの作成
ここでは、実行部分のみ掲載。（全てのコードは GitHub リポジトリとして公開している）
File: calc_obliquity.cpp
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
102
103
104
105
/***********************************************************
黄道傾斜角計算
DATE AUTHOR VERSION
2020.10.30 mk-mode.com 1.00 新規作成
Copyright(C) 2020 mk-mode.com All Rights Reserved.
----------------------------------------------------------
引数 : JST（日本標準時）
書式：最大23桁の数字
（先頭から、西暦年(4), 月(2), 日(2), 時(2), 分(2), 秒(2),
1秒未満(9)（小数点以下9桁（ナノ秒）まで））
無指定なら現在(システム日時)と判断。
***********************************************************/
#include "common.hpp"
#include "time.hpp"
#include "obliquity.hpp"
#include <cmath>
#include <cstdlib> // for EXIT_XXXX
#include <ctime>
#include <iostream>
#include <sstream>
#include <string>
int main(int argc, char* argv[]) {
// 定数
static constexpr double kPi = atan(1.0) * 4.0;
// 変数
std::string tm_str; // time string
unsigned int s_tm; // size of time string
int s_nsec; // size of nsec string
int ret; // return of functions
struct timespec jst; // JST
struct timespec utc; // UTC
struct tm t = {}; // for work
double jcn; // Julian Century Number
double eps; // 黄道傾斜角
namespace ns = calc_obliquity;
try {
// 日付取得
if (argc > 1) {
// コマンドライン引数より取得
tm_str = argv[1];
s_tm = tm_str.size();
if (s_tm > 23) {
std::cout << "[ERROR] Over 23-digits!" << std::endl;
return EXIT_FAILURE;
}
s_nsec = s_tm - 14;
std::istringstream is(tm_str);
is >> std::get_time(&t, "%Y%m%d%H%M%S");
jst.tv_sec = mktime(&t);
jst.tv_nsec = 0;
if (s_tm > 14) {
jst.tv_nsec = std::stod(
tm_str.substr(14, s_nsec) + std::string(9 - s_nsec, '0'));
}
} else {
// 現在日時の取得
ret = std::timespec_get(&jst, TIME_UTC);
if (ret != 1) {
std::cout << "[ERROR] Could not get now time!" << std::endl;
return EXIT_FAILURE;
}
}
// JST -> UTC
utc = ns::jst2utc(jst);
// ユリウス世紀数計算
ns::Time o_tm(utc);
jcn = o_tm.calc_t();
// 黄道傾斜角計算
ns::Obliquity o_ob;
eps = o_ob.calc_ob(jcn);
// 結果出力
std::cout << "JST: "
<< ns::gen_time_str(jst) << std::endl;
std::cout << "UTC: "
<< ns::gen_time_str(utc) << std::endl;
std::cout << " JD: "
<< std::fixed << std::setprecision(10) << o_tm.calc_jd()
<< " day" << std::endl;
std::cout << " T: "
<< std::fixed << std::setprecision(10) << jcn
<< " century (= Julian Century Number)" << std::endl;
std::cout << "EPS: "
<< std::fixed << std::setprecision(10) << eps
<< " rad." << std::endl;
std::cout << " "
<< std::fixed << std::setprecision(10) << eps * 180.0 / kPi
<< " deg." << std::endl;
} catch (...) {
std::cerr << "EXCEPTION!" << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
3. ソースコードのビルド（コンパイル＆リンク）
- リポジトリに
Makefileがあるので、それを使用してmakeするだけ。（リビルドする際はmake cleanをしてから） - 上記の
Makefile内では別途個別にインストールしたc++92コマンドを使用しているが、通常はc++であるので注意。
$ make
4. 動作確認
コマンドライン引数に JST（日本標準時）を「年・月・日・時・分・秒・ナノ秒」を最大23桁で指定して実行する。（JST（日本標準時）を指定しない場合は、システム日時を JST とみなす。 JST（日本標準時）を先頭から部分的に指定した場合は、指定していない部分を 0 とみなす）
出力結果の EPS の値が平均黄道傾斜角である（単位: ラジアン、度）。
$ ./calc_obliquity 20210113
JST: 2021-01-13 00:00:00.000
UTC: 2021-01-12 15:00:00.000
JD: 2459227.1250000000 day
T: 0.2103251198 century (= Julian Century Number)
EPS: 0.4090448419 rad.
23.4365430726 deg.
以上。
Comments