JPL 天文暦データのバイナリ化(Fortran 95 Ver.)!

更新日時:


以前、複数存在する DE430 のテキスト形式データファイルを1つのバイナリ形式のファイルにコンバートする方法についての記録しました。

但し、 FORTRAN 77 のソースコードでした。

今回、 Fortran 95 に移植しました。以下、その記録です。

0. 前提条件

  • LMDE 3 (Linux Mint Debian Edition 3; 64bit) での作業を想定。
  • GNU Fortran (GCC) 6.3.0 でのコンパイル作業を想定。(GCC でなくてもよいだろう)
  • 今回移植するのは、テキストファイルをバイナリ化する asc2eph.f のみ。

1. 使用するデータファイル

ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii/de430/ 内の以下のファイル。

  • ヘッダファイル: header.430_572
  • テキスト形式データファイル: ascp1550.430ascp2550.430
  • 検証用ファイル: testpo.430

2. 元の FORTRAN 77 ソースコード

ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran/ 内の以下のファイル。

  • コンバートプログラム: asc2eph.f (<= 今回の移植対象)
  • データ検証プログラム: testeph.f

3. コンバートプログラムの作成

元の FORTRAN 77 製コンバートプログラム asc2eph.f を文法等に注意しながら Fortran 95 用に書き換える。

