This subroutine reads the JPL planetary ephemeris
and gives the position and velocity of the point ntarg
with respect to ncent.
The numbering convention for ntarg and ncent is:
1 = mercury 8 = neptune
2 = venus 9 = pluto
3 = earth 10 = moon
4 = mars 11 = sun
5 = jupiter 12 = solar-system barycenter
6 = saturn 13 = earth-moon barycenter
7 = uranus 14 = nutations (longitude and obliq)
15 = librations, if on eph file
(if nutations are wanted, set ntarg = 14. for librations, set ntarg = 15. set ncent=0.)
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(jpl_ephemeris), | intent(inout) | :: | me | |||
| real(kind=wp), | intent(in) | :: | jd |
d.p. Julian ephemeris date at which interpolation is wanted. |
||
| integer, | intent(in) | :: | ntarg |
integer number of 'target' point. |
||
| integer, | intent(in) | :: | ncent |
integer number of 'center' point. |
||
| real(kind=wp), | intent(out), | dimension(6) | :: | rrd |
output 6-word d.p. array containing position and velocity
of point |
|
| logical, | intent(out) | :: | status_ok |
true if there were no problems |
subroutine get_state(me,jd,ntarg,ncent,rrd,status_ok) implicit none class(jpl_ephemeris),intent(inout) :: me real(wp),intent(in) :: jd !! d.p. Julian ephemeris date at which interpolation is wanted. integer,intent(in) :: ntarg !! integer number of 'target' point. integer,intent(in) :: ncent !! integer number of 'center' point. real(wp),dimension(6),intent(out) :: rrd !! output 6-word d.p. array containing position and velocity !! of point `ntarg` relative to `ncent`. !! the units are AU and AU/day (or km and km/sec if `me%km=.true.`). !! For librations the units are radians and radians !! per day. In the case of nutations the first four words of !! `rrd` will be set to nutations and rates, having units of !! radians and radians/day. logical,intent(out) :: status_ok !! true if there were no problems real(wp),dimension(2) :: et2 real(wp),dimension(6,13) :: pv real(wp),dimension(6,11) :: pvst real(wp),dimension(4) :: pnut integer,dimension(12) :: list integer :: i,j,k logical :: bsave status_ok = .false. if (me%initialized) then ! initialize et2 for 'state' and set up component count et2(1) = jd et2(2) = 0.0_wp list = 0 rrd = 0.0_wp if (ntarg /= ncent) then select case (ntarg) case (14) !nutation if (me%ipt(2,12)>0) then !me%ipt(35) list(11) = 2 call me%state(et2,list,pvst,pnut,status_ok) if (.not. status_ok) then write(error_unit,'(A)') 'error interpolating nutations in get_state.' return end if rrd(1:4) = pnut(1:4) rrd(5) = 0.0_wp rrd(6) = 0.0_wp return else rrd(1:4) = 0.0_wp write(error_unit,'(A)') 'error in get_state: the ephemeris file does not contain nutations.' return endif case (15) !librations if (me%ipt(2,13)>0) then !me%ipt(38) list(12) = 2 call me%state(et2,list,pvst,pnut,status_ok) if (.not. status_ok) then write(error_unit,'(A)') 'error interpolating librations in get_state.' return end if rrd = pvst(:,11) return else write(error_unit,'(A)') 'error in get_state: the ephemeris file does not contain librations.' return endif case default ! force barycentric output by 'state' bsave = me%bary me%bary = .true. ! set up proper entries in 'list' array for state call do i=1,2 k=ntarg if (i == 2) k = ncent if (k <= 10) list(k) = 2 if (k == 10) list(3) = 2 if (k == 3) list(10) = 2 if (k == 13) list(3) = 2 enddo ! make call to state call me%state(et2,list,pvst,pnut,status_ok) if (.not. status_ok) then write(error_unit,'(A)') 'error interpolating state in get_state.' return end if do i=1,10 do j = 1,6 pv(j,i) = pvst(j,i) end do enddo if (ntarg == 11 .or. ncent == 11) pv(:,11) = me%pvsun if (ntarg == 12 .or. ncent == 12) pv(:,12) = 0.0_wp if (ntarg == 13 .or. ncent == 13) pv(:,13) = pvst(:,3) if (ntarg*ncent == 30 .and. ntarg+ncent == 13) then pv(:,3) = 0.0_wp else if (list(3) == 2) pv(:,3) = pvst(:,3) - pvst(:,10)/(1.0_wp+me%emrat) if (list(10) == 2) pv(:,10) = pv(:,3) + pvst(:,10) end if rrd = pv(:,ntarg) - pv(:,ncent) me%bary = bsave end select end if status_ok = .true. else write(error_unit,'(A)') 'error in get_state: the ephemeris is not initialized.' end if end subroutine get_state