Fortran - 日・月の出・南中・入時刻の計算!

更新日時:


Fortran 95 で、日・月の出・南中・入時刻を計算してみました。(出・入はその時の方位角、南中はその時の高度も)

過去に Ruby で行ったことはありましたが。

0. 前提条件

  • LMDE 3 (Linux Mint Debian Edition 3; 64bit) での作業を想定。
  • GCC 6.3.0 (GFortran 6.3.0) でのコンパイルを想定。

1. 計算方法について

計算アルゴリズムは「日の出・日の入りの計算―天体の出没時刻の求め方」によるもの。

2. ソースコードの作成

以下は実行部分。

File: sun_moon.f95

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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
!*******************************************************************************
! 日・月の出・南中・入時刻、日・月の出・入方位角、日・月の南中高度の計算
!
!   date          name            version
!   2018.11.11    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
! ---
! 引数 : 99999999 [+-]999.99 [+-]999.99 [+]9999.99
!        第1: 日付 [必須]
!             計算対象の日付(グレゴリオ暦)を半角8桁数字で指定
!        第2: 緯度 [必須]
!             緯度を 度 で指定(度・分・秒は度に換算して指定すること)
!             (北緯はプラス、南緯はマイナス。桁数は特に制限なし)
!        第3: 経度 [必須]
!             経度を 度 で指定(度・分・秒は度に換算して指定すること)
!             (東経はプラス、西経はマイナス。桁数は特に制限なし)
!        第4: 標高 [必須]
!             標高をメートルで指定(マイナス値は指定不可)
!             (桁数は特に制限なし)
!*******************************************************************************
!
program sun_moon
  use const, only : SP, DP, SEC_DAY
  use time
  use delta_t
  use calc
  implicit none
  type(t_time) :: jst, utc
  type(t_time) :: tm_sr, tm_ss, tm_sm, tm_mr, tm_ms, tm_mm
  integer(SP)  :: utc_tai
  real(DP)     :: lat, lon, ht, dut1, dt, dt_d, dip, day_p
  real(DP)     :: ang_sr, ang_ss, ang_mr, ang_ms, ht_sm, ht_mm

  ! コマンドライン引数取得
  call get_arg(jst, lat, lon, ht)

  ! 初期処理
  call jst2utc(jst, utc)               ! JST -> UTC
  call utc2utc_tai(utc, utc_tai)       ! UTC -> UTC - TAI
  call utc2dut1(utc, dut1)             ! UTC -> DUT1
  call utc2dt(utc, utc_tai, dut1, dt)  ! UTC -> delta T
  dt_d = dt / SEC_DAY                  ! ΔTを日換算
  dip = 0.0353333_DP * sqrt(ht)        ! 地平線伏角
  day_p = day_progress(jst)            ! 2000年1月1日力学時正午からの経過日数(日)

  ! 各種計算
  call calc_sun_sr( 0, lon, lat, day_p, dt_d, dip, tm_sr, ang_sr)  ! 日の出
  call calc_sun_sr( 1, lon, lat, day_p, dt_d, dip, tm_ss, ang_ss)  ! 日の入
  call calc_sun_m(     lon, lat, day_p, dt_d, dip, tm_sm, ht_sm )  ! 日の南中
  call calc_moon_sr(0, lon, lat, day_p, dt_d, dip, tm_mr, ang_mr)  ! 月の出
  call calc_moon_sr(1, lon, lat, day_p, dt_d, dip, tm_ms, ang_ms)  ! 月の入
  call calc_moon_m(    lon, lat, day_p, dt_d, dip, tm_mm, ht_mm )  ! 月の南中

  ! 結果出力
  call display

  stop
