LAL  7.5.0.1-8083555
TerrestrialCoordinates.c
Go to the documentation of this file.
1 /*
2 * Copyright (C) 2007 Jolien Creighton, Reinhard Prix, Teviet Creighton, John Whelan
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with with program; see the file COPYING. If not, write to the
16 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17 * MA 02110-1301 USA
18 */
19 
20 
21 #include <math.h>
22 #include <lal/LALStdlib.h>
23 #include <lal/LALConstants.h>
24 #include <lal/Date.h>
25 #include <lal/Sort.h>
26 #include <lal/SkyCoordinates.h>
27 #include <string.h>
28 
29 #define LAL_EARTHFLAT (0.00335281)
30 #define LAL_HSERIES (0.0001) /* value H below which we expand sqrt(1+H)-1 */
31 #define LAL_BSERIES (0.001) /* value B below which we expand v */
32 
33 /**
34  * \author Creighton, T. D.
35  * \addtogroup TerrestrialCoordinates_c
36  * \brief Converts among equatorial, geographic, and horizon coordinates.
37  *
38  * The functions <tt>LALEquatorialToGeographic()</tt> and
39  * <tt>LALGeographicToEquatorial()</tt> convert between equatorial and
40  * geographic coordinate systems, reading coordinates in the first system
41  * from <tt>*input</tt> and storing the new coordinates in <tt>*output</tt>.
42  * The two pointers may point to the same object, in which case the
43  * conversion is done in place. The functions will also check
44  * <tt>input->system</tt> and set <tt>output->system</tt> as appropriate.
45  * Because the geographic coordinate system is not fixed, one must also
46  * specify the time of the transformation in <tt>*gpsTime</tt>.
47  *
48  * The function <tt>LALSystemToHorizon()</tt> transforms coordinates from
49  * either celestial equatorial coordinates or geographic coordinates to a
50  * horizon coordinate system, reading coordinates in the first system
51  * from <tt>*input</tt> and storing the horizon coordinates in
52  * <tt>*output</tt>, as above. The parameter <tt>*zenith</tt> specifies the
53  * direction of the vertical axis <em>in the original coordinate
54  * system</em>; the routine checks to see that <tt>input->system</tt> and
55  * <tt>zenith->system</tt> agree. Normally this routine is used to convert
56  * from \e geographic latitude and longitude to a horizon system, in
57  * which case <tt>*zenith</tt> simply stores the geographic (geodetic)
58  * coordinates of the observer; if converting from equatorial
59  * coordinates, <tt>zenith->longitude</tt> should store the local mean
60  * sidereal time of the horizon system.
61  *
62  * The function <tt>LALHorizonToSystem()</tt> does the reverse of the
63  * above, transforming coordinates from horizon coordinates to either
64  * equatorial or geographic coordinates as specified by
65  * <tt>zenith->system</tt>; the value of <tt>output->system</tt> is set to
66  * agree with <tt>zenith->system</tt>.
67  *
68  * Although it is conventional to specify an observation location by its
69  * \e geodetic coordinates, some routines may provide or require
70  * \e geocentric coordinates. The routines
71  * <tt>LALGeocentricToGeodetic()</tt> and <tt>LALGeodeticToGeocentric()</tt>
72  * perform this computation, reading and writing to the variable
73  * parameter structure <tt>*location</tt>. The function
74  * <tt>LALGeocentricToGeodetic()</tt> reads the fields <tt>location->x</tt>,
75  * \c y, \c z, and computes <tt>location->zenith</tt> and
76  * <tt>location->altitude</tt>. The function
77  * <tt>LALGeodeticToGeocentric()</tt> does the reverse, and also sets the
78  * fields <tt>location->position</tt> and <tt>location->radius</tt>.
79  *
80  * ### Algorithm ###
81  *
82  * These routines follow the formulae in Sec. 5.1 of \cite Lang_K1999 ,
83  * which we reproduce below.
84  *
85  * \par Geographic coordinates:
86  * Since geographic and equatorial
87  * coordinates share the same \f$z\f$-axis, the geographic latitude \f$\phi\f$ of
88  * a direction in space is the same as its declination \f$\delta\f$, and
89  * longitude \f$\lambda\f$ and right ascension \f$\alpha\f$ differ only through
90  * the rotation of the Earth:
91  * \f{equation}{
92  * \label{eq_lambda_geographic}
93  * \lambda = \alpha - \left(\frac{2\pi\,\mathrm{radians}}
94  * {24\,\mathrm{hours}}\right)\times\mathrm{GMST} \; ,
95  * \f}
96  * where GMST is Greenwich mean sidereal time. The conversion routines
97  * here simply use the functions in the date package to compute
98  * GMST for a given GPS time, and add it to the longitude. While this is
99  * simple enough, it does involve several function calls, so it is
100  * convenient to collect these into one routine.
101  *
102  * \par Horizon coordinates:
103  * We correct a typographical
104  * error on the second line of Eq. 5.45 of \cite Lang_K1999 , (it should
105  * have \f$\cos A\f$, not \f$\sin A\f$). We also note that while our latitudinal
106  * coordinate is just the altitude \f$a\f$ in this system, our longitudinal
107  * coordinate increases counterclockwise, and thus corresponds to the
108  * \e negative of the azimuth \f$A\f$ as defined by \cite Lang_K1999 .
109  * So we have:
110  * \f{eqnarray}{
111  * \label{eq_altitude_horizon}
112  * a & = & \arcsin(\sin\delta\sin\phi + \cos\delta\cos\phi\cos h) \; , \\
113  * \label{eq_azimuth_horizon}
114  * -A & = & \arctan\!2(\cos\delta\sin h, \sin\delta\cos\phi -
115  * \cos\delta\sin\phi\cos h) \; ,
116  * \f}
117  * where \f$\delta\f$ is the declination (geographic latitude) of the
118  * direction being transformed, \f$\phi\f$ is the geographic latitude of the
119  * observer's zenith (i.e.\ the observer's \e geodetic latitude), and
120  * \f$h\f$ is the <em>hour angle</em> of the direction being transformed. This
121  * is defined as:
122  * \f{eqnarray}{
123  * h & = & \lambda_\mathrm{zenith} - \lambda\\
124  * & = & \mathrm{LMST} - \alpha
125  * \f}
126  * where LMST is the local mean sidereal time at the point of
127  * observation. The inverse transformation is:
128  * \f{eqnarray}{
129  * \label{eq_delta_horizon}
130  * \delta & = & \arcsin(\sin a\sin\phi + \cos a\cos A\cos\phi) \; , \\
131  * \label{eq_h_horizon}
132  * h & = & \arctan\!2[\cos a\sin(-A), \sin a\cos\phi -
133  * \cos a\cos A\sin\phi] \; .
134  * \f}
135  * As explained in \ref CelestialCoordinates_c, the function
136  * \f$\arctan\!2(y,x)\f$ returns the argument of the complex number \f$x+iy\f$.
137  *
138  * \anchor inject_geodetic
139  * \image html inject_geodetic.png "The difference between geodetic and geocentric latitude."
140  *
141  * \par Geocentric coordinates:
142  * As shown in
143  * \ref inject_geodetic "this figure", the ellipticity of the Earth means that the
144  * vertical axis of a point on the Earth's surface does not pass through
145  * the geometric centre of the Earth. This means that the geodetic
146  * latitude of a location (defined as the latitude angle
147  * \f$\phi_\mathrm{geodetic}\f$ of that location's zenith direction) is
148  * typically some 10 arcminutes larger than its geocentric latitude
149  * (defined as the latitude angle \f$\phi_\mathrm{geographic}\f$ of the
150  * position vector from the geocentre through the location).
151  * Cartographers traditionally refer to locations by their geodetic
152  * coordinates, since these can be determined locally; however,
153  * geocentric coordinates are required if one wants to construct a
154  * uniform Cartesian system for the Earth as a whole.
155  *
156  * To transform from geodetic to geocentric coordinates, one first
157  * defines a "reference ellipsoid", the best-fit ellipsoid to the
158  * surface of the Earth. This is specified by the polar and equatorial
159  * radii of the Earth \f$r_p\f$ and \f$r_e\f$, or equivalently by \f$r_e\f$ and a
160  * flattening factor:
161  * \f[
162  * f \equiv 1 - \frac{r_p}{r_e} = 0.00335281 \; .
163  * \f]
164  * (This constant will eventually migrate into \ref LALConstants_h.)
165  * The surface of the ellipsoid is then specified by the equation
166  * \f[
167  * r = r_e ( 1 - f\sin^2\phi ) \; ,
168  * \f]
169  * where \f$\phi=\phi_\mathrm{geodetic}\f$ is the geodetic latitude. For
170  * points off of the reference ellipsoid, the transformation from
171  * geodetic coordinates \f$\lambda\f$, \f$\phi\f$ to geocentric Cartesian
172  * coordinates is:
173  * \f{eqnarray}{
174  * x & = & ( r_e C + h ) \cos\phi\cos\lambda \; , \\
175  * y & = & ( r_e C + h ) \cos\phi\sin\lambda \; , \\
176  * z & = & ( r_e S + h ) \sin\phi \; ,
177  * \f}
178  * where
179  * \f{eqnarray}{
180  * C & = & \frac{1}{\sqrt{\cos^2\phi + (1-f)^2\sin^2\phi}} \; , \\
181  * S & = & (1-f)^2 C \; ,
182  * \f}
183  * and \f$h\f$ is the perpendicular elevation of the location above the
184  * reference ellipsoid. The geocentric spherical coordinates are given
185  * simply by:
186  * \f{eqnarray}{
187  * r & = & \sqrt{ x^2 + y^2 + z^2 } \; , \\
188  * \lambda_\mathrm{geocentric} & = & \lambda \quad = \quad
189  * \lambda_\mathrm{geodetic} \; , \\
190  * \phi_\mathrm{geocentric} & = & \arcsin( z/r ) \; .
191  * \f}
192  * When computing \f$r\f$ we are careful to factor out the largest component
193  * before computing the sum of squares, to avoid floating-point overflow;
194  * however this should be unnecessary for radii near the surface of the
195  * Earth.
196  *
197  * The inverse transformation is somewhat trickier. Eq. 5.29
198  * of \cite Lang_K1999 conveniently gives the transformation in terms
199  * of a sequence of intermediate variables, but unfortunately these
200  * variables are not particularly computer-friendly, in that they are
201  * prone to underflow or overflow errors. The following equations
202  * essentially reproduce this sequence using better-behaved methods of
203  * calculation.
204  *
205  * Given geocentric Cartesian coordinates
206  * \f$x=r\cos\phi_\mathrm{geocentric}\cos\lambda\f$,
207  * \f$y=r\cos\phi_\mathrm{geocentric}\sin\lambda\f$, and
208  * \f$z=r\sin\phi_\mathrm{geocentric}\f$, one computes the following:
209  * \f{eqnarray}{
210  * \varpi & = & \sqrt{ \left(\frac{x}{r_e}\right)^2
211  * + \left(\frac{y}{r_e}\right)^2 } \;,\\
212  * E & = & (1-f)\left|\frac{z}{r_e}\right| - f(2-f) \;,\\
213  * F & = & (1-f)\left|\frac{z}{r_e}\right| + f(2-f) \;,\\
214  * P & = & \frac{4}{3}\left( EF + \varpi^2 \right) \quad = \quad
215  * \frac{4}{3}\left[ \varpi^2 + (1-f)^2\left(\frac{z}{r_e}\right)^2
216  * - f^2(2-f)^2 \right] \;,\\
217  * Q & = & 2\varpi(F^2 - E^2) \quad = \quad
218  * 8\varpi f(1-f)(2-f)\left|\frac{z}{r_e}\right| \;,\\
219  * D & = & P^3 + Q^2 \;,\\
220  * v & = & \left\{\begin{array}{lr}
221  * \left(\sqrt{D}+Q\right)^{1/3}
222  * - \left(\sqrt{D}-Q\right)^{1/3} &
223  * D\geq0 \\
224  * 2\sqrt{-P}\cos\left(\frac{1}{3}
225  * \arccos\left[\frac{Q}{-P\sqrt{-P}}\right]\right) &
226  * D\leq0 \end{array}\right.\\
227  * W & = & \sqrt{E^2 + \varpi v}\\
228  * G & = & \frac{1}{2}\left(E+W\right)\;,\\
229  * t & = & \sqrt{G^2+\frac{\varpi^2 F - \varpi vG}{W}}-G \;.
230  * \f}
231  * Once we have \f$t\f$ and \f$\varpi\f$, we can compute the geodetic longitude
232  * \f$\lambda\f$, latitude \f$\phi\f$, and elevation \f$h\f$:
233  * \f{eqnarray}{
234  * \lambda & = & \arctan\!2(y,x) \; , \\
235  * \phi & = & \mathrm{sgn}({z})\arctan\left[\frac{1}{2(1-f)}
236  * \left(\frac{(\varpi-t)(\varpi+t)}{\varpi t}\right)\right] \; , \\
237  * h & = & r_e(\varpi-t/\varpi)\cos\phi
238  * + [z-\mathrm{sgn}({z})r_e(1-f)]\sin\phi \; .
239  * \f}
240  *
241  * \anchor inject_geodeticsing
242  * \image html inject_geodeticsing.png "Singular surfaces in the geodetic coordinate system. The ellipticity of this spheroid has been exaggerated compared with the Earth."
243  *
244  * These formulae still leave certain areas where coordinate
245  * singularities or numerical cancelations can occur. Some of these have
246  * been dealt with in the code:
247  * <ul>
248  * <li> There is a coordinate singularity at \f$\varpi=0\f$, which we deal
249  * with by setting \f$\phi=\pm90^\circ\f$ and \f$\lambda\f$ arbitrarily to
250  * \f$0^\circ\f$. When \f$z=0\f$ as well, we arbitrarily choose the positive
251  * sign for \f$\phi\f$. As \f$\varpi\rightarrow0\f$, there are cancellations in
252  * the computation of \f$G\f$ and \f$t\f$, which call for special treatment
253  * (below).</li>
254  *
255  * <li> There is another coordinate singularity when \f$D\leq0\f$, where
256  * lines of constant geodetic latitude will cross each other, as shown in
257  * \ref inject_geodeticsing "this figure". That is, a given point within this
258  * region can be assigned a range of geodetic latitudes. The
259  * multi-valued region lies within an inner ellipsoid \f$P\leq0\f$, which in
260  * the case of the Earth has equatorial radius \f$r_0=r_ef(2-f)=42.6977\f$km
261  * and axial height \f$z_0=r_0/(1-f)=42.8413\f$km. The formula for \f$v\f$,
262  * above, has an analytic continuation to \f$D\leq0\f$, assigning consistent
263  * (though not unique) values of latitute and elevation to these points.</li>
264  *
265  * <li> Near the equator we have \f$Q\rightarrow0\f$, and the first
266  * expression for \f$v\f$ becomes a difference of nearly-equal numbers,
267  * leading to loss of precision. To deal with this, we write that
268  * expression for \f$v\f$ as:
269  * \f[
270  * \begin{array}{rcl@{\qquad}c@{\qquad}l}
271  * v &=& D^{1/6}\left[(1+B)^{1/3}-(1-B)^{1/3}\right]
272  * &\mbox{where}& B \;\;=\;\; \displaystyle \frac{Q}{\sqrt{D}} \\
273  * &\approx& D^{1/6}\left[\frac{2}{3}B+\frac{10}{81}B^3\right]
274  * &\mbox{as}& B \;\;\rightarrow\;\;0
275  * \end{array}
276  * \f]
277  *
278  * The switch from the "exact" formula to the series expansion is done
279  * for \f$B<10^{-3}\f$ (within about \f$2^\circ\f$ of the equator). This was
280  * found by experimentation to be the point where the inaccuracy of the
281  * series expansion is roughly equal to the imprecision in evaluating the
282  * "exact" formula. The resulting position errors are of order 1 part
283  * in \f$10^{12}\f$ or less; i.e.\ about 3--4 digits loss of precision.</li>
284  *
285  * <li> In some places we have expressions of the form \f$\sqrt{a^2+b}-a\f$,
286  * which becomes a difference of nearly equal numbers for \f$|b|\ll a^2\f$,
287  * resulting in loss of precision. There are three distinct lines or
288  * surfaces where this occurs:</li>
289  * <ol>
290  * <li> Near the ellipsoidal surface \f$P\approx0\f$, the expression \f$\sqrt{D}-Q\f$
291  * in the equation for \f$v\f$ becomes of this form.</li>
292  * <li> Near the polar axis \f$\varpi\approx0\f$, the expression for \f$t\f$
293  * becomes of this form.</li>
294  * <li> Near the polar axis \f$\varpi\approx0\f$ within the inner ellipsoid,
295  * we have \f$E<0\f$ and the expression for \f$G\f$ becomes of this form.
296  * </ol>
297  * In each case, we expand in the small parameter \f$H=b/a^2\f$, giving:
298  * \f[
299  * \sqrt{a^2+b}-a \;\;\approx\;\; a\left(\frac{1}{2} H
300  * - \frac{1}{8} H^2 + \frac{1}{16} H^3\right)
301  * \qquad\mbox{for}\qquad |H| = \left|\frac{b}{a^2}\right| \ll 1
302  * \f]
303  *
304  * We switch to the series expansion for \f$|H|<\times10^{-4}\f$, which again
305  * is the point where residual errors in the series expansion are about
306  * the same as the loss of precision without the expansion. This formula
307  * converges much better than the one for \f$v\f$ (above), resulting in
308  * negligible loss of precision in the final position.</li>
309  *
310  * <li> The only remaining known numerical singularities are at the
311  * poles of the inner ellipsoid, shown as red dots in
312  * \ref inject_geodeticsing "this figure". These points stubbornly resist
313  * high-precision calculation. However, they are extremely unlikely to
314  * come up in practice.</li>
315  * </ul>
316  *
317  * \par Ellipsoidal vs. orthometric elevation:
318  * In this module it
319  * is assumed that all elevations refer heights above the reference
320  * ellipsoid. This is the elevation computed by such techniques as GPS
321  * triangulation. However, the "true" orthometric elevation refers to
322  * the height above the mean sea level or \e geoid, a level surface
323  * in the Earth's gravitational potential. Thus, even if two points have
324  * the same ellipsoidal elevation, water will still flow from the higher
325  * to the lower orthometric elevation.
326  *
327  * The difference between the geoid and reference ellipsoid is called the
328  * "undulation of the geoid", and can vary by over a hundred metres
329  * over the Earth's surface. However, it can only be determined through
330  * painstaking measurements of local variations in the Earth's
331  * gravitational field. For this reason we will ignore the undulation of
332  * the geoid.
333  *
334  * ### Uses ###
335  *
336  * \code
337  * XLALGreenwichMeanSiderealTime()
338  * \endcode
339  *
340  */
341 /** @{ */
342 
343 /** \see See documentation in \ref TerrestrialCoordinates_c */
344 void
347  SkyPosition *input,
348  LIGOTimeGPS *gpsTime )
349 {
350  REAL8 gmst; /* siderial time (radians) */
351 
352  INITSTATUS(stat);
353 
354  /* Make sure parameter structures exist. */
355  ASSERT( input, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
356  ASSERT( output, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
357  ASSERT( gpsTime, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
358 
359  /* Make sure we're given the right coordinate system. */
360  if ( input->system != COORDINATESYSTEM_EQUATORIAL ) {
361  ABORT( stat, SKYCOORDINATESH_ESYS, SKYCOORDINATESH_MSGESYS );
362  }
363 
364  /* Compute the Greenwich mean sidereal time. */
365  gmst = fmod( XLALGreenwichMeanSiderealTime( gpsTime ), LAL_TWOPI );
366 
367  /* Add to longitude, and exit. */
369  output->longitude = fmod( input->longitude - gmst, LAL_TWOPI );
370  if ( output->longitude < 0.0 )
371  output->longitude += LAL_TWOPI;
372  output->latitude = input->latitude;
373  RETURN( stat );
374 }
375 
376 
377 /** \see See documentation in \ref TerrestrialCoordinates_c */
378 void
381  SkyPosition *input,
382  LIGOTimeGPS *gpsTime )
383 {
384  REAL8 gmst; /* siderial time (radians) */
385 
386  INITSTATUS(stat);
387 
388  /* Make sure parameter structures exist. */
389  ASSERT( input, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
390  ASSERT( output, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
391  ASSERT( gpsTime, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
392 
393  /* Make sure we're given the right coordinate system. */
394  if ( input->system != COORDINATESYSTEM_GEOGRAPHIC ) {
395  ABORT( stat, SKYCOORDINATESH_ESYS, SKYCOORDINATESH_MSGESYS );
396  }
397 
398  /* Compute the Greenwich mean sidereal time. */
399  gmst = fmod( XLALGreenwichMeanSiderealTime( gpsTime ), LAL_TWOPI );
400 
401  /* Subtract from longitude, and exit. */
403  output->longitude = fmod( input->longitude + gmst, LAL_TWOPI );
404  if ( output->longitude < 0.0 )
405  output->longitude += LAL_TWOPI;
406  output->latitude = input->latitude;
407  RETURN( stat );
408 }
409 
410 
411 /** \see See documentation in \ref TerrestrialCoordinates_c */
412 void
415  SkyPosition *input,
416  const SkyPosition *zenith )
417 {
418  REAL8 h, sinH, cosH; /* hour angle, and its sine and cosine */
419  REAL8 sinP, cosP; /* sin and cos of zenith latitude */
420  REAL8 sinD, cosD; /* sin and cos of position latitude (declination) */
421  REAL8 sina, sinA, cosA; /* sin and cos of altitude and -azimuth */
422 
423  INITSTATUS(stat);
424 
425  /* Make sure parameter structures exist. */
426  ASSERT( input, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
427  ASSERT( output, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
428  ASSERT( zenith, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
429 
430  /* Make sure we have consistent input coordinate systems. */
431  if ( input->system != zenith->system ) {
432  ABORT( stat, SKYCOORDINATESH_ESYS, SKYCOORDINATESH_MSGESYS );
433  }
434  if ( ( zenith->system != COORDINATESYSTEM_EQUATORIAL ) &&
435  ( zenith->system != COORDINATESYSTEM_GEOGRAPHIC ) ) {
436  ABORT( stat, SKYCOORDINATESH_ESYS, SKYCOORDINATESH_MSGESYS );
437  }
438 
439  /* Compute intermediates. */
440  h = zenith->longitude - input->longitude;
441  sinH = sin( h );
442  cosH = cos( h );
443  sinP = sin( zenith->latitude );
444  cosP = cos( zenith->latitude );
445  sinD = sin( input->latitude );
446  cosD = cos( input->latitude );
447 
448  /* Compute components. */
449  sina = sinD*sinP + cosD*cosP*cosH;
450  sinA = cosD*sinH;
451  cosA = sinD*cosP - cosD*sinP*cosH;
452 
453  /* Compute final results. */
455  output->latitude = asin( sina );
456  output->longitude = atan2( sinA, cosA );
457 
458  /* Optional phase correction. */
459  if ( output->longitude < 0.0 )
460  output->longitude += LAL_TWOPI;
461 
462  RETURN( stat );
463 }
464 
465 
466 /** \see See documentation in \ref TerrestrialCoordinates_c */
467 void
470  SkyPosition *input,
471  const SkyPosition *zenith )
472 {
473  REAL8 sinP, cosP; /* sin and cos of zenith latitude */
474  REAL8 sina, cosa; /* sin and cos of altitude */
475  REAL8 sinA, cosA; /* sin and cos of -azimuth */
476  REAL8 sinD, sinH, cosH; /* sin and cos of declination and hour angle */
477 
478  INITSTATUS(stat);
479 
480  /* Make sure parameter structures exist. */
481  ASSERT( input, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
482  ASSERT( output, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
483  ASSERT( zenith, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
484 
485  /* Make sure we're given the right coordinate system. */
486  if ( input->system != COORDINATESYSTEM_HORIZON ) {
487  ABORT( stat, SKYCOORDINATESH_ESYS, SKYCOORDINATESH_MSGESYS );
488  }
489  if ( ( zenith->system != COORDINATESYSTEM_EQUATORIAL ) &&
490  ( zenith->system != COORDINATESYSTEM_GEOGRAPHIC ) ) {
491  ABORT( stat, SKYCOORDINATESH_ESYS, SKYCOORDINATESH_MSGESYS );
492  }
493 
494  /* Compute intermediates. */
495  sinP = sin( zenith->latitude );
496  cosP = cos( zenith->latitude );
497  sina = sin( input->latitude );
498  cosa = cos( input->latitude );
499  sinA = sin( input->longitude );
500  cosA = cos( input->longitude );
501 
502  /* Compute components. */
503  sinD = sina*sinP + cosa*cosA*cosP;
504  sinH = cosa*sinA;
505  cosH = sina*cosP - cosa*cosA*sinP;
506 
507  /* Compute final results. */
508  output->system = zenith->system;
509  output->latitude = asin( sinD );
510  output->longitude = zenith->longitude - atan2( sinH, cosH );
511 
512  /* Optional phase correction. */
513  if ( output->longitude < 0.0 )
514  output->longitude += LAL_TWOPI;
515 
516  RETURN( stat );
517 }
518 
519 
520 /** \see See documentation in \ref TerrestrialCoordinates_c */
521 void
523 {
524  REAL8 c, s; /* position components in and orthogonal to the equator */
525  REAL8 cosP, sinP; /* cosine and sine of latitude */
526  REAL8 fFac; /* ( 1 - f )^2 */
527  REAL8 x, y; /* Cartesian coordinates */
528  REAL8 maxComp, r; /* max{x,y,z}, and sqrt(x^2+y^2+z^2) */
529 
530  INITSTATUS(stat);
531 
532  /* Make sure parameter structure exists. */
533  ASSERT( location, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
534 
535  /* Make sure the geodetic coordinates are in the right system. */
536  if ( location->geodetic.system != COORDINATESYSTEM_GEOGRAPHIC ) {
537  ABORT( stat, SKYCOORDINATESH_ESYS, SKYCOORDINATESH_MSGESYS );
538  }
539 
540  /* Compute intermediates. */
541  fFac = 1.0 - LAL_EARTHFLAT;
542  fFac *= fFac;
543  cosP = cos( location->geodetic.latitude );
544  sinP = sin( location->geodetic.latitude );
545  c = sqrt( 1.0 / ( cosP*cosP + fFac*sinP*sinP ) );
546  s = fFac*c;
547  c = ( LAL_REARTH_SI*c + location->elevation )*cosP;
548  s = ( LAL_REARTH_SI*s + location->elevation )*sinP;
549 
550  /* Compute Cartesian coordinates. */
551  location->x = x = c*cos( location->geodetic.longitude );
552  location->y = y = c*sin( location->geodetic.longitude );
553  location->z = s;
554 
555  /* Compute the radius. */
556  maxComp = x;
557  if ( y > maxComp )
558  maxComp = y;
559  if ( s > maxComp )
560  maxComp = s;
561  x /= maxComp;
562  y /= maxComp;
563  s /= maxComp;
564  r = sqrt( x*x + y*y + s*s );
565 
566  /* Compute the spherical coordinates, and exit. */
567  location->radius = maxComp*r;
568  location->geocentric.longitude = location->geodetic.longitude;
569  location->geocentric.latitude = asin( s / r );
571  RETURN( stat );
572 }
573 
574 
575 /** \see See documentation in \ref TerrestrialCoordinates_c */
576 void
578 {
579  REAL8 x, y, z; /* normalized geocentric coordinates */
580  REAL8 pi; /* axial distance */
581 
582  /* Declare some local constants. */
583  const REAL8 rInv = 1.0 / LAL_REARTH_SI;
584  const REAL8 f1 = 1.0 - LAL_EARTHFLAT;
585  const REAL8 f2 = LAL_EARTHFLAT*( 2.0 - LAL_EARTHFLAT );
586 
587  INITSTATUS(stat);
588 
589  /* Make sure parameter structure exists. */
590  ASSERT( location, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
591 
592  /* See if we've been given a special set of coordinates. */
593  x = rInv*location->x;
594  y = rInv*location->y;
595  z = rInv*location->z;
596  pi = sqrt( x*x + y*y );
598  location->geodetic.longitude = atan2( y, x );
599  location->radius = LAL_REARTH_SI*sqrt( x*x + y*y + z*z );
600  if ( pi == 0.0 ) {
601  if ( z >= 0.0 ) {
602  location->geodetic.latitude = LAL_PI_2;
603  location->elevation = z - f1;
604  } else {
605  location->geodetic.latitude = -LAL_PI_2;
606  location->elevation = -z - f1;
607  }
608  location->elevation *= LAL_REARTH_SI;
609  }
610 
611  /* Do the general transformation even if z=0. */
612  else {
613  REAL8 za, e, f, p, p3, q, q2, d, d2, b, v, w, g, h, gh, t, phi, tanP;
614  /* intermediate variables */
615 
616  /* See if we're inside the singular ellipsoid.
617  if ( pi <= 1.01*f2 ) {
618  REAL8 z1 = z*f1/f2;
619  REAL8 p1 = pi/f2;
620  if ( z1*z1 + p1*p1 < 1.02 ) {
621  ABORT( stat, SKYCOORDINATESH_ESING, SKYCOORDINATESH_MSGESING );
622  }
623  } */
624 
625  /* Compute intermediates variables. */
626  za = f1*fabs( z );
627  e = za - f2;
628  f = za + f2;
629  p = ( 4.0/3.0 )*( pi*pi + za*za - f2*f2 );
630  p3 = p*p*p;
631  q = 8.0*pi*f2*za;
632  q2 = q*q;
633  h = p3/q2;
634  d = p3 + q2;
635 
636  /* Compute v, using series expansion if necessary. */
637  if ( d >= 0.0 ) {
638  d2 = sqrt( d );
639  b = q/d2;
640  if ( fabs( h ) < LAL_HSERIES ) {
641  if ( h < 0.0 )
642  v = pow( d2 + q, 1.0/3.0 )
643  + pow( -0.5*q*h*( 1.0 - 0.25*h*( 1.0 - 0.5*h ) ), 1.0/3.0 );
644  else
645  v = pow( d2 + q, 1.0/3.0 )
646  - pow( 0.5*q*h*( 1.0 - 0.25*h*( 1.0 - 0.5*h ) ), 1.0/3.0 );
647  } else if ( b < LAL_BSERIES ) {
648  v = pow( d2, 1.0/3.0 )*( b*( 2.0/3.0 + b*b*10.0/81.0 ) );
649  } else if ( b < 1.0 ) {
650  v = pow( d2, 1.0/3.0 )*( pow( 1.0 + b, 1.0/3.0 ) -
651  pow( 1.0 - b, 1.0/3.0 ) );
652  } else {
653  v = pow( d2, 1.0/3.0 )*( pow( b + 1.0, 1.0/3.0 ) +
654  pow( b - 1.0, 1.0/3.0 ) );
655  }
656  } else {
657  v = 2.0*sqrt( -p )*cos( acos( q/( -p*sqrt( -p ) ) )/3.0 );
658  }
659 
660  /* Compute t, using series expansion if necessary. */
661  h = v*pi/( e*e );
662  w = fabs( e )*sqrt( 1.0 + h );
663  if ( e > 0.0 || h > LAL_HSERIES )
664  g = 0.5*( e + w );
665  else
666  g = -0.25*e*h*( 1.0 - 0.25*h*( 1.0 - 0.5*h ) );
667  gh = fabs( pi*( f*pi - v*g )/w );
668  h = gh/( g*g );
669  if ( fabs( h ) > LAL_HSERIES )
670  t = g*( sqrt( 1.0 + h ) - 1.0 );
671  else
672  t = 0.5*g*h*( 1.0 - 0.25*h*( 1.0 - 0.5*h ) );
673 
674  /* Compute latitude, longitude, and elevation. */
675  tanP = ( pi - t )*( pi + t )/( 2.0*f1*pi*t );
676  phi = atan( tanP );
677  location->geodetic.latitude = phi;
678  if ( z < 0.0 )
679  location->geodetic.latitude *= -1.0;
680  location->elevation = ( pi - t/pi )*cos( phi );
681  location->elevation += ( fabs( z ) - f1 )*sin( phi );
682  location->elevation *= LAL_REARTH_SI;
683  }
684 
685  /* Transformation complete. */
686  RETURN( stat );
687 }
688 /** @} */
#define ABORT(statusptr, code, mesg)
#define ASSERT(assertion, statusptr, code, mesg)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
#define LAL_EARTHFLAT
#define LAL_HSERIES
#define LAL_BSERIES
static double f(double theta, double y, double xi)
Definition: XLALMarcumQ.c:258
#define LAL_PI_2
pi/2
Definition: LALConstants.h:181
#define LAL_REARTH_SI
Earth equatorial radius, m.
Definition: LALConstants.h:416
#define LAL_TWOPI
2*pi is circumference of a circle divided by its radius
Definition: LALConstants.h:180
double REAL8
Double precision real floating-point number (8 bytes).
static const INT4 r
Definition: Random.c:82
static const INT4 q
Definition: Random.c:81
#define SKYCOORDINATESH_ENUL
Unexpected null pointer in arguments.
#define SKYCOORDINATESH_ESYS
Wrong coordinate system in input.
@ COORDINATESYSTEM_GEOGRAPHIC
The Earth-fixed geographic coordinate system.
@ COORDINATESYSTEM_HORIZON
A horizon coordinate system.
@ COORDINATESYSTEM_EQUATORIAL
The sky-fixed equatorial coordinate system.
void LALGeocentricToGeodetic(LALStatus *stat, EarthPosition *location)
void LALEquatorialToGeographic(LALStatus *stat, SkyPosition *output, SkyPosition *input, LIGOTimeGPS *gpsTime)
void LALHorizonToSystem(LALStatus *stat, SkyPosition *output, SkyPosition *input, const SkyPosition *zenith)
void LALGeographicToEquatorial(LALStatus *stat, SkyPosition *output, SkyPosition *input, LIGOTimeGPS *gpsTime)
void LALGeodeticToGeocentric(LALStatus *stat, EarthPosition *location)
void LALSystemToHorizon(LALStatus *stat, SkyPosition *output, SkyPosition *input, const SkyPosition *zenith)
REAL8 XLALGreenwichMeanSiderealTime(const LIGOTimeGPS *gpstime)
Convenience wrapper, calling XLALGreenwichSiderealTime() with the equation of equinoxes set to 0.
This structure stores the location of a point on (or near) the surface of the Earth in both geodetic ...
SkyPosition geocentric
The geographic coordinates of the direction from the centre of the Earth through the point; that is,...
REAL8 radius
The distance of the point from the geocentre, in metres.
REAL8 elevation
The vertical distance of the point above the reference ellipsoid, in metres.
REAL8 z
The Earth-fixed geocentric Cartesian coordinates of the point, in metres.
SkyPosition geodetic
The geographic coordinates of the upward vertical direction from the point; that is,...
LAL status structure, see The LALStatus structure for more details.
Definition: LALDatatypes.h:947
Epoch relative to GPS epoch, see LIGOTimeGPS type for more details.
Definition: LALDatatypes.h:458
This structure stores the two spherical coordinates of a sky position; ie a generic latitude and long...
REAL8 longitude
The longitudinal coordinate (in radians), as defined above.
REAL8 latitude
The latitudinal coordinate (in radians), as defined above.
CoordinateSystem system
The coordinate system in which latitude/longitude are expressed.
void output(int gps_sec, int output_type)
Definition: tconvert.c:440