Fortran - 太陽・月の視位置計算(海保略算式版)!

更新日時:


Fortran 95 で、海上保安庁・海洋情報部の「コンピュータによる天体の位置計算式」を利用して、太陽や月の視位置等を計算してみました。

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

0. 前提条件

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

1. 海保略算式について

  • コンピュータによる天体の位置計算式」内の PDF ドキュメントや、(当記事冒頭に記載の)当ブログ過去記事を参照のこと。
  • 当然ながら、用意されている係数データファイルの年しか値を計算できない。

2. ソースコードの作成

  • 1ファイルで全てをまかなえるようにしている。(モジュール化はしていない)

File: apparent_sun_moon_jcg.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
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
!*******************************************************************************
! 海上保安庁の天測暦より太陽・月の視位置を計算
! * 視黄経・視黄緯の計算を追加したもの
!
!   Date          Author          Version
!   2018.11.07    mk-mode.com     1.00 新規作成
!   2018.11.10    mk-mode.com     1.01 テキストファイル OPEN/READ 時のエラー処理
!                                      を変更
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
! ---
! 引数: UST(協定世界時)
!       * 書式:YYYYMMDDHHMMSS (年・月・日・時・分・秒)
!       * 無指定なら現在(システム日時)を UTC とみなす
!       * 対応範囲は計算用係数データの用意されている年のみ(2008年〜)
! ---
! * 係数ファイルは、扱いやすくするために1つのファイルにまとめた上、整形している
!   (1行は、年, a, b, 係数)
!*******************************************************************************
!
program jcg_ephemeris
  implicit none

  ! SP: 単精度(4), DP: 倍精度(8)
  integer,      parameter :: SP = kind(1.0)
  integer(SP),  parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
  real(DP),     parameter :: PI         = atan(1.0_DP) * 4.0_DP  ! 円周率
  real(DP),     parameter :: PI_180     = PI / 180.0_DP          ! PI / 180
  character(*), parameter :: F_DELTA_T  = "DELTA_T.txt"          ! ΔT一覧ファイル
  character(*), parameter :: F_SUN_RA   = "SUN_RA.txt"           ! 係数ファイル
  character(*), parameter :: F_SUN_DEC  = "SUN_DEC.txt"          ! 係数ファイル
  character(*), parameter :: F_SUN_DIST = "SUN_DIST.txt"         ! 係数ファイル
  character(*), parameter :: F_MOON_RA  = "MOON_RA.txt"          ! 係数ファイル
  character(*), parameter :: F_MOON_DEC = "MOON_DEC.txt"         ! 係数ファイル
  character(*), parameter :: F_MOON_HP  = "MOON_HP.txt"          ! 係数ファイル
  character(*), parameter :: F_R        = "R.txt"                ! 係数ファイル
  character(*), parameter :: F_EPS      = "EPS.txt"              ! 係数ファイル
  integer(SP),  parameter :: UID_DT     = 10                     ! ΔT一覧ファイル読み込み用
  integer(SP),  parameter :: UID_C      = 11                     ! 係数ファイル読み込み用
  ! 構造型
  type :: t_time   ! 日時
    integer(SP) :: year   = 0
    integer(SP) :: month  = 0
    integer(SP) :: day    = 0
    integer(SP) :: hour   = 0
    integer(SP) :: minute = 0
    integer(SP) :: second = 0
  end type t_time
  type :: t_dt     ! ΔT
    integer(SP) :: year = 0
    integer(SP) :: dt   = 0
  end type t_dt
  type :: t_sun    ! 係数(太陽)
    integer(SP) :: year  = 0
    integer(SP) :: a     = 0
    integer(SP) :: b     = 0
    real(DP)    :: c(18) = 0.0_DP
  end type t_sun
  type :: t_moon   ! 係数(月)
    integer(SP) :: year  = 0
    integer(SP) :: a     = 0
    integer(SP) :: b     = 0
    real(DP)    :: c(30) = 0.0_DP
  end type t_moon
  type :: t_r      ! 係数(R, ε)
    integer(SP) :: year  = 0
    integer(SP) :: a     = 0
    integer(SP) :: b     = 0
    real(DP)    :: c(8)  = 0.0_DP
  end type t_r
  type(t_dt), allocatable :: dts(:)                 ! ΔT一覧
  type(t_dt)   :: tmp(20)                           ! ΔT一覧は最大20年分対応
  integer(SP)  :: len_d                             ! ΔT件数
  type(t_time) :: utc                               ! UTC(協定世界時)
  logical      :: stat                              ! コマンドライン引数取得ステータス
  integer(SP)  :: t, dt                             ! 通日 T, ΔT
  real(DP)     :: f                                 ! 世界時 UT(時・分・秒) の端数
  real(DP)     :: tm, tm_r                          ! 計算用時刻引数(tm_r は R 用)
  type(t_sun)  :: c_sun_ra, c_sun_dec, c_sun_dist   ! 該当区間係数配列
  type(t_moon) :: c_moon_ra, c_moon_dec, c_moon_hp  ! 該当区間係数配列
  type(t_r)    :: c_r, c_eps                        ! 該当区間係数配列
  real(DP)     :: sun_ra, sun_dec, sun_dist         ! 計算結果
  real(DP)     :: moon_ra, moon_dec, moon_hp        ! 計算結果
  real(DP)     :: r, eps                            ! 計算結果
  real(DP)     :: sun_h, moon_h, sun_sd, moon_sd    ! 計算結果
  real(DP)     :: sun_lmd, sun_bet                  ! 計算結果
  real(DP)     :: moon_lmd, moon_bet                ! 計算結果
  real(DP)     :: lmd_s_m                           ! 計算結果

  ! コマンドライン引数(UTC)取得
  call get_arg(utc, stat)
  if (.not. stat) stop

  ! ΔT一覧取得
  ! (+ΔT一覧配列用メモリ確保)
  call get_dt_list(tmp, len_d)
  allocate(dts(len_d))
  dts(1:len_d) = tmp(1:len_d)

  ! 計算対象年がΔT一覧に含まれなければ終了
  if (utc%year < dts(1)%year .or. utc%year > dts(len_d)%year) then
    print '("* Year must be between ", I0, " and ", I0, "!")', &
      & dts(1)%year, dts(len_d)%year
    stop
  end if

  ! 事前計算
  call calc_t(utc, t)               ! 通日 T
  call calc_f(utc, f)               ! 世界時 UT(時・分・秒) の端数
  call get_dt(dts, utc%year, dt)    ! ΔT
  call calc_tm(t, f, dt, tm, tm_r)  ! 計算用時刻引数

  ! 対象区間の係数データ取得
  call get_coeff_sun( F_SUN_RA,   utc%year, tm,   c_sun_ra  )
  call get_coeff_sun( F_SUN_DEC,  utc%year, tm,   c_sun_dec )
  call get_coeff_sun( F_SUN_DIST, utc%year, tm,   c_sun_dist)
  call get_coeff_moon(F_MOON_RA,  utc%year, tm,   c_moon_ra )
  call get_coeff_moon(F_MOON_DEC, utc%year, tm,   c_moon_dec)
  call get_coeff_moon(F_MOON_HP,  utc%year, tm,   c_moon_hp )
  call get_coeff_r(   F_R,        utc%year, tm_r, c_r       )
  call get_coeff_r(   F_EPS,      utc%year, tm,   c_eps     )

  ! 各種計算
  call calc_sun(tm, c_sun_ra,   sun_ra, .true.)
  call calc_sun(tm, c_sun_dec,  sun_dec)
  call calc_sun(tm, c_sun_dist, sun_dist)
  call calc_moon(tm, c_moon_ra,  moon_ra, .true.)
  call calc_moon(tm, c_moon_dec, moon_dec)
  call calc_moon(tm, c_moon_hp,  moon_hp)
  call calc_r(tm_r, c_r,   r,  .true.)
  call calc_r(tm,   c_eps, eps)
  call calc_h(f, sun_ra,  r, sun_h )
  call calc_h(f, moon_ra, r, moon_h)
  call calc_sd_sun(sun_dist, sun_sd)
  call calc_sd_moon(moon_hp, moon_sd)
  call calc_lmd(sun_ra,  sun_dec,  eps, sun_lmd)
  call calc_lmd(moon_ra, moon_dec, eps, moon_lmd)
  call calc_bet(sun_ra,  sun_dec,  eps, sun_bet)
  call calc_bet(moon_ra, moon_dec, eps, moon_bet)
  call calc_lmd_s_m(sun_lmd, moon_lmd, lmd_s_m)

  ! ΔT一覧配列用メモリ解放
  deallocate(dts)

  ! 結果出力
  call display

  stop
