Skip to content

Instantly share code, notes, and snippets.

@komasaru
Created October 14, 2018 01:49
Show Gist options
  • Save komasaru/447a920345e1579102d6f0ff712d45cd to your computer and use it in GitHub Desktop.
Save komasaru/447a920345e1579102d6f0ff712d45cd to your computer and use it in GitHub Desktop.
Fortran 95 source code to convert JPL source cord from ascii to binary
!***********************************************************************
! 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment