Quick CIRS to observed place transformation.
Use of this routine is appropriate when efficiency is important and where many star positions are all to be transformed for one date. The star-independent astrometry parameters can be obtained by calling APIO[13] or APCO[13].
Status: support routine.
This routine returns zenith distance rather than altitude in order to reflect the fact that no allowance is made for depression of the horizon.
The accuracy of the result is limited by the corrections for refraction, which use a simple Atan(z) + Btan^3(z) model. Providing the meteorological parameters are known accurately and there are no gross local effects, the predicted observed coordinates should be within 0.05 arcsec (optical) or 1 arcsec (radio) for a zenith distance of less than 70 degrees, better than 30 arcsec (optical or radio) at 85 degrees and better than 20 arcmin (optical) or 30 arcmin (radio) at the horizon.
Without refraction, the complementary routines ATIOQ and ATOIQ are self-consistent to better than 1 microarcsecond all over the celestial sphere. With refraction included, consistency falls off at high zenith distances, but is still better than 0.05 arcsec at 85 degrees.
It is advisable to take great care with units, as even unlikely values of the input parameters are accepted and processed in accordance with the models used.
The CIRS RA,Dec is obtained from a star catalog mean place by allowing for space motion, parallax, the Sun's gravitational lens effect, annual aberration and precession-nutation. For star positions in the ICRS, these effects can be applied by means of the ATCI13 (etc.) routines. Starting from classical "mean place" systems, additional transformations will be needed first.
"Observed" Az,El means the position that would be seen by a perfect geodetically aligned theodolite. This is obtained from the CIRS RA,Dec by allowing for Earth orientation and diurnal aberration, rotating from equator to horizon coordinates, and then adjusting for refraction. The HA,Dec is obtained by rotating back into equatorial coordinates, and is the position that would be seen by a perfect equatorial with its polar axis aligned to the Earth's axis of rotation. Finally, the RA is obtained by subtracting the HA from the local ERA.
The star-independent CIRS-to-observed-place parameters in ASTROM may be computed with APIO[13] or APCO[13]. If nothing has changed significantly except the time, APER[13] may be used to perform the requisite adjustment to the ASTROM array.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(in) | :: | ri | CIRS right ascension |
||
| real(kind=wp), | intent(in) | :: | di | CIRS declination |
||
| real(kind=wp), | intent(in), | dimension(30) | :: | astrom | star-independent astrometry parameters: (1) PM time interval (SSB, Julian years) (2-4) SSB to observer (vector, au) (5-7) Sun to observer (unit vector) (8) distance from Sun to observer (au) (9-11) v: barycentric observer velocity (vector, c) (12) sqrt(1-|v|^2): reciprocal of Lorenz factor (13-21) bias-precession-nutation matrix (22) longitude + s' (radians) (23) polar motion xp wrt local meridian (radians) (24) polar motion yp wrt local meridian (radians) (25) sine of geodetic latitude (26) cosine of geodetic latitude (27) magnitude of diurnal aberration vector (28) "local" Earth rotation angle (radians) (29) refraction constant A (radians) (30) refraction constant B (radians) |
|
| real(kind=wp), | intent(out) | :: | aob | observed azimuth (radians: N=0,E=90) |
||
| real(kind=wp), | intent(out) | :: | zob | observed zenith distance (radians) |
||
| real(kind=wp), | intent(out) | :: | hob | observed hour angle (radians) |
||
| real(kind=wp), | intent(out) | :: | dob | observed declination (CIO-based, radians) |
||
| real(kind=wp), | intent(out) | :: | rob | observed right ascension (CIO-based, radians) |
subroutine ATIOQ ( ri, di, astrom, aob, zob, hob, dob, rob )
implicit none
real(wp),intent(in) :: ri !! CIRS right ascension
real(wp),intent(in) :: di !! CIRS declination
real(wp),dimension(30),intent(in) :: astrom !! star-independent astrometry parameters:
!! (1) PM time interval (SSB, Julian years)
!! (2-4) SSB to observer (vector, au)
!! (5-7) Sun to observer (unit vector)
!! (8) distance from Sun to observer (au)
!! (9-11) v: barycentric observer velocity (vector, c)
!! (12) sqrt(1-|v|^2): reciprocal of Lorenz factor
!! (13-21) bias-precession-nutation matrix
!! (22) longitude + s' (radians)
!! (23) polar motion xp wrt local meridian (radians)
!! (24) polar motion yp wrt local meridian (radians)
!! (25) sine of geodetic latitude
!! (26) cosine of geodetic latitude
!! (27) magnitude of diurnal aberration vector
!! (28) "local" Earth rotation angle (radians)
!! (29) refraction constant A (radians)
!! (30) refraction constant B (radians)
real(wp),intent(out) :: aob !! observed azimuth (radians: N=0,E=90)
real(wp),intent(out) :: zob !! observed zenith distance (radians)
real(wp),intent(out) :: hob !! observed hour angle (radians)
real(wp),intent(out) :: dob !! observed declination (CIO-based, radians)
real(wp),intent(out) :: rob !! observed right ascension (CIO-based, radians)
! Minimum sine and cosine of altitude for refraction purposes
real(wp),parameter :: selmin = 0.05_wp
real(wp),parameter :: celmin = 1.0e-6_wp
real(wp) :: v(3), x, y, z, xhd, yhd, zhd, f, &
xhdt, yhdt, zhdt, xaet, yaet, zaet, azobs, &
r, tz, w, del, cosdel, xaeo, yaeo, zaeo, &
zdobs, hmobs, dcobs, raobs
! CIRS RA,Dec to Cartesian -HA,Dec.
call S2C ( ri-astrom(28), di, v )
x = v(1)
y = v(2)
z = v(3)
! Polar motion.
xhd = x + astrom(23)*z
yhd = y - astrom(24)*z
zhd = z - astrom(23)*x + astrom(24)*y
! Diurnal aberration.
f = ( 1.0_wp - astrom(27)*yhd )
xhdt = f * xhd
yhdt = f * ( yhd + astrom(27) )
zhdt = f * zhd
! Cartesian -HA,Dec to Cartesian Az,El (S=0,E=90).
xaet = astrom(25)*xhdt - astrom(26)*zhdt
yaet = yhdt
zaet = astrom(26)*xhdt + astrom(25)*zhdt
! Azimuth (N=0,E=90).
if ( xaet/=0.0_wp .or. yaet/=0.0_wp ) then
azobs = atan2 ( yaet, -xaet )
else
azobs = 0.0_wp
end if
! ----------
! Refraction
! ----------
! Cosine and sine of altitude, with precautions.
r = max ( sqrt ( xaet*xaet + yaet*yaet ), celmin)
z = max ( zaet, selmin )
! A*tan(z)+B*tan^3(z) model, with Newton-Raphson correction.
tz = r/z
w = astrom(30)*tz*tz
del = ( astrom(29) + w ) * tz / &
( 1.0_wp + ( astrom(29) + 3.0_wp*w ) / ( z*z ) )
! Apply the change, giving observed vector.
cosdel = 1.0_wp - del*del/2.0_wp
f = cosdel - del*z/r
xaeo = xaet*f
yaeo = yaet*f
zaeo = cosdel*zaet + del*r
! Observed ZD.
zdobs = atan2 ( sqrt ( xaeo*xaeo + yaeo*yaeo ), zaeo )
! Az/El vector to HA,Dec vector (both right-handed).
v(1) = astrom(25)*xaeo + astrom(26)*zaeo
v(2) = yaeo
v(3) = - astrom(26)*xaeo + astrom(25)*zaeo
! To spherical -HA,Dec.
call C2S ( v, hmobs, dcobs )
! Right ascension (with respect to CIO).
raobs = astrom(28) + hmobs
! Return the results.
aob = ANP(azobs)
zob = zdobs
hob = -hmobs
dob = dcobs
rob = ANP(raobs)
end subroutine ATIOQ