calcelements Subroutine

private pure subroutine calcelements(np, jd, itbl, z)

Calculate current elements z(jd) for planet j from jpl data

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: np

planet number (1-9)
real(kind=wp), intent(in) :: jd

julian date
integer, intent(in) :: itbl

which table to use (1-2)
real(kind=wp), intent(out), dimension(8) :: z

elements for jd

  • a = semi major axis (au)
  • e = eccentricity (rad)
  • i = inclination (rad)
  • l = mean longitude (rad)
  • w = longitude of perihelion (rad)
  • o = longitude of ascending mode (rad)
  • ma = mean anomaly (rad)
  • ea = eccentric anomaly (rad)

Calls

proc~~calcelements~~CallsGraph proc~calcelements standish_module::calcelements proc~kepler standish_module::kepler proc~calcelements->proc~kepler

Called by

proc~~calcelements~~CalledByGraph proc~calcelements standish_module::calcelements proc~helio standish_module::helio proc~helio->proc~calcelements proc~standish_module_test standish_module::standish_module_test proc~standish_module_test->proc~helio proc~standish_rv_func standish_module::standish_ephemeris%standish_rv_func proc~standish_module_test->proc~standish_rv_func proc~standish_rv_func->proc~helio

Source Code

    pure subroutine calcelements (np, jd, itbl, z)

    implicit none

    integer, intent (in) :: np        !! planet number (1-9)
    integer, intent (in) :: itbl      !! which table to use (1-2)
    real(wp), intent (in) :: jd       !! julian date
    real(wp), dimension(8), intent(out) :: z  !! elements for `jd`
                                              !!
                                              !!  * a = semi major axis (au)
                                              !!  * e = eccentricity (rad)
                                              !!  * i = inclination (rad)
                                              !!  * l = mean longitude (rad)
                                              !!  * w = longitude of perihelion (rad)
                                              !!  * o = longitude of ascending mode (rad)
                                              !!  * ma = mean anomaly (rad)
                                              !!  * ea = eccentric anomaly (rad)

    integer :: i   !! counter
    real(wp) :: t  !! centuries since epoch
    real(wp) :: tz !! perturbation term

    t = (jd-epoch) * day2century

    do i = 1, 6      ! a,e,i,l,w,o
        z (i) = eph_set(itbl)%o(i, np) + eph_set(itbl)%o(i+6, np) * t
        ! if (i>2) z(i) = modulo(z(i), twopi)  !optional scaling
    end do

    if (itbl==2) then
        ! perturbation term nonzero for planets 5-9 if table 2 used
        tz =  eph_set(itbl)%o(13, np) * t ** 2 + &
              eph_set(itbl)%o(14, np) * cos(eph_set(itbl)%o(16, np)*t) + &
              eph_set(itbl)%o(15, np) * sin(eph_set(itbl)%o(16, np)*t)
    else
        tz = zero
    end if

    z (7) = modulo ((z(4)-z(5)+tz), twopi)  ! mean anomaly
    z (8) = kepler (z(7), z(2))             ! eccentric anomaly

    end subroutine calcelements