| home | send comment | send link | add bookmark |

Coordinate conversion: celestial to horizon

by Stephen R. Schmitt


Equatorial coordinates

By extending the lines of latitude and longitude outward from the Earth and onto the inside of the celestial sphere we get the equatorial coordinate system. The coordinates of stars, planets, and other celestial objects corresponding to latitude and longitude are declination (DEC) and right ascension (RA).

The declination of an object is its angle in degrees, minutes, and seconds of arc above or below the celestial equator. The right ascension is the angle between an object and the location of the vernal equinox (First Point in Aries) measured eastward along the celestial equator in hours, minutes, and seconds of sidereal time. Since the location of the vernal equinox changes due to the precession of the Earth's axis of rotation, coordinates must be given with reference to a date or epoch.

Right ascension is given in time units. One hour corresponds to 1/24 of a circle, or 15° of arc. As the Earth rotates, the sky moves to the West by about 1 hour of right ascension during each hour of clock time or exactly one hour of sidereal time. The Earth makes one full revolution in about 23 hours and 56 minutes of clock time or 24 hours of sidereal time. Sidereal time corresponds to the right ascension of the zenith, the point in the sky directly overhead.

For example, the coordinates of the star Regulus (Leo α) for epoch J2000 are:

RA:   10h 08m 22.3s
DEC: +11° 58' 02"
When the local sidereal time is 10h 08m 22.3s, it would be on the local meridian.

Horizon coordinates: azimuth and altitude

This is a local coordinate system to use for locating objects in the night sky as seen from a point on the Earth's surface. Azimuth is the angle of a celestial object around the sky from north. It is measure along the horizon in from North 0° through East 90°, South 180°, West 270° and back to North. Altitude is the complement of the zenith angle, which is the angle from the local meridian to the hour circle of object being observed. An object directly overhead would have an altitude of 90°. An object with a calculated altitude of 0° may not appear exactly on the horizon due to the refraction of light through the atmosphere.

Coordinate transformation

The azimuth (AZ) and altitude (ALT) of an object in the sky can be calculated easily using the date, universal time (UT), and the latitude (LAT) and longitude (LON) of the observing site and the right ascension (RA) and declination (DEC) of the object. All coordinates are expressed in degrees in the range 0° to 360°, so that trigonometric functions can be used for coordinate conversion.

Local Mean Sidereal Time

The mean sidereal time (MST) is calculated from a polynomial function of UT since epoch J2000. This formula gives MST, the sidereal time at the Greenwich meridian (at longitude 0°) in degrees. To get local mean sidereal time (LMST), add longitude if East or subtract longitude if West.

MST = f(UT)
LMST = MST + LON 

Hour Angle

The hour angle is the angle between an observer's meridian projected onto the celestial sphere and the right ascension of a celestial body. It is used in coordinate conversion.

HA = LMST - RA 

HA and DEC to ALT and AZ

Using the RA, DEC and HA for the object, and the Latitude (LAT) of the observing site, the following formulas give the ALT and AZ of the object at the time and longitude that was used to calculate HA.

sin(ALT) = sin(DEC)·sin(LAT) + cos(DEC)·cos(LAT)·cos(HA)

           sin(DEC) - sin(ALT)·sin(LAT)
cos(A)   = ----------------------------
                cos(ALT)·cos(LAT)
If sin(HA) is negative, then AZ = A, otherwise AZ = 360 - A

Zeno source code

Zeno 1.2 is an interpreter for the Zeno programming language. It is an easy to learn and is suitable for educational purposes.


type horizon   : record
                 alt : real
                 az  : real
                 end record

const PI : real := 2*arccos(0)

program

    var tm, ra, dec, lat, lon : real
    var h : horizon

    % get UTC from MS Windows operating system
    % tm := systemtime

    % or select UTC
    tm := datetime( 2004, 4, 7, 1, 0, 0 )

    % Pleiades M45, 03h47.0m, +24°07'
    put "Pleiades (M45) "...
    ra := ra2real( 3, 47.0 )
    dec := dms2real( 24, 07, 0 )
    
    % Boston, MA US  42°21' N, 71°04' W
    put "from Boston, MA"
    lat := dms2real( 42, 21, 0 )
    lon := dms2real( -71, 4, 0 )
    
    coord_to_horizon( tm, ra, dec, lat, lon, h )
    put "alt = ", h.alt, "°, az = ", h.az, "°"
    put "UTC: ", datetimestr( tm )

