Convert star catalog coordinates to position+velocity vector.
Status: support routine.
The star data accepted by this routine are "observables" for an imaginary observer at the solar-system barycenter. Proper motion and radial velocity are, strictly, in terms of barycentric coordinate time, TCB. For most practical applications, it is permissible to neglect the distinction between TCB and ordinary "proper" time on Earth (TT/TAI). The result will, as a rule, be limited by the intrinsic accuracy of the proper-motion and radial- velocity data; moreover, the pv-vector is likely to be merely an intermediate result, so that a change of time unit would cancel out overall.
In accordance with normal star-catalog conventions, the object's right ascension and declination are freed from the effects of secular aberration. The frame, which is aligned to the catalog equator and equinox, is Lorentzian and centered on the SSB.
The resulting position and velocity pv-vector is with respect to the same frame and, like the catalog coordinates, is freed from the effects of secular aberration. Should the "coordinate direction", where the object was located at the catalog epoch, be required, it may be obtained by calculating the magnitude of the position vector PV(1-3,1) dividing by the speed of light in au/day to give the light-time, and then multiplying the space velocity PV(1-3,2) by this light-time and adding the result to PV(1-3,1).
Summarizing, the pv-vector returned is for most stars almost identical to the result of applying the standard geometrical "space motion" transformation. The differences, which are the subject of the Stumpff paper referenced below, are:
(i) In stars with significant radial velocity and proper motion, the constantly changing light-time distorts the apparent proper motion. Note that this is a classical, not a relativistic, effect.
(ii) The transformation complies with special relativity.
Care is needed with units. The star coordinates are in radians and the proper motions in radians per Julian year, but the parallax is in arcseconds; the radial velocity is in km/s, but the pv-vector result is in au and au/day.
The RA proper motion is in terms of coordinate angle, not true angle. If the catalog uses arcseconds for both RA and Dec proper motions, the RA proper motion will need to be divided by cos(Dec) before use.
Straight-line motion at constant speed, in the inertial frame, is assumed.
An extremely small (or zero or negative) parallax is interpreted to mean that the object is on the "celestial sphere", the radius of which is an arbitrary (large) value (see the constant PXMIN). When the distance is overridden in this way, the status, initially zero, has 1 added to it.
If the space velocity is a significant fraction of c (see the constant VMAX), it is arbitrarily set to zero. When this action occurs, 2 is added to the status.
The relativistic adjustment involves an iterative calculation. If the process fails to converge within a set number (IMAX) of iterations, 4 is added to the status.
The inverse transformation is performed by the routine PVSTAR.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(in) | :: | ra | right ascension (radians) [see Note 1] |
||
| real(kind=wp), | intent(in) | :: | dec | declination (radians) [see Note 1] |
||
| real(kind=wp), | intent(in) | :: | pmr | RA proper motion (radians/year) [see Note 1] |
||
| real(kind=wp), | intent(in) | :: | pmd | Dec proper motion (radians/year) [see Note 1] |
||
| real(kind=wp), | intent(in) | :: | px | parallax (arcseconds) [see Note 1] |
||
| real(kind=wp), | intent(in) | :: | rv | radial velocity (km/s, positive = receding) [see Note 1] |
||
| real(kind=wp), | intent(out), | dimension(3,2) | :: | pv | pv-vector (au, au/day) [see Note 2] |
|
| integer, | intent(out) | :: | j | status [see Note 2]: * 0 = no warnings * 1 = distance overridden (Note 6) * 2 = excessive velocity (Note 7) * 4 = solution didn't converge (Note 8) * else = binary logical OR of the above |
subroutine STARPV ( ra, dec, pmr, pmd, px, rv, pv, j )
implicit none
real(wp),intent(in) :: ra !! right ascension (radians) [see Note 1]
real(wp),intent(in) :: dec !! declination (radians) [see Note 1]
real(wp),intent(in) :: pmr !! RA proper motion (radians/year) [see Note 1]
real(wp),intent(in) :: pmd !! Dec proper motion (radians/year) [see Note 1]
real(wp),intent(in) :: px !! parallax (arcseconds) [see Note 1]
real(wp),intent(in) :: rv !! radial velocity (km/s, positive = receding) [see Note 1]
real(wp),dimension(3,2),intent(out) :: pv !! pv-vector (au, au/day) [see Note 2]
integer,intent(out) :: j !! status [see Note 2]:
!! * 0 = no warnings
!! * 1 = distance overridden (Note 6)
!! * 2 = excessive velocity (Note 7)
!! * 4 = solution didn't converge (Note 8)
!! * else = binary logical OR of the above
! Smallest allowed parallax
real(wp),parameter :: pxmin = 1.0e-7_wp
! Largest allowed speed (fraction of c)
real(wp),parameter :: vmax = 0.5_wp
! Julian years to days
real(wp),parameter :: y2d = 365.25_wp
! Radians to arcseconds
real(wp),parameter :: dr2as = 206264.8062470963551564734_wp
! Astronomical unit (m, IAU 2012)
real(wp),parameter :: aum = 149597870.7e3_wp
! Speed of light (au per day)
real(wp),parameter :: c = d2s*cmps/aum
! Maximum number of iterations for relativistic solution
integer,parameter :: imax = 100
integer :: i
integer :: iwarn
real(wp) :: w, r, rd, rad, decd, v, x(3), usr(3), ust(3), &
vsr, vst, betst, betsr, bett, betr, od, odel, &
dd, ddel, odd, oddel, d, del, ur(3), ut(3)
! Distance (au).
if ( px>=pxmin ) then
w = px
iwarn = 0
else
w = pxmin
iwarn = 1
end if
r = dr2as / w
! Radial velocity (au/day).
rd = d2s * rv * 1.0e3_wp / aum
! Proper motion (radian/day).
rad = pmr / y2d
decd = pmd / y2d
! To pv-vector (au,au/day).
call S2PV ( ra, dec, r, rad, decd, rd, pv )
! If excessive velocity, arbitrarily set it to zero.
call PM ( pv(1,2), v )
if ( v/c > vmax ) then
call ZP ( pv(1,2) )
iwarn = iwarn + 2
end if
! Isolate the radial component of the velocity (au/day).
call PN ( pv(1,1), w, x )
call PDP ( x, pv(1,2), vsr )
call SXP ( vsr, x, usr )
! Isolate the transverse component of the velocity (au/day).
call PMP ( pv(1,2), usr, ust )
call PM ( ust, vst )
! Special-relativity dimensionless parameters.
betsr = vsr / c
betst = vst / c
! Determine the inertial-to-observed relativistic correction terms.
od = 0.0_wp
odel = 0.0_wp
odd = 0.0_wp
oddel = 0.0_wp
bett = betst
betr = betsr
do i=1,imax
d = 1.0_wp + betr
w = betr*betr + bett*bett
del = - w / ( sqrt(1.0_wp-w) + 1.0_wp )
betr = d*betsr + del
bett = d*betst
if ( i > 1 ) then
dd = abs(d-od)
ddel = abs(del-odel)
if ( i>2 .and. &
dd>=odd .and. &
ddel>=oddel ) exit
if ( i >= imax ) iwarn = iwarn + 4
odd = dd
oddel = ddel
end if
od = d
odel = del
end do
! Replace observed radial velocity with inertial value.
if ( betsr /= 0.0_wp ) then
w = d + del/betsr
else
w = 1.0_wp
end if
call SXP ( w, usr, ur )
! Replace observed tangential velocity with inertial value.
call SXP ( d, ust, ut )
! Combine the two to obtain the inertial space velocity.
call PPP ( ur, ut, pv(1,2) )
! Return the status.
j = iwarn
end subroutine STARPV