Project Celestial Coordinates to a Sky Dome in GLSL

Projecting celestial coordinates onto a 3D hemispherical sky dome requires transforming spherical angles into Cartesian unit vectors. This article outlines the mathematical formulas used to convert celestial coordinates—both horizontal (Altitude/Azimuth) and equatorial (Right Ascension/Declination)—into local 3D coordinates, followed by optimized GLSL shader implementations for real-time rendering.

Coordinate Systems Overview

Rendering a celestial sky dome typically involves two primary systems:

1. Equatorial Coordinates: Right Ascension (\(\alpha\) or \(\text{RA}\)) and Declination (\(\delta\) or \(\text{Dec}\)). These represent fixed positions on the celestial sphere.
2. Horizontal Coordinates: Azimuth (\(A\)) and Altitude (\(a\) or Elevation). These represent the observer's local view, where Altitude ranges from \(0^\circ\) to \(90^\circ\) (horizon to zenith) and Azimuth ranges from \(0^\circ\) to \(360^\circ\).

In graphics pipelines, the sky dome is typically modeled as a unit hemisphere where the vertical axis represents the observer's zenith (commonly the \(+Y\) axis).

Direct Horizontal Projection Formula

If coordinates are already in Altitude (\(a\)) and Azimuth (\(A\)), the projection to a Cartesian unit vector \(\vec{P} = (x, y, z)\) on a hemisphere where \(+Y\) is Up, \(+Z\) is North, and \(+X\) is East is:

\[x = \cos(a) \sin(A)\] \[y = \sin(a)\] \[z = \cos(a) \cos(A)\]

For a standard OpenGL camera space where North is \(-Z\):

\[x = \cos(a) \sin(A)\] \[y = \sin(a)\] \[z = -\cos(a) \cos(A)\]

Converting Equatorial (RA/Dec) to Horizontal (Alt/Az)

To place stars defined by Right Ascension (\(\alpha\)) and Declination (\(\delta\)) onto a local sky dome, transform them using the observer's geographic latitude (\(\phi\)) and Local Sidereal Time (\(\text{LST}\)).

1. Calculate the Local Hour Angle (\(H\))

\[H = \text{LST} - \alpha\]

2. Calculate Altitude (\(a\))

\[\sin(a) = \sin(\phi) \sin(\delta) + \cos(\phi) \cos(\delta) \cos(H)\] \[a = \arcsin(\sin(a))\]

3. Calculate Azimuth (\(A\))

\[\cos(A) = \frac{\sin(\delta) - \sin(\phi) \sin(a)}{\cos(\phi) \cos(a)}\]

Using the two-argument arctangent for quadrant correctness:

\[A = \text{atan2}(-\sin(H) \cos(\delta), \cos(\phi) \sin(\delta) - \sin(\phi) \cos(\delta) \cos(H))\]

GLSL Implementation

The following GLSL functions demonstrate projecting both coordinate types onto a hemispherical unit dome.

Projecting Altitude and Azimuth

// Converts Altitude and Azimuth (in radians) to a unit vector on the hemisphere
// +Y = Zenith, +X = East, -Z = North
vec3 horizontalToDome(float altitude, float azimuth) {
    float cosAlt = cos(altitude);
    return vec3(
        cosAlt * sin(azimuth),
        sin(altitude),
        -cosAlt * cos(azimuth)
    );
}

Full Equatorial (RA/Dec) Projection Matrix

Instead of computing trigonometric conversions per vertex or fragment, you can project equatorial coordinates into a rotating equatorial Cartesian frame and multiply by an observer orientation matrix.

// Converts RA and Dec directly to a celestial unit vector (+Y = Celestial North Pole)
vec3 celestialToCartesian(float ra, float dec) {
    float cosDec = cos(dec);
    return vec3(
        cosDec * cos(ra),
        sin(dec),
        cosDec * sin(ra)
    );
}

// Transform equatorial vector to local horizontal sky dome
vec3 equatorialToLocalDome(vec3 eqDir, float lat, float lst) {
    // Rotate by Local Sidereal Time around the celestial pole (Y axis)
    float cosLST = cos(lst);
    float sinLST = sin(lst);
    mat3 rotLST = mat3(
        cosLST,  0.0, sinLST,
        0.0,     1.0, 0.0,
        -sinLST, 0.0, cosLST
    );

    vec3 hourAngleDir = rotLST * eqDir;

    // Rotate by observer colatitude (90 deg - latitude) around X axis to align Zenith to +Y
    float angle = 1.57079632679 - lat; // PI/2 - latitude
    float cosLat = cos(angle);
    float sinLat = sin(angle);
    mat3 rotLat = mat3(
        1.0, 0.0,     0.0,
        0.0, cosLat, -sinLat,
        0.0, sinLat,  cosLat
    );

    return rotLat * hourAngleDir;
}

Mapping Dome Direction to 2D Planar Dome Textures

To project a hemisphere direction vector vec3 dir (where dir.y >= 0.0) onto a 2D circular fish-eye texture map:

vec2 domeToUV(vec3 dir) {
    // Equidistant fisheye projection
    float r = atan(length(dir.xz), dir.y) / 1.57079632679; // Normalized radius [0, 1]
    vec2 normXZ = normalize(dir.xz);
    return normXZ * (r * 0.5) + vec2(0.5);
}