LAL  7.1.7.1-9b0afb2

Detailed Description

Converts among equatorial, geographic, and horizon coordinates.

Author
Creighton, T. D.

The functions LALEquatorialToGeographic() and LALGeographicToEquatorial() convert between equatorial and geographic coordinate systems, reading coordinates in the first system from *input and storing the new coordinates in *output. The two pointers may point to the same object, in which case the conversion is done in place. The functions will also check input->system and set output->system as appropriate. Because the geographic coordinate system is not fixed, one must also specify the time of the transformation in *gpsTime.

The function LALSystemToHorizon() transforms coordinates from either celestial equatorial coordinates or geographic coordinates to a horizon coordinate system, reading coordinates in the first system from *input and storing the horizon coordinates in *output, as above. The parameter *zenith specifies the direction of the vertical axis in the original coordinate system; the routine checks to see that input->system and zenith->system agree. Normally this routine is used to convert from geographic latitude and longitude to a horizon system, in which case *zenith simply stores the geographic (geodetic) coordinates of the observer; if converting from equatorial coordinates, zenith->longitude should store the local mean sidereal time of the horizon system.

The function LALHorizonToSystem() does the reverse of the above, transforming coordinates from horizon coordinates to either equatorial or geographic coordinates as specified by zenith->system; the value of output->system is set to agree with zenith->system.

Although it is conventional to specify an observation location by its geodetic coordinates, some routines may provide or require geocentric coordinates. The routines LALGeocentricToGeodetic() and LALGeodeticToGeocentric() perform this computation, reading and writing to the variable parameter structure *location. The function LALGeocentricToGeodetic() reads the fields location->x, y, z, and computes location->zenith and location->altitude. The function LALGeodeticToGeocentric() does the reverse, and also sets the fields location->position and location->radius.

Algorithm

These routines follow the formulae in Sec. 5.1 of [17] , which we reproduce below.

Geographic coordinates:
Since geographic and equatorial coordinates share the same \(z\)-axis, the geographic latitude \(\phi\) of a direction in space is the same as its declination \(\delta\), and longitude \(\lambda\) and right ascension \(\alpha\) differ only through the rotation of the Earth:

\begin{equation} \label{eq_lambda_geographic} \lambda = \alpha - \left(\frac{2\pi\,\mathrm{radians}} {24\,\mathrm{hours}}\right)\times\mathrm{GMST} \; , \end{equation}

where GMST is Greenwich mean sidereal time. The conversion routines here simply use the functions in the date package to compute GMST for a given GPS time, and add it to the longitude. While this is simple enough, it does involve several function calls, so it is convenient to collect these into one routine.
Horizon coordinates:
We correct a typographical error on the second line of Eq. 5.45 of [17] , (it should have \(\cos A\), not \(\sin A\)). We also note that while our latitudinal coordinate is just the altitude \(a\) in this system, our longitudinal coordinate increases counterclockwise, and thus corresponds to the negative of the azimuth \(A\) as defined by [17] . So we have:

\begin{eqnarray} \label{eq_altitude_horizon} a & = & \arcsin(\sin\delta\sin\phi + \cos\delta\cos\phi\cos h) \; , \\ \label{eq_azimuth_horizon} -A & = & \arctan\!2(\cos\delta\sin h, \sin\delta\cos\phi - \cos\delta\sin\phi\cos h) \; , \end{eqnarray}

where \(\delta\) is the declination (geographic latitude) of the direction being transformed, \(\phi\) is the geographic latitude of the observer's zenith (i.e. the observer's geodetic latitude), and \(h\) is the hour angle of the direction being transformed. This is defined as:

\begin{eqnarray} h & = & \lambda_\mathrm{zenith} - \lambda\\ & = & \mathrm{LMST} - \alpha \end{eqnarray}

where LMST is the local mean sidereal time at the point of observation. The inverse transformation is:

\begin{eqnarray} \label{eq_delta_horizon} \delta & = & \arcsin(\sin a\sin\phi + \cos a\cos A\cos\phi) \; , \\ \label{eq_h_horizon} h & = & \arctan\!2[\cos a\sin(-A), \sin a\cos\phi - \cos a\cos A\sin\phi] \; . \end{eqnarray}

As explained in Module CelestialCoordinates.c, the function \(\arctan\!2(y,x)\) returns the argument of the complex number \(x+iy\).