contains
  ! コマンドライン引数取得
  ! * YYYYMMDDHHMMSS 形式
  ! * 14桁超入力された場合は、15桁目以降の部分は切り捨てる
  ! * コマンドライン引数がなければ、システム日付を UTC とする
  ! * 日時の整合性チェックは行わない
  !
  ! :param(inout) type(t_date) utc: UTC(協定世界時)
  ! :param(out)   logical     stat: T: 正常, F: 異常
  subroutine get_arg(utc, stat)
    implicit none
    type(t_time), intent(inout) :: utc
    logical,      intent(out)   :: stat
    character(14) :: gc
    integer(SP)   :: dt(8)

    stat = .true.
    if (iargc() == 0) then
      call date_and_time(values=dt)
      utc = t_time(dt(1), dt(2), dt(3), dt(5), dt(6), dt(7))
    else
      call getarg(1, gc)
      if (len(trim(gc)) /= 14) then
        print *, "*** Format must be 'YYYYMMDDHHMMSS'"
        stat = .false.
        return
      end if
      read (gc, '(I4(5I2))') utc
    end if
  end subroutine get_arg

  ! ΔT(UT1(世界時1) - TT(地球時))一覧の取得
  !
  ! :param(out) type(t_dt) dts(50): ΔT一覧
  ! :param(out) integer(4)   len_h: データ件数
  subroutine get_dt_list(dts, len_d)
    implicit none
    type(t_dt),  intent(out) :: dts(20)
    integer(SP), intent(out) :: len_d
    integer(SP)   :: ios, i, year, dt

    ! ΔT一覧 TXT ファイル OPEN
    open (unit   = UID_DT,      &
        & iostat = ios,         &
        & file   = F_DELTA_T,   &
        & action = "read",      &
        & form   = "formatted", &
        & status = "old")
    if (ios /= 0) then
      print '("[ERROR:", I0 ,"] Failed to open file: ", A)', ios, F_DELTA_T
      stop
    end if

    ! ΔT一覧 TXT ファイル READ
    i = 0
    do
      read (UID_DT, fmt=*, iostat = ios) year, dt
      if (ios < 0) then
        exit
      else if (ios > 0) then
        print '("[ERROR:", I0 ,"] Failed to read file: ", A)', ios, F_DELTA_T
      end if
      i = i + 1
      dts(i) = t_dt(year, dt)
    end do
    len_d = i

    ! ΔT一覧 TXT ファイル CLOSE
    close(UID_DT)
  end subroutine get_dt_list

  ! 対象年の ΔT(UT1(世界時1) - TT(地球時))取得
  !
  ! :param(in)  type(t_dt)  dts: ΔT一覧
  ! :param(in)  integer(4) year: 対象の西暦年
  ! :param(out) integer(4)   dt: 対象のΔT
  subroutine get_dt(dts, year, dt)
    implicit none
    type(t_dt),  intent(in)  :: dts(:)
    integer(SP), intent(in)  :: year
    integer(SP), intent(out) :: dt
    integer(SP) :: i

    dt = 0
    do i = 1, size(dts)
      if (dts(i)%year == year) then
        dt = dts(i)%dt
        exit
      end if
    end do
  end subroutine get_dt

  ! 通日 T の計算
  ! * 通日 T は2月0日を第0日とした通算日数で、次式により求める。
  !     T = 30 * P + Q * (S - Y) + P * (1 - Q) + 日
  !   但し、
  !     P = 月 - 1, Q = [(月 + 7) / 10]
  !     Y = [(年 / 4) - [(年 / 4)] + 0.77]
  !     S = [P * 0.55 - 0.33]
  !   で、[] は整数部のみを抜き出すことを意味する。
  !
  ! :param(in)  type(t_time) utc: UTC
  ! :param(out) integer(4)     t: 通日 T
  subroutine calc_t(utc, t)
    implicit none
    type(t_time), intent(in)  :: utc
    integer(SP),  intent(out) :: t
    integer(SP) :: p, q, y, s

    p = utc%month - 1
    q = int((utc%month + 7) / 10)
    y = int(utc%year / 4.0_DP - int(utc%year / 4.0_DP) + 0.77_DP)
    s = int(p * 0.55_DP - 0.33_DP)
    t = 30 * p + q * (s - y) + p * (1 - q) + utc%day
  end subroutine calc_t

  ! 世界時(時・分・秒)の端数計算
  ! * 次式により求める。
  !     F = 時 / 24 + 分 / 1440 + 秒 / 86400
  !
  ! :param(in)  type(t_time) utc: UTC
  ! :param(out) integer(4)     f: 世界時の端数
  subroutine calc_f(utc, f)
    implicit none
    type(t_time), intent(in)  :: utc
    real(DP),     intent(out) :: f

    f = utc%hour   /    24.0_DP &
    & + utc%minute /  1440.0_DP &
    & + utc%second / 86400.0_DP
  end subroutine calc_f

  ! 計算用時刻引数 tm の計算
  ! * 次式により求める。
  !   (R 計算用は tm_r, その他は tm)
  !     tm   = T + F + ΔT / 86400
  !     tm_r = T + F
  !
  ! :param(in)  integer(4)  t: 通日 T
  ! :param(in)  real(8)     f: 世界時の端数
  ! :param(in)  integer(4) dt: ΔT
  ! :param(out) real(4)    tm: 時刻引数
  ! :param(out) real(4)    tm: 時刻引数(R 用)
  subroutine calc_tm(t, f, dt, tm, tm_r)
    implicit none
    integer(SP), intent(in)  :: t, dt
    real(DP),    intent(in)  :: f
    real(DP),    intent(out) :: tm, tm_r

    tm_r = t + f
    tm   = tm_r + dt / 86400.0_DP
  end subroutine calc_tm

  ! 対象年の係数データ取得
  ! * SUN(R.A., Dec., Dist.) 用
  !
  ! :param(in)  character(*) f_name: 係数ファイルの名称
  ! :param(in)  integer(4)     year: 対象の西暦年
  ! :param(in)  real(8)          tm: 時刻引数
  ! :param(out) type(t_sun)    item: 係数一覧
  subroutine get_coeff_sun(f_name, year, tm, item)
    implicit none
    character(*), intent(in)  :: f_name
    integer(SP),  intent(in)  :: year
    real(DP),     intent(in)  :: tm
    type(t_sun),  intent(out) :: item
    integer(SP) :: ios, y, a, b
    real(DP)    :: c(18)

    ! 係数 TXT ファイル OPEN
    open (unit   = UID_C,        &
        & iostat = ios,          &
        & file   = trim(f_name), &
        & action = "read",       &
        & form   = "formatted",  &
        & status = "old")
    if (ios /= 0) then
      print '("[ERROR:", I0 ,"] Failed to open file: ", A)', ios, f_name
      stop
    end if

    ! 係数 TXT ファイル READ
    do
      read (UID_C, fmt=*, iostat = ios) y, a, b, c
      if (ios < 0) then
        exit
      else if (ios > 0) then
        print '("[ERROR:", I0 ,"] Failed to read file: ", A)', ios, f_name
      end if
      if (y /= year) cycle
      if (int(tm) < a .or. b < int(tm)) cycle
      item = t_sun(y, a, b, c)
      exit
    end do

    ! 係数 TXT ファイル CLOSE
    close(UID_C)
  end subroutine get_coeff_sun

  ! 対象年の係数データ取得
  ! * MOON(R.A., Dec., H.P..) 用
  !
  ! :param(in)  character(*) f_name: 係数ファイルの名称
  ! :param(in)  integer(4)     year: 対象の西暦年
  ! :param(in)  real(8)         tm: 時刻引数
  ! :param(out) type(t_moon)  item: 係数一覧
  subroutine get_coeff_moon(f_name, year, tm, item)
    implicit none
    character(*), intent(in)  :: f_name
    integer(SP),  intent(in)  :: year
    real(DP),     intent(in)  :: tm
    type(t_moon), intent(out) :: item
    integer(SP) :: ios, y, a, b
    real(DP)    :: c(30)

    ! 係数 TXT ファイル OPEN
    open (unit   = UID_C,        &
        & iostat = ios,          &
        & file   = trim(f_name), &
        & action = "read",       &
        & form   = "formatted",  &
        & status = "old")
    if (ios /= 0) then
      print '("[ERROR:", I0 ,"] Failed to open file: ", A)', ios, f_name
      stop
    end if

    ! 係数 TXT ファイル READ
    do
      read (UID_C, fmt=*, iostat = ios) y, a, b, c
      if (ios < 0) then
        exit
      else if (ios > 0) then
        print '("[ERROR:", I0 ,"] Failed to read file: ", A)', ios, f_name
      end if
      if (y /= year) cycle
      if (int(tm) < a .or. b < int(tm)) cycle
      item = t_moon(y, a, b, c)
      exit
    end do

    ! 係数 TXT ファイル CLOSE
    close(UID_C)
  end subroutine get_coeff_moon

  ! 対象年の係数データ取得
  ! * R, EPS 用
  !
  ! :param(in)  character(*) f_name: 係数ファイルの名称
  ! :param(in)  integer(4)     year: 対象の西暦年
  ! :param(in)  real(8)          tm: 時刻引数
  ! :param(out) type(t_r)      item: 係数一覧
  subroutine get_coeff_r(f_name, year, tm, item)
    implicit none
    character(*), intent(in)  :: f_name
    integer(SP),  intent(in)  :: year
    real(DP),     intent(in)  :: tm
    type(t_r),    intent(out) :: item
    integer(SP) :: ios, y, a, b
    real(DP)    :: c(8)

    ! 係数 TXT ファイル OPEN
    open (unit   = UID_C,        &
        & iostat = ios,          &
        & file   = trim(f_name), &
        & action = "read",       &
        & form   = "formatted",  &
        & status = "old")
    if (ios /= 0) then
      print '("[ERROR:", I0 ,"] Failed to open file: ", A)', ios, f_name
      stop
    end if

    ! 係数 TXT ファイル READ
    do
      read (UID_C, fmt=*, iostat = ios) y, a, b, c
      if (ios < 0) then
        exit
      else if (ios > 0) then
        print '("[ERROR:", I0 ,"] Failed to read file: ", A)', ios, f_name
      end if
      if (y /= year) cycle
      if (int(tm) < a .or. b < int(tm)) cycle
      item = t_r(y, a, b, c)
      exit
    end do

    ! 係数 TXT ファイル CLOSE
    close(UID_C)
  end subroutine get_coeff_r

  ! SUN(R.A., Dec., Dist.) の計算
  ! * θ, 係数配列から次式により所要値を計算する。
  !     f(t) = C_0 + C_1 * cos(θ) + C_2 * cos(2θ) + ... + C_N * cos(Nθ)
  !
  ! :param(in)  real(8)       tm: 時刻引数
  ! :param(in)  type(t_sun) coef: 係数
  ! :param(out) real(8)      val: 計算結果
  ! :param(in)  logical, optional is_ra: R.A. => T, その他 => F
  subroutine calc_sun(tm, coef, val, is_ra)
    implicit none
    real(DP),    intent(in)  :: tm
    type(t_sun), intent(in)  :: coef
    logical,     intent(in), optional :: is_ra
    real(DP),    intent(out) :: val
    integer(SP) :: i
    real(DP)    :: th

    val = 0.0
    do i = 1, size(coef%c)
      th = theta(real(coef%a, DP), real(coef%b, DP), tm) * (i - 1) * PI_180
      val = val + coef%c(i) * cos(th)
    end do
    if (.not. present(is_ra)) return
    do while (val >= 24.0_DP)
      val = val - 24.0_DP
    end do
    do while (val <= 0.0_DP)
      val = val + 24.0_DP
    end do
  end subroutine calc_sun

  ! MOON(R.A., Dec., H.P.) の計算
  ! * θ, 係数配列から次式により所要値を計算する。
  !     f(t) = C_0 + C_1 * cos(θ) + C_2 * cos(2θ) + ... + C_N * cos(Nθ)
  !
  ! :param(in)  real(8)        tm: 時刻引数
  ! :param(in)  type(t_moon) coef: 係数
  ! :param(out) real(8)       val: 計算結果
  ! :param(in)  logical, optional is_ra: R.A. => T, その他 => F
  subroutine calc_moon(tm, coef, val, is_ra)
    implicit none
    real(DP),     intent(in)  :: tm
    type(t_moon), intent(in)  :: coef
    logical,      intent(in), optional :: is_ra
    real(DP),     intent(out) :: val
    integer(SP) :: i
    real(DP)    :: th

    val = 0.0
    do i = 1, size(coef%c)
      th = theta(real(coef%a, DP), real(coef%b, DP), tm) * (i - 1) * PI_180
      val = val + coef%c(i) * cos(th)
    end do
    if (.not. present(is_ra)) return
    do while (val >= 24.0_DP)
      val = val - 24.0_DP
    end do
    do while (val <= 0.0_DP)
      val = val + 24.0_DP
    end do
  end subroutine calc_moon

  ! R, EPS の計算
  ! * θ, 係数配列から次式により所要値を計算する。
  !     f(t) = C_0 + C_1 * cos(θ) + C_2 * cos(2θ) + ... + C_N * cos(Nθ)
  !
  ! :param(in)  real(8)     tm: 時刻引数
  ! :param(in)  type(t_r) coef: 係数
  ! :param(out) real(8)    val: 計算結果
  ! :param(in)  logical, optional is_r: R => T, EPS => F
  subroutine calc_r(tm, coef, val, is_r)
    implicit none
    real(DP),  intent(in)  :: tm
    type(t_r), intent(in)  :: coef
    logical,   intent(in), optional :: is_r
    real(DP),  intent(out) :: val
    integer(SP) :: i
    real(DP)    :: th

    val = 0.0
    do i = 1, size(coef%c)
      th = theta(real(coef%a, DP), real(coef%b, DP), tm) * (i - 1) * PI_180
      val = val + coef%c(i) * cos(th)
    end do
    if (.not. present(is_r)) return
    do while (val >= 24.0_DP)
      val = val - 24.0_DP
    end do
    do while (val <= 0.0_DP)
      val = val + 24.0_DP
    end do
  end subroutine calc_r

  ! θ の計算
  ! * θ を次式により計算する。(単位: °)
  !     θ = cos^(-1)((2 * t - (a + b)) / (b - a))
  !   但し、0°<= θ <= 180°
  !
  ! :param(in)  real(8)   t_a: 区間 a
  ! :param(in)  real(8)   t_b: 区間 b
  ! :param(in)  real(8)     t: 時刻引数
  ! :return     real(8) theta: θ
  real(DP) function theta(t_a, t_b, t)
    implicit none
    real(DP), intent(in) :: t_a, t_b, t
    real(DP) :: a, b

    a = t_a
    b = t_b
    ! 年末のΔT秒分も計算可能とするための応急処置
    if (t_b < t) b = t
    theta = (2 * t - (a + b)) / (b - a)
    theta = acos(theta) * 180 / PI
  end function theta

  ! グリニジ時角の計算
  ! * 次式によりグリニジ時角を計算する。
  !     h = E + UT
  !   (但し、E = R - R.A.)
  !
  ! :param(in)  real(8)  f: 世界時 UT(時・分・秒) の端数
  ! :param(in)  real(8) ra: R.A.
  ! :param(in)  real(8)  r: R
  ! :param(out) real(8)  h: グリニジ時角
  subroutine calc_h(f, ra,  r, h)
    implicit none
    real(DP), intent(in)  :: f, ra, r
    real(DP), intent(out) :: h

    h = r - ra + f * 24.0_DP
  end subroutine calc_h

  ! 視半径(太陽)の計算
  ! * 次式により視半径を計算する。
  !     S.D.= 16.02 ′/ Dist.
  !
  ! :param(in)  real(8) dist: SUN(Dist.)
  ! :param(out) real(8)   sd: SUN(S.D.)
  subroutine calc_sd_sun(dist, sd)
    implicit none
    real(DP), intent(in)  :: dist
    real(DP), intent(out) :: sd

    sd = 16.02_DP / dist
  end subroutine calc_sd_sun

  ! 視半径(月)の計算
  ! * 次式により視半径を計算する。
  !     S.D.= sin^(-1) (0.2725 * sin(H.P.))
  !
  ! :param(in)  real(8) hp: MOON(H.P.)
  ! :param(out) real(8) sd: MOON(S.D.)
  subroutine calc_sd_moon(hp, sd)
    implicit none
    real(DP), intent(in)  :: hp
    real(DP), intent(out) :: sd

    sd = asin(0.2725_DP * sin(hp * PI_180)) * 60.0_DP * 180.0_DP / PI
  end subroutine calc_sd_moon

  ! 視黄経の計算
  ! * 次式により視黄経を計算する
  !     λ = arctan(sin δ sin ε + cos δ sin α cos ε / cos δ cos α)
  !   (α : 視赤経、δ : 視赤緯、 ε : 黄道傾斜角)
  !
  ! :param(in)  real(8)  ra: R.A.(視赤経)
  ! :param(in)  real(8) dec: Dec.(視赤緯)
  ! :param(in)  real(8) eps: ε(黄道傾角)
  ! :param(out) real(8) lmd: λ(視黄経)
  subroutine calc_lmd(ra, dec, eps, lmd)
    implicit none
    real(DP), intent(in)  :: ra, dec, eps
    real(DP), intent(out) :: lmd
    real(DP) :: a, d, e, l_a, l_b

    a = ra * 15.0_DP * PI_180
    d = dec * PI_180
    e = eps * PI_180
    l_a = sin(d) * sin(e) + cos(d) * sin(a) * cos(e)
    l_b = cos(d) * cos(a)
    lmd = atan2(l_a, l_b) * 180.0_DP / PI
    if (lmd < 0.0_DP) lmd = lmd + 360.0_DP
  end subroutine calc_lmd

  ! 視黄緯の計算
  ! * 次式により視黄経を計算する
  !     β  = arcsin(sin δ cos ε − cos δ sin α sin ε)
  !   (α : 視赤経、δ : 視赤緯、 ε : 黄道傾斜角)
  !
  ! :param(in)  real(8)  ra: R.A.(視赤経)
  ! :param(in)  real(8) dec: Dec.(視赤緯)
  ! :param(in)  real(8) eps: ε(黄道傾角)
  ! :param(out) real(8) bet: β(視黄緯)
  subroutine calc_bet(ra, dec, eps, bet)
    implicit none
    real(DP), intent(in)  :: ra, dec, eps
    real(DP), intent(out) :: bet
    real(DP) :: a, d, e

    a = ra * 15.0_DP * PI_180
    d = dec * PI_180
    e = eps * PI_180
    bet = sin(d) * cos(e) - cos(d) * sin(a) * sin(e)
    bet = asin(bet) * 180.0_DP / PI
  end subroutine calc_bet

  ! 視黄経差(太陽 - 月)の計算
  ! * SUN_LAMBDA - MOON_LAMBDA
  !
  ! :param(in)  real(8)  sun_lmd: 視黄経(太陽)
  ! :param(in)  real(8) moon_lmd: 視黄経(月)
  ! :param(out) real(8)  lmd_s_m: 視黄経差(太陽 - 月)
  subroutine calc_lmd_s_m(lmd_s, lmd_m, lmd_s_m)
    implicit none
    real(DP), intent(in)  :: lmd_s, lmd_m
    real(DP), intent(out) :: lmd_s_m

    lmd_s_m = lmd_s - lmd_m
  end subroutine calc_lmd_s_m

  ! 結果出力
  ! * param なし
  subroutine display
    implicit none

    print '("[ ", I4, "-", I0.2, "-", I0.2, " ", &
      & I0.2, ":", I0.2, ":", I0.2, " UTC ]")', utc
    print '("SUN    R.A. = ", F12.8, " h  (= ", A, ")")', &
        & sun_ra,   hour2hms(sun_ra)
    print '("SUN    DEC. = ", F12.8, " ° (= ", A, ")")', &
        & sun_dec,  deg2dms(sun_dec)
    print '("SUN   DIST. = ", F12.8, " AU")', sun_dist
    print '("SUN     hG. = ", F12.8, " h  (= ", A, ")")', &
        & sun_h,    hour2hms(sun_h)
    print '("SUN    S.D. = ", F12.8, " ′ (= ", A, ")")', &
        & sun_sd,   deg2dms(sun_sd / 60.0_DP)
    print '("MOON   R.A. = ", F12.8, " h  (= ", A, ")")', &
        & moon_ra,  hour2hms(moon_ra)
    print '("MOON   DEC. = ", F12.8, " ° (= ", A, ")")', &
        & moon_dec, deg2dms(moon_dec)
    print '("MOON   H.P. = ", F12.8, " ° (= ", A, ")")', &
        & moon_hp,  deg2dms(moon_hp)
    print '("MOON    hG. = ", F12.8, " h  (= ", A, ")")', &
        & moon_h,   hour2hms(moon_h)
    print '("MOON   S.D. = ", F12.8, " ′ (= ", A, ")")', &
        & moon_sd,  deg2dms(moon_sd / 60.0_DP)
    print '("         R  = ", F12.8, " h  (= ", A, ")")', &
        & r,        hour2hms(r)
    print '("       EPS. = ", F12.8, " ° (= ", A, ")")', &
        & eps,      deg2dms(eps)
    print '(A)', "---"
    print '("SUN  LAMBDA = ", F12.8, " ° (= ", A, ")")', &
        & sun_lmd,  deg2dms(sun_lmd)
    print '("SUN    BETA = ", F12.8, " ° (= ", A, ")")', &
        & sun_bet,  deg2dms(sun_bet)
    print '("MOON LAMBDA = ", F12.8, " ° (= ", A, ")")', &
        & moon_lmd, deg2dms(moon_lmd)
    print '("MOON   BETA = ", F12.8, " ° (= ", A, ")")', &
        & moon_bet, deg2dms(moon_bet)
    print '("DIFF LAMBDA = ", F12.8, " ° (= ", A, ")")', &
        & lmd_s_m,  deg2dms(lmd_s_m)
  end subroutine display

  ! 99.999h -> 99h99m99s 変換
  !
  ! :param(in) real(8)      hour: 時間
  ! :return    character(19) hms: "999 h 99 m 99.999 s" 形式
  function hour2hms(hour) result(hms)
    implicit none
    real(DP), intent(in) :: hour
    character(19) :: hms
    integer(SP)   :: h, m
    real(DP)      :: h_r, m_r, s

    h   = int(hour)
    h_r = hour - real(h, DP)
    m   = int(h_r * 60.0_DP)
    m_r = h_r * 60.0_DP - real(m, DP)
    s   = m_r * 60.0_DP
    if (m < 0) m = isign(m, 1)
    if (s < 0.0_DP) s = dsign(s, 1.0_DP)
    write (hms, '(I3, " h ", I2, " m ", F6.3, " s")') h, m, s
  end function hour2hms

  ! 99.999° -> 99°99′99″ 変換
  !
  ! :param(in) real(8)       deg: 角度
  ! :return    character(24) dms: "999 ° 99 ′ 99.999 ″" 形式
  function deg2dms(deg) result(dms)
    implicit none
    real(DP), intent(in) :: deg
    character(24) :: dms
    character(1)  :: pm
    integer(SP)   :: d, m
    real(DP)      :: d_0, d_r, m_r, s

    d_0 = deg
    pm = " "
    if (d_0 < 0.0_DP) pm = "-"
    if (d_0 < 0.0_DP) d_0 = dsign(d_0, 1.0_DP)
    d = int(d_0)
    d_r = d_0 - real(d, DP)
    m   = int(d_r * 60.0_DP)
    m_r = d_r * 60.0_DP - real(m, DP)
    s   = m_r * 60.0_DP
    if (m < 0) m = isign(m, 1)
    if (s < 0.0_DP) s = dsign(s, 1.0_DP)
    write (dms, '(A1, I3, " °", I3, " ′", F6.3, " ″")') pm, d, m, s
  end function deg2dms
end program jcg_ephemeris

定数ファイルや、定数ファイル生成用の Ruby スクリプトとともに GitHub にアップしている。

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

(1ファイルなので Makefile は使用しない)

$ make -o apparent_sun_moon_jcg apparent_sun_moon_jcg.f95

4. 準備

まず、ΔT(地球時(TT)と世界時1(UT1))の値を記述したファイル DELTA_T.txt が存在することを確認し、不足していれば、値を追加する。(隔年の値は、各年の解説ファイル(na99-rei.pdf)に記載されている)

次に、海保サイトからダウンロードした各年の係数ファイル(na99-data.tx という名称のファイル)を text ディレクトリに配置し、別途作成の Ruby スクリプトで整形し直す。(このスクリプトは、今回必要な太陽・月・R,黄道傾角(ε)にしか対応していない)

$ ./gen_coeffs_f95.rb

同じディレクトリ内に以下の8個のテキストファイルが生成されるはず。

  • SUN_RA.txt
  • SUN_DEC.txt
  • SUN_DIST.txt
  • MOON_RA.txt
  • MOON_DEC.txt
  • MOON_HP.txt
  • R.txt
  • EPS.txt

5. 実行方法

$ ./apparent_sun_moon_jcg [YYYYMMDDHHMMSS]
  • UTC(協定世界時)は「年・月・日・時・分・秒」を14桁で指定する。
  • UTC(協定世界時)を指定しない場合は、システム日時を UTC とみなす。

