standish_module_test Subroutine

    subroutine standish_module_test()

    implicit none

    type(standish_ephemeris) :: eph
    real(wp)              :: et         !! ephemeris time (sec)
    real(wp)              :: jd         !! julian date
    real(wp),dimension(6) :: rv         !! state of targ w.r.t. obs
    logical               :: status_ok  !! true if there were no problems
    integer :: itbl

    !>
    !  State of Earth w.r.t. Sun from SPICE
    !  See: http://wgc.jpl.nasa.gov:8080/webgeocalc/#NewCalculation
    real(wp),dimension(6),parameter :: rv_from_spice = [-26502576.96907434_wp,&
                                                        132754176.58943012_wp,&
                                                        57555793.70952077_wp,&
                                                        -29.78644078_wp,&
                                                        -5.02614568_wp,&
                                                        -2.17905509_wp]

    write(*,*) ''
    write(*,*) '---------------'
    write(*,*) ' standish_module_test'
    write(*,*) '---------------'
    write(*,*) ''

    et = 0.0_wp ! J2000

    ! get earth w.r.t. sun in J2000 frame:
    call eph%get_rv(et,body_earth_moon_barycenter,body_sun,rv,status_ok)

    if (status_ok) then
        write(*,*) ''
        write(*,*) 'State of Earth wrt Sun @ J2000'
        write(*,'(A,1X,*(E26.16,1X))') 'standish:', rv
        write(*,'(A,1X,*(E26.16,1X))') 'SPICE:   ', rv_from_spice
        write(*,'(A,1X,*(E26.16,1X))') 'diff:    ', rv - rv_from_spice
        write(*,*) ''
    else
        error stop 'error calling standish_ephemeris'
    end if

    ! low-level routine tests:
    write(*,*) 'helio:'
    jd = et_to_jd(et)
    call helio (3, jd, rv, itbl)
    write(*,*) ''
    write(*,*) 'State of Earth wrt Sun [ecliptic]'
    write(*,'(*(E26.16,1X))') rv
    write(*,*) ''

    end subroutine standish_module_test