contains
  ! コマンドライン引数取得
  !
  ! :param(out) type(t_time) jst: 日付
  ! :param(out) real(8)      lat: 緯度
  ! :param(out) real(8)      lon: 経度
  ! :param(out) real(8)       ht: 標高
  subroutine get_arg(jst, lat, lon ,ht)
    implicit none
    type(t_time), intent(out) :: jst
    real(DP),     intent(out) :: lat, lon ,ht
    character(32) :: buf
    integer(SP)   :: dt(8), ios, y, m, d

    if (iargc() /= 4) then
      print *, "Arguments: YYYYMMDD LATITUDE LONGITUDE HEIGHT"
      stop
    end if
    call getarg(1, buf)
    read (buf, '(I4I2I2)', iostat=ios) y, m, d
    if (ios /= 0) then
      print *, "DATE is invalid!"
      stop
    end if
    jst = t_time(y, m, d, 0, 0, 0)
    call getarg(2, buf)
    read (buf, *, iostat=ios) lat
    if (ios /= 0) then
      print *, "LATITUDE is invalid!"
      stop
    end if
    call getarg(3, buf)
    read (buf, *, iostat=ios) lon
    if (ios /= 0) then
      print *, "LONGITUDE is invalid!"
      stop
    end if
    call getarg(4, buf)
    read (buf, *, iostat=ios) ht
    if (ios /= 0) then
      print *, "HEIGHT is invalid!"
      stop
    end if
    if (ht < 0.0_DP) then
      print *, "HEIGHT is invalid!"
      stop
    end if
  end subroutine get_arg

  ! 結果出力
  subroutine display
    implicit none
    character(1) :: s_lat, s_lon

    s_lat = "N"
    s_lon = "E"
    if (lat < 0.0_DP) then
      s_lat = "S"
      lat = dsign(lat, 1.0_DP)
    end if
    if (lon < 0.0_DP) then
      s_lon = "W"
      lon = dsign(lon, 1.0_DP)
    end if
    print '("[", I4, "-", I0.2, "-", I0.2, "JST ", &
        & F0.4, A, X, F0.4, A, X, F0.2, "m]")', &
        & jst%year, jst%month, jst%day, lat, s_lat, lon, s_lon, ht
    print '("日の出 ", I0.2, ":", I0.2, ":", I0.2, &
        & " (方位角 ", F6.2, ")")', &
        & tm_sr%hour, tm_sr%minute, tm_sr%second, ang_sr
    print '("日南中 ", I0.2, ":", I0.2, ":", I0.2, &
        & " ( 高度 ", F6.2, ")")', &
        & tm_sm%hour, tm_sm%minute, tm_sm%second, ht_sm
    print '("日の入 ", I0.2, ":", I0.2, ":", I0.2, &
        & " (方位角 ", F6.2, ")")', &
        & tm_ss%hour, tm_ss%minute, tm_ss%second, ang_ss
    if (ang_mr < 0.0_DP) then
      print '(A)', "月の出 --:--:-- (方位角 ---.--)"
    else
      print '("月の出 ", I0.2, ":", I0.2, ":", I0.2, &
          & " (方位角 ", F6.2, ")")', &
          & tm_mr%hour, tm_mr%minute, tm_mr%second, ang_mr
    end if
    if (ht_mm < 0.0_DP) then
      print '(A)', "月南中 --:--:-- ( 高度 ---.--)"
    else
      print '("月南中 ", I0.2, ":", I0.2, ":", I0.2, &
          & " ( 高度 ", F6.2, ")")', &
          & tm_mm%hour, tm_mm%minute, tm_mm%second, ht_mm
    end if
    if (ang_ms < 0.0_DP) then
      print '(A)', "月の入 --:--:-- (方位角 ---.--)"
    else
      print '("月の入 ", I0.2, ":", I0.2, ":", I0.2, &
          & " (方位角 ", F6.2, ")")', &
          & tm_ms%hour, tm_ms%minute, tm_ms%second, ang_ms
    end if
  end subroutine display
end program sun_moon

モジュールや定数ファイルを含めた全てのファイルは GitHub にアップしている。(当然ながら、モジュール部分の方が計算の本質となっている)

3. ソースコードのビルド

全てのファイルを用意してから、以下を実行。

$ make
  • やり直す場合は、 make clean してから。

4. 準備

  • うるう年ファイル LEAP_SEC.txt, DUT1 ファイル DUT1.txt は適宜最新のものに更新すること。

5. 実行方法

日付、緯度、経度、標高を 99999999 [+-]999.99 [+-]999.99 [+]9999.99 形式指定して実行。

$ .sun_moon 99999999 [+-]999.99 [+-]999.99 [+]9999.99
  • 第1: 日付 [必須]
     計算対象の日付(グレゴリオ暦)を半角8桁数字で指定
  • 第2: 緯度 [必須]
     緯度を 度 で指定(度・分・秒は度に換算して指定すること)
     (北緯はプラス、南緯はマイナス。桁数は特に制限なし)
  • 第3: 経度 [必須]
     経度を 度 で指定(度・分・秒は度に換算して指定すること)
     (東経はプラス、西経はマイナス。桁数は特に制限なし)
  • 第4: 標高 [必須]
     標高をメートルで指定(マイナス値は指定不可)
     (桁数は特に制限なし)

以下、実行例。(2019年1月24日、松江市役所の緯度・経度、標高0m で計算)

$ ./sun_moon 20190124 35.4681 133.0486 0
[2019-01-24JST 35.4681N 133.0486E .00m]
日の出 07:13:01 (方位角 113.34)
日南中 12:19:43 ( 高度  35.26)
日の入 17:26:45 (方位角 246.79)
月の出 21:04:39 (方位角  80.76)
月南中 02:44:31 ( 高度  66.04)
月の入 09:26:40 (方位角 282.21)

以上。

 

Sponsored Link

 

コメントする