end program

%
% compute horizon coordinates from ra, dec, lat, lon, and utc
% ra, dec, lat, lon in  degrees
% utc is a time number in seconds
%
% results returned in h : horizon record structure
%
procedure coord_to_horizon( utc, ra, dec, lat, lon : real, var h : horizon )

    var lmst, ha, sin_alt, cos_az, alt, az : real

    % compute hour angle in degrees
    ha := mean_sidereal_time( utc, lon ) - ra
    if (ha < 0) then
        ha := ha + 360
    end if

    % convert degrees to radians
    ha  := ha*PI/180
    dec := dec*PI/180
    lat := lat*PI/180

    % compute altitude in radians
    sin_alt := sin(dec)*sin(lat) + cos(dec)*cos(lat)*cos(ha)
    alt := arcsin(sin_alt)
    
    % compute azimuth in radians
    % divide by zero error at poles or if alt = 90 deg
    cos_az := (sin(dec) - sin(alt)*sin(lat))/(cos(alt)*cos(lat))
    az  := arccos(cos_az)

    % convert radians to degrees
    h.alt := alt*180/PI
    h.az  := az*180/PI

    % choose hemisphere
    if (sin(ha) > 0) then
        h.az  := 360 - h.az
    end if

end procedure

%
% "mean_sidereal_time" returns the Mean Sidereal Time in units of degrees. 
% Use lon := 0 to get the Greenwich MST. 
% East longitudes are positive; West longitudes are negative
%
% returns: time in degrees
%
function mean_sidereal_time( now, lon : real ) : real

    var year   : int := dateyear( now )
    var month  : int := datemonth( now )
    var day    : int := dateday( now )
    var hour   : int := timehour( now )
    var minute : int := timeminute( now )
    var second : int := timesecond( now )

    if (month = 1) or (month = 2) then
        year  := year - 1
        month := month + 12
    end if

    var a, b, c, d : int

    a := floor( year/100 )
    b := 2 - a + floor( a/4 )
    c := floor( 365.25*year )
    d := floor( 30.6001*( month + 1 ) )

    var jd, jt, mst : real

    % days since J2000.0
    jd := b + c + d - 730550.5 + day 
        + (hour + minute/60.0 + second/3600.0)/24.0
    
    % julian centuries since J2000.0
    jt := jd/36525.0

    % mean sidereal time
    mst := 280.46061837 + 360.98564736629*jd 
         + 0.000387933*jt*jt - jt*jt*jt/38710000 + lon

    if (mst > 0.0) then
        while (mst > 360.0)
            mst := mst - 360.0
        end while
    else
        while (mst < 0.0)
            mst := mst + 360.0
        end while
    end if
        
    return mst

end function

% convert right ascension (hours, minutes) to degrees as real
function ra2real( hr : int, min : real ) : real

    assert (hr >= 0) and ( hr <= 24 )
    assert (min >= 0) and (min < 60)
    
    return 15*(hr + min/60)

end function

% convert angle (deg, min, sec) to degrees as real
function dms2real( deg, min, sec : int ) : real

    var rv : real

    assert (min >= 0) and (min < 60)
    assert (sec >= 0) and (sec < 60)
   
    if deg < 0 then
        rv := deg - min/60 - sec/3600
    else
        rv := deg + min/60 + sec/3600
    end if

    return rv

end function


% generate date time string from time value
function datetimestr( tm : real ) : string

    var str : string := ""
    var yy, mm, dd, hr, mi, se : int
    yy := dateyear( tm )
    mm := datemonth( tm )
    dd := dateday( tm )
    hr := timehour( tm )
    mi := timeminute( tm )
    se := timesecond( tm )    
    str := d2(mm) & "/" & d2(dd) & "/" & intstr(yy,4) & " " 
         & d2(hr) & ":" & d2(mi) & ":" & d2(se)
    return str

end function

% format two digits with leading zero if needed
function d2( n : int ) : string

    assert (0 <= n) and (n < 100)
    if n < 10 then
        return "0" & intstr( n, 1 )
    else
        return intstr( n, 2 )
    end if

end function

Sample output

Pleiades (M45) from Boston, MA
alt = 21.0656°, az = 283.967°
UTC: 04/07/2004 01:00:00


Copyright © 2005, Stephen R. Schmitt