inject_geodetic.png
The difference between geodetic and geocentric latitude.
Geocentric coordinates:
As shown in this figure, the ellipticity of the Earth means that the vertical axis of a point on the Earth's surface does not pass through the geometric centre of the Earth. This means that the geodetic latitude of a location (defined as the latitude angle \(\phi_\mathrm{geodetic}\) of that location's zenith direction) is typically some 10 arcminutes larger than its geocentric latitude (defined as the latitude angle \(\phi_\mathrm{geographic}\) of the position vector from the geocentre through the location). Cartographers traditionally refer to locations by their geodetic coordinates, since these can be determined locally; however, geocentric coordinates are required if one wants to construct a uniform Cartesian system for the Earth as a whole.

To transform from geodetic to geocentric coordinates, one first defines a "reference ellipsoid", the best-fit ellipsoid to the surface of the Earth. This is specified by the polar and equatorial radii of the Earth \(r_p\) and \(r_e\), or equivalently by \(r_e\) and a flattening factor:

\[ f \equiv 1 - \frac{r_p}{r_e} = 0.00335281 \; . \]

(This constant will eventually migrate into Header LALConstants.h.) The surface of the ellipsoid is then specified by the equation

\[ r = r_e ( 1 - f\sin^2\phi ) \; , \]

where \(\phi=\phi_\mathrm{geodetic}\) is the geodetic latitude. For points off of the reference ellipsoid, the transformation from geodetic coordinates \(\lambda\), \(\phi\) to geocentric Cartesian coordinates is:

\begin{eqnarray} x & = & ( r_e C + h ) \cos\phi\cos\lambda \; , \\ y & = & ( r_e C + h ) \cos\phi\sin\lambda \; , \\ z & = & ( r_e S + h ) \sin\phi \; , \end{eqnarray}

where

\begin{eqnarray} C & = & \frac{1}{\sqrt{\cos^2\phi + (1-f)^2\sin^2\phi}} \; , \\ S & = & (1-f)^2 C \; , \end{eqnarray}

and \(h\) is the perpendicular elevation of the location above the reference ellipsoid. The geocentric spherical coordinates are given simply by:

\begin{eqnarray} r & = & \sqrt{ x^2 + y^2 + z^2 } \; , \\ \lambda_\mathrm{geocentric} & = & \lambda \quad = \quad \lambda_\mathrm{geodetic} \; , \\ \phi_\mathrm{geocentric} & = & \arcsin( z/r ) \; . \end{eqnarray}

When computing \(r\) we are careful to factor out the largest component before computing the sum of squares, to avoid floating-point overflow; however this should be unnecessary for radii near the surface of the Earth.

The inverse transformation is somewhat trickier. Eq. 5.29 of [17] conveniently gives the transformation in terms of a sequence of intermediate variables, but unfortunately these variables are not particularly computer-friendly, in that they are prone to underflow or overflow errors. The following equations essentially reproduce this sequence using better-behaved methods of calculation.

Given geocentric Cartesian coordinates \(x=r\cos\phi_\mathrm{geocentric}\cos\lambda\), \(y=r\cos\phi_\mathrm{geocentric}\sin\lambda\), and \(z=r\sin\phi_\mathrm{geocentric}\), one computes the following:

