AE2HD Subroutine

public subroutine AE2HD(az, el, phi, ha, dec)

Horizon to equatorial coordinates: transform azimuth and altitude to hour angle and declination.

Status: support routine.

Notes

1. All the arguments are angles in radians.

2. The sign convention for azimuth is north zero, east +pi/2.

3. HA is returned in the range +/-pi. Declination is returned in the range +/-pi/2.

4. The latitude PHI is pi/2 minus the angle between the Earth's rotation axis and the adopted zenith. In many applications it will be sufficient to use the published geodetic latitude of the site. In very precise (sub-arcsecond) applications, PHI can be corrected for polar motion.

5. The azimuth AZ must be with respect to the rotational north pole, as opposed to the ITRS pole, and an azimuth with respect to north on a map of the Earth's surface will need to be adjusted for polar motion if sub-arcsecond accuracy is required.

6. Should the user wish to work with respect to the astronomical zenith rather than the geodetic zenith, PHI will need to be adjusted for deflection of the vertical (often tens of arcseconds), and the zero point of HA will also be affected.

7. The transformation is the same as Ve = Ry(phi-pi/2)Rz(pi)Vh, where Ve and Vh are lefthanded unit vectors in the (ha,dec) and (az,el) systems respectively and Rz and Ry are rotations about first the z-axis and then the y-axis. (n.b. Rz(pi) simply reverses the signs of the x and y components.) For efficiency, the algorithm is written out rather than calling other utility functions. For applications that require even greater efficiency, additional savings are possible if constant terms such as functions of latitude are computed once and for all.

8. Again for efficiency, no range checking of arguments is carried out.

Last revision: 2018 January 2

Arguments

    subroutine AE2HD (az, el, phi, ha, dec)

    implicit none

    real(wp),intent(in) :: az !! azimuth
    real(wp),intent(in) :: el !! elevation
    real(wp),intent(in) :: phi !! observatory latitude
    real(wp),intent(out) :: ha !! hour angle
    real(wp),intent(out) :: dec !! declination

    real(wp) :: sa, ca, se, ce, sp, cp, x, y, z, r

    !  Useful trig functions.
    sa = sin(az)
    ca = cos(az)
    se = sin(el)
    ce = cos(el)
    sp = sin(phi)
    cp = cos(phi)

    !  Az,Alt unit vector.
    x = - ca*ce*sp + se*cp
    y = - sa*ce
    z = ca*ce*cp + se*sp

    !  To spherical.
    r = sqrt(x*x + y*y)
    if ( r==0.0_wp ) then
       ha = 0.0_wp
    else
       ha = atan2(y,x)
    end if
    dec = atan2(z,r)

    end subroutine AE2HD