state Subroutine

private subroutine state(me, et2, list, pv, pnut, status_ok)

This subroutine reads and interpolates the JPL planetary ephemeris file.

Note

The ephemeris is assumed to have been initialized.

Type Bound

jpl_ephemeris

Arguments

Type IntentOptional Attributes Name
class(jpl_ephemeris), intent(inout) :: me
real(kind=wp), intent(in), dimension(2) :: et2

2-word Julian ephemeris epoch at which interpolation is wanted. any combination of et2(1)+et2(2) which falls within the time span on the file is a permissible epoch.

a. for ease in programming, the user may put the entire epoch in et2(1) and set et2(2)=0. b. for maximum interpolation accuracy, set et2(1) = the most recent midnight at or before interpolation epoch and set et2(2) = fractional part of a day elapsed between et2(1) and epoch. c. as an alternative, it may prove convenient to set et2(1) = some fixed epoch, such as start of integration, and et2(2) = elapsed interval between then and epoch.
integer, intent(in), dimension(12) :: list

12-word integer array specifying what interpolation is wanted for each of the bodies on the file. list(i) = 0 : no interpolation for body i, list(i) = 1 : position only, list(i) = 2 : position and velocity.

The designation of the astronomical bodies by i is: i = 1 : mercury, i = 2 : venus, i = 3 : earth-moon barycenter, i = 4 : mars, i = 5 : jupiter, i = 6 : saturn, i = 7 : uranus, i = 8 : neptune, i = 9 : pluto, i =10 : geocentric moon, i =11 : nutations in longitude and obliquity, i =12 : lunar librations (if on file).
real(kind=wp), intent(out), dimension(6,11) :: pv

dp 6 x 11 array that will contain requested interpolated quantities (other than nutation, stored in pnut). the body specified by list(i) will have its state in the array starting at pv(1,i). (on any given call, only those words in pv which are affected by the first 10 list entries, and by list(12) if librations are on the file, are set. the rest of the pv array is untouched.) the order of components starting in pv(1,i) is: x,y,z,dx,dy,dz.

All output vectors are referenced to the earth mean equator and equinox of J2000 if the DE number is 200 or greater; of B1950 if the DE number is less than 200.

The moon state is always geocentric; the other nine states are either heliocentric or solar-system barycentric, depending on the setting of the bary variable in the class.

Lunar librations, if on file, are put into pv(k,11) if list(12) is 1 or 2.
real(kind=wp), intent(out), dimension(4) :: pnut

dp 4-word array that will contain nutations and rates, depending on the setting of list(11). the order of quantities in pnut is:

  • d psi (nutation in longitude),
  • d epsilon (nutation in obliquity),
  • d psi dot,
  • d epsilon dot.
logical, intent(out) :: status_ok

true if there were no problems

Calls

proc~~state~~CallsGraph proc~state jpl_ephemeris_module::jpl_ephemeris%state proc~interp jpl_ephemeris_module::jpl_ephemeris%interp proc~state->proc~interp proc~split jpl_ephemeris_module::split proc~state->proc~split

Called by

proc~~state~~CalledByGraph proc~state jpl_ephemeris_module::jpl_ephemeris%state proc~get_state jpl_ephemeris_module::jpl_ephemeris%get_state proc~get_state->proc~state proc~ephemeris_test jpl_ephemeris_module::ephemeris_test proc~ephemeris_test->proc~get_state proc~get_rv_from_jpl_ephemeris jpl_ephemeris_module::jpl_ephemeris%get_rv_from_jpl_ephemeris proc~get_rv_from_jpl_ephemeris->proc~get_state