\begin{eqnarray} \varpi & = & \sqrt{ \left(\frac{x}{r_e}\right)^2 + \left(\frac{y}{r_e}\right)^2 } \;,\\ E & = & (1-f)\left|\frac{z}{r_e}\right| - f(2-f) \;,\\ F & = & (1-f)\left|\frac{z}{r_e}\right| + f(2-f) \;,\\ P & = & \frac{4}{3}\left( EF + \varpi^2 \right) \quad = \quad \frac{4}{3}\left[ \varpi^2 + (1-f)^2\left(\frac{z}{r_e}\right)^2 - f^2(2-f)^2 \right] \;,\\ Q & = & 2\varpi(F^2 - E^2) \quad = \quad 8\varpi f(1-f)(2-f)\left|\frac{z}{r_e}\right| \;,\\ D & = & P^3 + Q^2 \;,\\ v & = & \left\{\begin{array}{lr} \left(\sqrt{D}+Q\right)^{1/3} - \left(\sqrt{D}-Q\right)^{1/3} & D\geq0 \\ 2\sqrt{-P}\cos\left(\frac{1}{3} \arccos\left[\frac{Q}{-P\sqrt{-P}}\right]\right) & D\leq0 \end{array}\right.\\ W & = & \sqrt{E^2 + \varpi v}\\ G & = & \frac{1}{2}\left(E+W\right)\;,\\ t & = & \sqrt{G^2+\frac{\varpi^2 F - \varpi vG}{W}}-G \;. \end{eqnarray}

Once we have \(t\) and \(\varpi\), we can compute the geodetic longitude \(\lambda\), latitude \(\phi\), and elevation \(h\):

\begin{eqnarray} \lambda & = & \arctan\!2(y,x) \; , \\ \phi & = & \mathrm{sgn}({z})\arctan\left[\frac{1}{2(1-f)} \left(\frac{(\varpi-t)(\varpi+t)}{\varpi t}\right)\right] \; , \\ h & = & r_e(\varpi-t/\varpi)\cos\phi + [z-\mathrm{sgn}({z})r_e(1-f)]\sin\phi \; . \end{eqnarray}

inject_geodeticsing.png
Singular surfaces in the geodetic coordinate system. The ellipticity of this spheroid has been exaggerated compared with the Earth.

These formulae still leave certain areas where coordinate singularities or numerical cancelations can occur. Some of these have been dealt with in the code:

Ellipsoidal vs. orthometric elevation:
In this module it is assumed that all elevations refer heights above the reference ellipsoid. This is the elevation computed by such techniques as GPS triangulation. However, the "true" orthometric elevation refers to the height above the mean sea level or geoid, a level surface in the Earth's gravitational potential. Thus, even if two points have the same ellipsoidal elevation, water will still flow from the higher to the lower orthometric elevation.

The difference between the geoid and reference ellipsoid is called the "undulation of the geoid", and can vary by over a hundred metres over the Earth's surface. However, it can only be determined through painstaking measurements of local variations in the Earth's gravitational field. For this reason we will ignore the undulation of the geoid.

Uses

Prototypes

void LALEquatorialToGeographic (LALStatus *stat, SkyPosition *output, SkyPosition *input, LIGOTimeGPS *gpsTime)
 
void LALGeographicToEquatorial (LALStatus *stat, SkyPosition *output, SkyPosition *input, LIGOTimeGPS *gpsTime)
 
void LALSystemToHorizon (LALStatus *stat, SkyPosition *output, SkyPosition *input, const SkyPosition *zenith)
 
void LALHorizonToSystem (LALStatus *stat, SkyPosition *output, SkyPosition *input, const SkyPosition *zenith)
 
void LALGeodeticToGeocentric (LALStatus *stat, EarthPosition *location)
 
void LALGeocentricToGeodetic (LALStatus *stat, EarthPosition *location)
 

Function Documentation

◆ LALEquatorialToGeographic()

void LALEquatorialToGeographic ( LALStatus stat,
SkyPosition output,
SkyPosition input,
LIGOTimeGPS gpsTime 
)
See also
See documentation in Module TerrestrialCoordinates.c

Definition at line 345 of file TerrestrialCoordinates.c.

◆ LALGeographicToEquatorial()

void LALGeographicToEquatorial ( LALStatus stat,
SkyPosition output,
SkyPosition input,
LIGOTimeGPS gpsTime 
)
See also
See documentation in Module TerrestrialCoordinates.c

Definition at line 379 of file TerrestrialCoordinates.c.

◆ LALSystemToHorizon()

void LALSystemToHorizon ( LALStatus stat,
SkyPosition output,
SkyPosition input,
const SkyPosition zenith 
)
See also
See documentation in Module TerrestrialCoordinates.c

Definition at line 413 of file TerrestrialCoordinates.c.

◆ LALHorizonToSystem()

void LALHorizonToSystem ( LALStatus stat,
SkyPosition output,
SkyPosition input,
const SkyPosition zenith 
)
See also
See documentation in Module TerrestrialCoordinates.c

Definition at line 468 of file TerrestrialCoordinates.c.

◆ LALGeodeticToGeocentric()

void LALGeodeticToGeocentric ( LALStatus stat,
EarthPosition location 
)
See also
See documentation in Module TerrestrialCoordinates.c

Definition at line 522 of file TerrestrialCoordinates.c.

◆ LALGeocentricToGeodetic()

void LALGeocentricToGeodetic ( LALStatus stat,
EarthPosition location 
)
See also
See documentation in Module TerrestrialCoordinates.c

Definition at line 577 of file TerrestrialCoordinates.c.