ephemeris_test Subroutine

public subroutine ephemeris_test()

Ephemeris test routine.

Ephemeris files

This routine requires the DE405 and DE421 JPL binary ephemeris files to be present in the ./eph directory. These can be built by using the instructions here.

For example (on Linux):

  wget ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran/*
  wget ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii/de405/*
  #edit asc2eph.f file to set NRECL = 4:
  sed -i '_original' '/^C.*PARAMETER ( NRECL = 4 )/s/^C//' asc2eph.f
  gfortran asc2eph.f -o asc2eph
  cat header.405 ascp*.405 | ./asc2eph
  mkdir Fortran-Astrodynamics-Toolkit/eph
  mv JPLEPH Fortran-Astrodynamics-Toolkit/eph/JPLEPH.405

Arguments

None

Calls

proc~~ephemeris_test~~CallsGraph proc~ephemeris_test jpl_ephemeris_module::ephemeris_test proc~get_constants jpl_ephemeris_module::jpl_ephemeris%get_constants proc~ephemeris_test->proc~get_constants proc~get_state jpl_ephemeris_module::jpl_ephemeris%get_state proc~ephemeris_test->proc~get_state proc~initialize_ephemeris jpl_ephemeris_module::jpl_ephemeris%initialize_ephemeris proc~ephemeris_test->proc~initialize_ephemeris proc~state jpl_ephemeris_module::jpl_ephemeris%state proc~get_state->proc~state proc~interp jpl_ephemeris_module::jpl_ephemeris%interp proc~state->proc~interp proc~split jpl_ephemeris_module::split proc~state->proc~split

Source Code

    subroutine ephemeris_test()

    implicit none

    character(len=6),dimension(nmax) :: nams
    real(wp) :: jd
    real(wp),dimension(6) :: rv,rv1,rv2,diffrv
    real(wp),dimension(3) :: ss
    real(wp),dimension(nmax) :: vals
    integer :: nvs,ntarg,nctr,i,j
    type(jpl_ephemeris) :: eph405, eph421
    logical :: status_ok_405,status_ok_421

    character(len=*),parameter :: ephemeris_file_405 = './eph/JPLEPH.405' !! JPL DE405 ephemeris file
    character(len=*),parameter :: ephemeris_file_421 = './eph/JPLEPH.421' !! JPL DE421 ephemeris file

    write(*,*) ''
    write(*,*) '---------------'
    write(*,*) ' ephemeris_test'
    write(*,*) '---------------'
    write(*,*) ''

    !initialize:
    call eph405%initialize(filename=ephemeris_file_405,status_ok=status_ok_405)
    call eph421%initialize(filename=ephemeris_file_421,status_ok=status_ok_421)

    if (status_ok_405) then

        !get some constants from the file:
        call eph405%get_constants(nams,vals,ss,nvs)

        write(*,'(A)') ''
        write(*,'(A)') 'Ephemeris initialized'
        write(*,'(A,1X,F15.3,1X,A,1X,F15.3)') 'JD range: ',ss(1), 'to ', ss(2)
        write(*,'(A)') ''
        do i=1,nvs
            write(*,'(A,1X,D25.16)') nams(i), vals(i)
        end do

        jd = 2451536.5d0       ! julian date

        if (jd < ss(1) .or. jd > ss(2)) then

            write(*,'(A)') ''
            write(*,*) 'error: jed out of bounds.'
            write(*,*) 'jed   = ', jd
            write(*,*) 'ss(1) = ', ss(1)
            write(*,*) 'ss(2) = ', ss(2)

        else

            !test DE405:
            do j=1,2

                if (j==1) then
                    ntarg = 3      !earth
                    nctr  = 11     !sun
                else
                    ntarg = 10     !moon
                    nctr  = 3      !earth
                end if

                write(*,*) ''
                write(*,*) 'DE405'
                write(*,*) 'state of "'//trim(list_of_bodies(ntarg))//&
                            '" wrt "'//trim(list_of_bodies(nctr))//'"'

                do i=1,10
                    call eph405%get_state( jd, ntarg, nctr, rv, status_ok_405)
                    write(*,'(F15.2,1X,*(E25.16,1X))') jd, norm2(rv(1:3)), rv
                    jd = jd + 10.0_wp
                end do

            end do
        end if

    else
        write(*,*) 'Error opening DE405 ephemeris file'
    end if

    if (status_ok_405 .and. status_ok_421) then

        !compare DE405 with DE421
        do j=1,2

            if (j==1) then
                ntarg = 3      !earth
                nctr  = 11     !sun
            else
                ntarg = 10     !moon
                nctr  = 3      !earth
            end if

            write(*,*) ''
            write(*,*) 'DE421 - DE405 Difference'
            write(*,*) 'state of "'//trim(list_of_bodies(ntarg))//&
                        '" wrt "'//trim(list_of_bodies(nctr))//'"'

            do i=1,10
                call eph405%get_state( jd, ntarg, nctr, rv1, status_ok_405)
                call eph421%get_state( jd, ntarg, nctr, rv2, status_ok_421)
                diffrv = rv2 - rv1
                write(*,'(F15.2,1X,*(E25.16,1X))') jd, norm2(diffrv(1:3)), norm2(diffrv(4:6))
                jd = jd + 10.0_wp
            end do

        end do

    else
        write(*,*) 'Error opening DE421 ephemeris file'
    end if

    !cleanup:
    call eph405%close()
    call eph421%close()

    end subroutine ephemeris_test