Skip to content

Instantly share code, notes, and snippets.

@komasaru komasaru/asc2eph.f95
Created Oct 14, 2018

Embed
What would you like to do?
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
You can’t perform that action at this time.