File: asc2eph.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
!***********************************************************************
! ASC2EPH creates a binary format JPL Planetary Ephemeris file from
! one or more ascii text files.
!
! date          name            version
! 2018.08.22    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!***********************************************************************
!
program main
  implicit none
  ! Unix プラットフォームでなければ nrecl = 1
  integer, parameter :: nrecl  = 4
  integer, parameter :: oldmax = 400
  integer, parameter :: nmax   = 1000
  character(len= 6) :: cnam(nmax)
  character(len= 6) :: ttl(14, 3)
  character(len=12) :: header
  real(8) :: au, cval(nmax), db(3000), db2z
  real(8) :: emrat, ss(3), t1, t2
  integer :: i, j, irecsz, in, out
  integer :: k, kk, kp2, ksize
  integer :: ipt(3, 12)  ! Pointers to number of coefficients for bodies
  integer :: lpt(3)      ! Pointer to number of coefficients for lunar librations
  integer :: rpt(3)      ! Pointer to number of coefficients for lunar euler angel rates
  integer :: tpt(3)      ! Pointer to number of coefficients for TT-TDB
  integer :: n, ncon, nrout, ncoeff, nrw, numde
  logical :: first = .true.


  ! By default, the output ephemeris will span the same interval as the
  ! input ascii data file(s).  The user may reset these to other JED's.
  !
  db2z = -99999999_8
  t1   = -99999999_8
  t2   =  99999999_8
  if (nrecl /= 1 .and. nrecl /= 4) then
    write (*,*) '*** ERROR: User did not set NRECL ***'
    stop
  end if

  ! Write a fingerprint to the screen.
  !
  write (*,*) ' JPL ASCII-TO-DIRECT-I/O program. ' // &
            & ' Last modified 15-Aug-2013.'

  ! Read the size and number of main ephemeris records.
  ! (from header.xxx)
  !
  read (*,'(6X,I6)') ksize
  write (*,'(/(A,I6))') 'KSIZE = ', ksize
  irecsz = nrecl * ksize

  ! Now for the alphameric heading records (GROUP 1010)
  !
  call nxtgrp(header)
  if (header /= 'GROUP   1010') then
    call errprt(1010, 'NOT HEADER')
  end if
  read (*,'(14A6)') ttl
  write (*,'(/(14A6))') ttl

  ! Read start, end and record span  (GROUP 1030)
  !
  call nxtgrp(header)
  if (header /= 'GROUP   1030') then
    call errprt(1030, 'NOT HEADER')
  end if
  read (*,'(3D12.0)') ss

  ! Read number of constants and names of constants (GROUP 1040/4).
  !
  call nxtgrp(header)
  if (header /= 'GROUP   1040 ') then
    call errprt(1040, 'NOT HEADER')
  end if
  read (*,'(I6)') n
  read (*,'(10A8)') (cnam(i), i=1,n)
  ncon = n

  ! Read number of values and values (GROUP 1041/4)
  !
  call nxtgrp(header)
  if (header /= 'GROUP   1041') then
    call errprt(1041, 'NOT HEADER')
  end if
  read (*,'(I6)') n
  read (*,*) (cval(i), i=1,n)
  do i = 1, n
    if (cnam(i) == 'AU    ') au    = cval(i)
    if (cnam(i) == 'EMRAT ') emrat = cval(i)
    if (cnam(i) == 'DENUM ') numde = cval(i)
  end do
  write (*,'(500(/2(A8,D24.16)))') (cnam(i), cval(i), i=1,n)

  ! Zero out pointer arrays
  !
  do i = 1, 3
    do j = 1, 12
      ipt(i, j) = 0
    end do
    lpt(i) = 0
    rpt(i) = 0
    tpt(i) = 0
  end do

  ! Read pointers needed by INTERP (GROUP 1050)
  !
  call nxtgrp(header)
  if (header /= 'GROUP   1050') then
    call errprt(1050, 'NOT HEADER')
  end if
  write (*,'(/)')
  do i = 1, 3
    read  (*,'(15I6)') (ipt(i, j), j=1,12), lpt(i), rpt(i), tpt(i)
    write (*,'(15I5)') (ipt(i, j), j=1,12), lpt(i), rpt(i), tpt(i)
  end do

  ! Open direct-access output file ('JPLEPH')

  open (unit   = 12,             &
      & file   = 'JPLEPH',       &
      & access = 'DIRECT',       &
      & form   = 'UNFORMATTED',  &
      & recl   = irecsz,         &
      & status = 'NEW')

  ! Read and write the ephemeris data records (GROUP 1070).
  !
  call nxtgrp(header)
  if (header /= 'GROUP   1070') then
    call errprt(1070, 'NOT HEADER')
  end if
  nrout = 0
  in    = 0
  out   = 0
  do
    read (*,'(2I6)') nrw, ncoeff
    if (nrw /= 0) exit
  end do
  do k = 1, ncoeff, 3
    kp2 = min(k + 2, ncoeff)
    read (*,*, iostat = in) (db(kk), kk=k,kp2)
    if (IN /= 0) stop ' Error reading 1st set of coeffs'
  end do

  do
    if ((in /= 0) .or. (db(2) >= t2)) exit
    if (2 * ncoeff /= ksize) then
      call errprt(ncoeff, ' 2*NCOEFF not equal to KSIZE')
    end if

    ! Skip this data block if the end of the interval is less
    ! than the specified start time or if the it does not begin
    ! where the previous block ended.
    !
    if  ((db(2) >= t1) .and. (db(1) >= db2z)) then
      if (first) then
        ! Don't worry about the intervals overlapping
        ! or abutting if this is the first applicable
        ! interval.
        !
        db2z  = db(1)
        first = .false.
      end if
      if (db(1) /= db2z) then
        ! Beginning of current interval is past the end
        ! of the previous one.
        !
        call errprt(nrw, 'Records do not overlap or abut')
      end if
      db2z  = db(2)
      nrout = nrout + 1
      write (12, rec=nrout+2, iostat=out) (db(k), k=1,ncoeff)
      if (out /= 0) then
        call errprt(nrout, 'th record not written because of error')
      end if

      ! Save this block's starting date, its interval span, and its end 
      ! date.
      !
      if (nrout == 1) then
        ss(1) = db(1)
        ss(3) = db(2) - db(1)
      end if
      ss(2) = db(2)

      ! Update the user as to our progress every 100th block.
      !
      if (mod(nrout, 100) == 1) then
        if (db(1) >= t1) then
          write (*,'(I6,A,F12.2)') &
            & nrout, ' EPHEMERIS RECORDS WRITTEN.  LAST JED = ', db(2)
        else
          write (*,*) ' Searching for first requested record...'
        end if
      end if
    end if
    read (*,'(2I6)', iostat=in) nrw, ncoeff
    if (in == 0)then
      do K = 1, ncoeff, 3
        kp2 = min(k + 2, ncoeff)
        read (*,*, iostat=in) (db(kk), kk=k,kp2)
        if (in /= 0) stop ' Error reading nth set of coeffs'
      end do
    end if
  end do
  write (*,'(I6,A,F12.2)') &
    & nrout, ' EPHEMERIS RECORDS WRITTEN.  LAST JED = ', db(2)

  ! Write header records onto output file.
  nrout = 1
  if (ncon <= oldmax) then
    write (12, rec=1, iostat=out) &
      & ttl, (cnam(i), i=1,oldmax), &
      & ss, ncon, au, emrat, ipt, numde, lpt, rpt, tpt
  else
    k = oldmax + 1
    write (12, rec=1, iostat=out) &
      & ttl, (cnam(i), i=1,oldmax), &
      & ss, ncon, au, emrat, ipt, numde, lpt, (cnam(j), j=k,ncon), rpt, tpt
  end if
  if (out /= 0 ) then
    call errprt(nrout, 'st record not written because of error')
  end if
  nrout = 2
  if (ncon <= oldmax) then
    write (12, rec=2, iostat=out) (cval(i), i=1,oldmax)
  else
    write (12, rec=2, iostat=out) (cval(i), i=1,ncon)
  end if
  if (out /= 0) then
    call errprt(nrout, 'nd record not written because of error')
  end if

  ! We're through.  Wrap it up.
  close (12)
  stop ' OK'

contains

  subroutine errprt(i, msg)
    implicit none
    integer, intent(IN) :: i
    character(len=*), intent(IN) :: msg

    write (*,'(A,I8,2X,A50)') 'ERROR #', i, msg

    stop ' ERROR '
  end subroutine errprt

  subroutine nxtgrp(header)
    implicit none
    character(len=*), intent(OUT) :: header
    character(len=12) :: blank

    ! Start with nothing.
    !
    header = ' '

    ! The next non-blank line we encounter is a header record.
    ! The group header and data are seperated by a blank line.
    !
    do
      if (header /= ' ') exit
      read (*,'(A)') header
    end do

    ! Found the header.  Read the blank line so we can get at the data.
    !
    if (header /= 'GROUP   1070') read (*,'(A)') blank

    return
  end subroutine nxtgrp
end program main

4. コンバートプログラムのコンパイル

$ gfortran -o asc2eph asc2eph.f95

データ検証プログラムもコンパイルしておく。(こちらは FORTRAN 77 のまま)

$ gfortran -o testeph testeph.f

5. テキスト形式データのマージ&コンバート

asc2eph をヘッダファイル・テキスト形式データファイルと同じディレクトリに配置して以下のように実行する。(必要なテキスト形式データファイルだけをマージするようにしてもよいし、全てのテキスト形式データファイルをマージするようにしてもよい)

$ cat header.430_572 ascp1950.430 ascp2050.430 | ./asc2eph
  JPL ASCII-TO-DIRECT-I/O program.  Last modified 15-Aug-2013.

KSIZE =   2036

