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
| !*******************************************************************************
! 章動の計算
! : IAU2000A 章動理論(MHB2000, IERS2003)による
! 黄経における章動(Δψ), 黄道傾斜における章動(Δε) の計算
!
! * IAU SOFA(International Astronomical Union, Standards of Fundamental Astronomy)
! の提供する C ソースコード "nut00a.c" で実装されているアルゴリズムを使用する。
! * 参考サイト
! - [SOFA Library Issue 2012-03-01 for ANSI C: Complete List](http://www.iausofa.org/2012_0301_C/sofa/)
! - [USNO Circular 179](http://aa.usno.navy.mil/publications/docs/Circular_179.php)
! - [IERS Conventions Center](http://62.161.69.131/iers/conv2003/conv2003_c5.html)
!
! Date Author Version
! 2018.10.18 mk-mode.com 1.00 新規作成
! 2018.11.09 mk-mode.com 1.01 時刻の取扱変更(マイクロ秒 => ミリ秒)
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
! ---
! 引数 : 日時(TT（地球時）)
! * 書式：YYYYMMDD[HHMMSS[MMM]]
! (最後の MMM はミリ秒)
! * 無指定なら現在(システム日時)を地球時とみなす。
!*******************************************************************************
!
program nutation_model
use const
use nutation
use time
implicit none
type(t_time) :: tt
real(DP) :: jd, jc
real(DP) :: dpsi_ls, dpsi_pl, dpsi, dpsi_d, dpsi_s
real(DP) :: deps_ls, deps_pl, deps, deps_d, deps_s
! コマンドライン引数(TT)取得
call get_arg(tt)
if (tt%year == 0) stop
! 章動計算
call gc2jd(tt, jd)
call jd2jc(jd, jc)
call calc_lunisolar(jc, dpsi_ls, deps_ls)
call calc_planetary(jc, dpsi_pl, deps_pl)
dpsi = dpsi_ls + dpsi_pl
deps = deps_ls + deps_pl
dpsi_d = dpsi * R2D
deps_d = deps * R2D
dpsi_s = dpsi_d * D2S
deps_s = deps_d * D2S
! 結果出力
print '(" [", A, " TT]")', date_fmt(tt)
print '(" DeltaPsi = ", E22.14e2, " rad")', dpsi
print '(" = ", E22.14e2, " °")', dpsi_d
print '(" = ", E22.14e2, " ″")', dpsi_s
print '(" DeltaEps = ", E22.14e2, " rad")', deps
print '(" = ", E22.14e2, " °")', deps_d
print '(" = ", E22.14e2, " ″")', deps_s
stop
contains
! コマンドライン引数取得
! * YYYYMMDD[HHMMSS[MMM]] 形式
! * 17桁超入力された場合は、21桁目以降の部分は切り捨てる
! * コマンドライン引数がなければ、システム日付を TT とする
! * 日時の整合性チェックは行わない
!
! :param(inout) type(t_time) tt
subroutine get_arg(tt)
implicit none
type(t_time), intent(inout) :: tt
character(20) :: gc
integer(SP) :: dt(8), len_gc
if (iargc() == 0) then
call date_and_time(values=dt)
tt = t_time(dt(1), dt(2), dt(3), dt(5), dt(6), dt(7), dt(8))
else
call getarg(1, gc)
len_gc = len(trim(gc))
if (len_gc /= 8 .and. len_gc /= 14 .and. len_gc /= 17) then
print *, "Format: YYYYMMDD[HHMMSS[MMM]]"
return
end if
read (gc, FMT_DT_0) tt
end if
end subroutine get_arg
end program nutation_model
|