Source Code

    subroutine state(me,et2,list,pv,pnut,status_ok)

    implicit none

    class(jpl_ephemeris),intent(inout)   :: me
    real(wp),dimension(2),intent(in)     :: et2     !! 2-word Julian ephemeris epoch at which interpolation
                                                    !! is wanted.  any combination of `et2(1)+et2(2)` which falls
                                                    !! within the time span on the file is a permissible epoch.
                                                    !!
                                                    !! ***a.*** for ease in programming, the user may put the
                                                    !!    entire epoch in `et2(1)` and set `et2(2)=0`.
                                                    !! ***b.*** for maximum interpolation accuracy, set `et2(1)` =
                                                    !!    the most recent midnight at or before interpolation
                                                    !!    epoch and set `et2(2)` = fractional part of a day
                                                    !!    elapsed between `et2(1)` and epoch.
                                                    !! ***c.*** as an alternative, it may prove convenient to set
                                                    !!    `et2(1)` = some fixed epoch, such as start of integration,
                                                    !!    and `et2(2)` = elapsed interval between then and epoch.

    integer,dimension(12),intent(in)     :: list    !! 12-word integer array specifying what interpolation
                                                    !! is wanted for each of the bodies on the file.
                                                    !!    ***list(i) = 0*** : no interpolation for body i,
                                                    !!    ***list(i) = 1*** : position only,
                                                    !!    ***list(i) = 2*** : position and velocity.
                                                    !!
                                                    !! The designation of the astronomical bodies by i is:
                                                    !!    ***i = 1*** : mercury,
                                                    !!    ***i = 2*** : venus,
                                                    !!    ***i = 3*** : earth-moon barycenter,
                                                    !!    ***i = 4*** : mars,
                                                    !!    ***i = 5*** : jupiter,
                                                    !!    ***i = 6*** : saturn,
                                                    !!    ***i = 7*** : uranus,
                                                    !!    ***i = 8*** : neptune,
                                                    !!    ***i = 9*** : pluto,
                                                    !!    ***i =10*** : geocentric moon,
                                                    !!    ***i =11*** : nutations in longitude and obliquity,
                                                    !!    ***i =12*** : lunar librations (if on file).

    real(wp),dimension(6,11),intent(out) :: pv      !! dp 6 x 11 array that will contain requested interpolated
                                                    !! quantities (other than nutation, stored in `pnut`).
                                                    !! the body specified by `list(i)` will have its
                                                    !! state in the array starting at `pv(1,i)`.
                                                    !! (on any given call, only those words in `pv` which are
                                                    !! affected by the first 10 `list` entries, and by `list(12)`
                                                    !! if librations are on the file, are set.
                                                    !! the rest of the `pv` array is untouched.)
                                                    !! the order of components starting in `pv(1,i)` is:
                                                    !! `x`,`y`,`z`,`dx`,`dy`,`dz`.
                                                    !!
                                                    !! All output vectors are referenced to the earth mean
                                                    !! equator and equinox of J2000 if the DE number is 200 or
                                                    !! greater; of B1950 if the DE number is less than 200.
                                                    !!
                                                    !! The moon state is always geocentric; the other nine states
                                                    !! are either heliocentric or solar-system barycentric,
                                                    !! depending on the setting of the `bary` variable in the class.
                                                    !!
                                                    !! Lunar librations, if on file, are put into `pv(k,11)` if
                                                    !! `list(12)` is `1` or `2`.

    real(wp),dimension(4),intent(out)    :: pnut    !! dp 4-word array that will contain nutations and rates,
                                                    !! depending on the setting of `list(11)`.  the order of
                                                    !! quantities in `pnut` is:
                                                    !!
                                                    !! * `d psi`  (nutation in longitude),
                                                    !! * `d epsilon` (nutation in obliquity),
                                                    !! * `d psi dot`,
                                                    !! * `d epsilon dot`.
    logical,intent(out) :: status_ok !! true if there were no problems

    real(wp),dimension(2)    :: t
    real(wp),dimension(4)    :: pjd
    real(wp) :: aufac,s,tmp1,tmp2
    integer :: istat,i,j,k,nr

    status_ok = .true.

    s = et2(1)-0.5_wp
    call split(s,pjd(1:2))
    call split(et2(2),pjd(3:4))
    pjd(1) = pjd(1)+pjd(3)+0.5_wp
    pjd(2) = pjd(2)+pjd(4)
    call split(pjd(2),pjd(3:4))
    pjd(1) = pjd(1)+pjd(3)

    ! error return for epoch out of range

    if (pjd(1)+pjd(4)<me%ss(1) .or. pjd(1)+pjd(4)>me%ss(2)) then
        write(error_unit,'(A,F12.2,A,2F22.2)') &
                'Error: requested jed,',&
                et2(1)+et2(2),&
                ' not within ephemeris limits,',&
                me%ss(1),me%ss(2)
        status_ok = .false.
        return
    end if

    ! calculate record # and relative time in interval

    nr = int((pjd(1)-me%ss(1))/me%ss(3))+3
    if (pjd(1)==me%ss(2)) nr = nr-1

    tmp1 = dble(nr-3)*me%ss(3) + me%ss(1)
    tmp2 = pjd(1) - tmp1
    t(1) = (tmp2 + pjd(4))/me%ss(3)

    ! read correct record if not in core

    if (nr/=me%nrl) then
        me%nrl = nr
        read(me%nrfile,rec=nr,iostat=istat) (me%buf(k),k=1,me%ncoeffs)
        if (istat/=0) then
            write(error_unit,'(2F12.2,A80)') et2,'Error return in state'
            status_ok = .false.
            return
        end if
    endif

    if (me%km) then
        t(2) = me%ss(3)*86400.0_wp
        aufac = 1.0_wp
    else
        t(2) = me%ss(3)
        aufac = 1.0_wp/me%au
    endif

    ! interpolate ssbary sun

    call me%interp(me%buf(me%ipt(1,11)),t,me%ipt(2,11),3,me%ipt(3,11),2,me%pvsun)

    me%pvsun = me%pvsun * aufac

    ! check and interpolate whichever bodies are requested

    do i=1,10
        if (list(i)==0) cycle
        call me%interp(me%buf(me%ipt(1,i)),t,me%ipt(2,i),3,me%ipt(3,i),list(i),pv(1,i))
        do j=1,6
            if (i<=9 .and. .not.me%bary) then
                pv(j,i) = pv(j,i)*aufac - me%pvsun(j)
            else
                pv(j,i) = pv(j,i)*aufac
            endif
        enddo
    end do

    ! do nutations if requested (and if on file)

    if (list(11)>0 .and. me%ipt(2,12)>0) &
        call me%interp(me%buf(me%ipt(1,12)),t,me%ipt(2,12),2,me%ipt(3,12),list(11),pnut)

    ! get librations if requested (and if on file)

    if (list(12)>0 .and. me%ipt(2,13)>0) &
        call me%interp(me%buf(me%ipt(1,13)),t,me%ipt(2,13),3,me%ipt(3,13),list(12),pv(1,11))

    end subroutine state