TerrestrialCoordinates.c File Reference

Author:: Creighton, T.Converts among equatorial, geographic, and horizon coordinates.

#include <math.h>
#include <lal/LALStdlib.h>
#include <lal/LALConstants.h>
#include <lal/Date.h>
#include <lal/Sort.h>
#include <lal/SkyCoordinates.h>
#include <string.h>

Include dependency graph for TerrestrialCoordinates.c:

Go to the source code of this file.

Defines

#define LAL_EARTHFLAT   (0.00335281)

#define LAL_HSERIES   (0.0001)

#define LAL_BSERIES   (0.001)

Functions

NRCSID (TERRESTRIALCOORDINATESC,"$Id: TerrestrialCoordinates.c,v 1.20 2008/05/22 19:27:42 teviet Exp $")

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)

Detailed Description

Author:: Creighton, T.Converts among equatorial, geographic, and horizon coordinates.

Date:

Date: 2008/05/22 19:27:42

Id: TerrestrialCoordinates.c,v 1.20 2008/05/22 19:27:42 teviet Exp

Description

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 Lang_K1999, which we reproduce below.

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.

error on the second line of Eq.(5.45) of Lang_K1999 (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 Lang_K1999. 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 \nonumber\\ & = & \mathrm{LMST} - \alpha \nonumber \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 CelestialCoordinates.c, the function $\arctan\!2(y,x)$ returns the argument of the complex number $x+iy$ .

Fig. 1: The difference between geodetic and geocentric latitude.

Geocentric coordinates: As shown in Fig. 1, 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 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 Lang_K1999 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 } \;,\nonumber\\ E & = & (1-f)\left|\frac{z}{r_e}\right| - f(2-f) \;,\nonumber\\ F & = & (1-f)\left|\frac{z}{r_e}\right| + f(2-f) \;,\nonumber\\ 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] \;,\nonumber\\ Q & = & 2(F^2 - E^2) \quad = \quad 8f(1-f)(2-f)\left|\frac{z}{r_e}\right| \;,\nonumber\\ D & = & P^3 + \varpi^2 Q^2 \;,\nonumber\\ v & = & \left\{\begin{array}{lr} \left(\sqrt{D}+\varpi Q\right)^{1/3} - \left(\sqrt{D}-\varpi Q\right)^{1/3} & D\geq0 \\ 2\sqrt{-P}\cos\left(\frac{1}{3} \arccos\left[\frac{Q\varpi}{P\sqrt{-P}}\right]\right) & D\leq0 \end{array}\right.\nonumber\\ G & = & \mbox{$\frac{1}{2}$}\left(E+\sqrt{E^2 + \varpi v}\right)\;,\nonumber\\ H & = & \frac{\varpi^2 F - \varpi vG}{G^2(2G-E)} \;,\nonumber\\ t & = & G\left(\sqrt{1+H}-1\right) \;.\nonumber \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{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}$

These formulae, however, introduce certain concerns of numerical precision that have been only partially dealt with in this code. Specifically:

There is a coordinate singularity at $\varpi=0$ , which we deal with by setting $\phi=\pm90^\circ$ and $\lambda$ arbitrarily to $0^\circ$ . When as well, we arbitrarily choose the positive sign for $\phi$ . However, the computation of in particular has tricky cancelations as $\varpi\rightarrow0$ , which may give rise to numerical errors. These have not yet been thoroughly explored.

There is another coordinate singularity when $D\rightarrow0$ , which defines an ellipsoid with equatorial radius km and axial height km. Within this ellipsoid, lines of constant latitude begin to cross one another. The listed solution is an analytic continuation of the exterior solution which assigns these points a unique, if arbitrary, geodetic latitude. This solution has some peculiar behaviour, such as giving points in the equatorial plane a positive latitude. In practice, however, users will rarely be interested coordinate transformations deep within the Earth's core.

The equations for and have square and cube roots of expressions involving squares and cubes of numbers. For formal robustness one should factor out the leading-order dependence, so that one is assured of taking squares and cubes of numbers near unity. However, we are using REAL8 precision, and have already normalized our distances by the Earth's radius , so the point is almost certainly irrelevant.

The expression for may go to zero, leading to precision errors in ; the number of digits of precision lost is on the order of the number of leading zeros after the decimal place in . I arbitrarily state that we should not lose more than 4 of our 16 decimal places of precision, meaning that we should series-expand the square root for $H<10^{-4}$ . To get our 12 places of precision back, we need an expansion to :
$t \approx G\left(\frac{1}{2}H - \frac{3}{8}H^2 + \frac{5}{16}H^3\right) \; .$

When computing $\phi$ , we first compute $t-\varpi$ , $t+\varpi$ , $t^{-1}$ , and $\varpi^{-1}$ , sort them by order of magnitude, and alternately multiply large and small terms. We note that if the argument of the $\arctan$ function is large we have

$\arctan(x) = \mathrm{sgn}(x)\frac{\pi}{2} - \arctan\left(\frac{1}{x}\right) \; ,$

but the atan() function in the C math library should be smart enough to do this itself.

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

LALGPStoUTC()
XLALGreenwichMeanSiderealTime()
LALDHeapSort()

Definition in file TerrestrialCoordinates.c.

Define Documentation

#define LAL_EARTHFLAT (0.00335281)

Definition at line 677 of file TerrestrialCoordinates.c.

#define LAL_HSERIES (0.0001)

Definition at line 678 of file TerrestrialCoordinates.c.

#define LAL_BSERIES (0.001)

Definition at line 679 of file TerrestrialCoordinates.c.

Function Documentation

NRCSID	(	TERRESTRIALCOORDINATESC	,
		"$Id: TerrestrialCoordinates.	c,
		v 1.20 2008/05/22 19:27:42 teviet Exp $"
	)

void LALEquatorialToGeographic	(	LALStatus *	stat,
		SkyPosition *	output,
		SkyPosition *	input,
		LIGOTimeGPS *	gpsTime
	)

Definition at line 685 of file TerrestrialCoordinates.c.

void LALGeographicToEquatorial	(	LALStatus *	stat,
		SkyPosition *	output,
		SkyPosition *	input,
		LIGOTimeGPS *	gpsTime
	)

Definition at line 719 of file TerrestrialCoordinates.c.

void LALSystemToHorizon	(	LALStatus *	stat,
		SkyPosition *	output,
		SkyPosition *	input,
		const SkyPosition *	zenith
	)

Definition at line 753 of file TerrestrialCoordinates.c.

void LALHorizonToSystem	(	LALStatus *	stat,
		SkyPosition *	output,
		SkyPosition *	input,
		const SkyPosition *	zenith
	)

Definition at line 808 of file TerrestrialCoordinates.c.

void LALGeodeticToGeocentric	(	LALStatus *	stat,
		EarthPosition *	location
	)

Definition at line 862 of file TerrestrialCoordinates.c.

void LALGeocentricToGeodetic	(	LALStatus *	stat,
		EarthPosition *	location
	)

Definition at line 917 of file TerrestrialCoordinates.c.

Generated on Sat Aug 30 03:14:16 2008 for LAL by

1.5.2


Defines
#define	LAL_EARTHFLAT (0.00335281)
#define	LAL_HSERIES (0.0001)
#define	LAL_BSERIES (0.001)
Functions
	NRCSID (TERRESTRIALCOORDINATESC,"$Id: TerrestrialCoordinates.c,v 1.20 2008/05/22 19:27:42 teviet Exp $")
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)