JPL Planetary Ephemeris DE430/LE430
Start Epoch: JED=  2287184.5 1549-DEC-21 00:00:00
Final Epoch: JED=  2688976.5 2650-JAN-25 00:00:00

  DENUM   0.4300000000000000D+03  LENUM   0.4300000000000000D+03
  TDATEF  0.2013032920043800D+14  TDATEB  0.2013032919100700D+14
  JDEPOC  0.2440400500000000D+07  CENTER  0.0000000000000000D+00
  CLIGHT  0.2997924580000000D+06  BETA    0.1000000000000000D+01
  GAMMA   0.1000000000000000D+01  AU      0.1495978707000000D+09
  EMRAT   0.8130056907419062D+02  GM1     0.4912480450364760D-10
  GM2     0.7243452332644120D-09  GMB     0.8997011390199871D-09
  GM4     0.9549548695550771D-10  GM5     0.2825345840833870D-06

          ===< 中略 >===

  MA1093  0.2871640482670196D-15  MA1107  0.9549128215887647D-16
  MA1171  0.6248568708463761D-16  MA1467  0.1115280133034817D-15



    3  171  231  309  342  366  387  405  423  441  753  819  899    00
   14   10   13   11    8    7    6    6    6   13   11   10   10    00
    4    2    2    1    1    1    1    1    1    8    2    4    4    00
     1 EPHEMERIS RECORDS WRITTEN.  LAST JED =   2396784.50
   101 EPHEMERIS RECORDS WRITTEN.  LAST JED =   2399984.50

       ===< 中略 >===

  3401 EPHEMERIS RECORDS WRITTEN.  LAST JED =   2505584.50
  3425 EPHEMERIS RECORDS WRITTEN.  LAST JED =   2506352.50
STOP  OK
  • 成功すると、同じディレクトリ内に JPLEPH というバイナリファイルが作成される。(全てのテキスト形式データファイルをマージ&コンバートした場合は ftp://ssd.jpl.nasa.gov/pub/eph/planets/Linux/de430/linux_p1550p2650.430 と同じサイズになるはず)
  • JPLEPH が既に存在している場合はコンパイルエラーになるので、コンパイルし直す際は JPLEPH ファイルを削除してから。

6. バイナリ形式データの検証

testephJPLEPHtestpo.430 と同じディレクトリに配置して以下のように実行する。

$ ./testeph < testpo.430
  JPL TEST-EPHEMERIS program.  Last modified March 2013.

    2433264.50    2506352.50         32.00
  DENUM   0.4300000000000000D+03
  LENUM   0.4300000000000000D+03
  TDATEF  0.2013032920043800D+14
  TDATEB  0.2013032919100700D+14

===< 中略 >===

   line -- jed --   t#   c#   x#   --- jpl value ---   --- user value --    -- difference --

   100 2436294.5    4    9    3        -13.1973964533829        -13.1973964533829  0.17764E-14
   200 2439338.5    7    6    1        -27.5380053144894        -27.5380053144894  0.71054E-14
   300 2442382.5    4   11    6         -0.0039421481584         -0.0039421481584  0.26021E-17

===< 中略 >===

  2200 2500212.5    5    3    1         -1.9119755198047         -1.9119755198047  0.00000E+00
  2300 2503256.5   15    0    4          0.0003062237721          0.0003062237721  0.00000E+00
  2400 2506300.5    6    5    3          5.3862187413309          5.3862187413309  0.88818E-15

difference の列が 0.00000E+00 もしくは 0.88818E-15 のように限りなく小さい値になっていれば、データに誤りがない(テキストデータが正常にバイナリデータにコンバートされている)ということではないでしょうか。(本来は全ての difference 値が 0.00000E+00 になるべきであろう。この誤差は testpo.430 読み込み時の誤差によるもの(浮動小数点数が受ける影響によるもの))

また、全てのテキスト形式データファイルをマージ&コンバートした場合は、当方の環境では testeph 実行後に「浮動小数点」に関する例外が発生した。(原因は不明。これも浮動小数点数が受ける影響によるものかもしれない)

検証用フィル testpo.430 の仕様については、当ブログ過去記事「月・惑星の暦 JPL DE430 について!」を参照のこと。

7. 参考サイト

JPL 本家サイトと FTP サイト。


以上。

 

Sponsored Link

 

コメントする