以下、実行例。

$ ./apparent_sun_moon_jcg 20190121000000
[ 2019-01-21 00:00:00 UTC ]
SUN    R.A. =  20.18948049 h  (=  20 h 11 m 22.130 s)
SUN    DEC. = -20.01113167 ° (= - 20 °  0 ′40.074 ″)
SUN   DIST. =   0.98398646 AU
SUN     hG. = -12.18482730 h  (= -12 h 11 m  5.378 s)
SUN    S.D. =  16.28071178 ′ (=    0 ° 16 ′16.843 ″)
MOON   R.A. =   7.97458715 h  (=   7 h 58 m 28.514 s)
MOON   DEC. =  20.71934023 ° (=   20 ° 43 ′ 9.625 ″)
MOON   H.P. =   1.02075727 ° (=    1 °  1 ′14.726 ″)
MOON    hG. =   0.03006605 h  (=   0 h  1 m 48.238 s)
MOON   S.D. =  16.68856406 ′ (=    0 ° 16 ′41.314 ″)
         R  =   8.00465319 h  (=   8 h  0 m 16.752 s)
       EPS. =  23.43561057 ° (=   23 ° 26 ′ 8.198 ″)
---
SUN  LAMBDA = 300.63618339 ° (=  300 ° 38 ′10.260 ″)
SUN    BETA =   0.00013479 ° (=    0 °  0 ′ 0.485 ″)
MOON LAMBDA = 117.53324938 ° (=  117 ° 31 ′59.698 ″)
MOON   BETA =   0.06976027 ° (=    0 °  4 ′11.137 ″)
DIFF LAMBDA = 183.10293401 ° (=  183 °  6 ′10.562 ″)

以上。

 

Sponsored Link

 

コメントする