LAL  7.1.7.1-9b0afb2
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
346  SkyPosition *output,
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
380  SkyPosition *output,
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
414  SkyPosition *output,
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
469  SkyPosition *output,
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 /** @} */
void LALGeodeticToGeocentric(LALStatus *stat, EarthPosition *location)
Epoch relative to GPS epoch, see LIGOTimeGPS type for more details.
Definition: LALDatatypes.h:458
void LALEquatorialToGeographic(LALStatus *stat, SkyPosition *output, SkyPosition *input, LIGOTimeGPS *gpsTime)
REAL8 XLALGreenwichMeanSiderealTime(const LIGOTimeGPS *gpstime)
Convenience wrapper, calling XLALGreenwichSiderealTime() with the equation of equinoxes set to 0...
REAL8 radius
The distance of the point from the geocentre, in metres.
The sky-fixed equatorial coordinate system.
REAL8 elevation
The vertical distance of the point above the reference ellipsoid, in metres.
REAL8 latitude
The latitudinal coordinate (in radians), as defined above.
#define SKYCOORDINATESH_ESYS
Wrong coordinate system in input.
SkyPosition geocentric
The geographic coordinates of the direction from the centre of the Earth through the point; that is...
This structure stores the two spherical coordinates of a sky position; ie a generic latitude and long...
CoordinateSystem system
The coordinate system in which latitude/longitude are expressed.
void LALGeographicToEquatorial(LALStatus *stat, SkyPosition *output, SkyPosition *input, LIGOTimeGPS *gpsTime)
#define LAL_BSERIES
void LALGeocentricToGeodetic(LALStatus *stat, EarthPosition *location)
static double f(double theta, double y, double xi)
Definition: XLALMarcumQ.c:258
#define RETURN(statusptr)
#define INITSTATUS(statusptr)
The Earth-fixed geographic coordinate system.
double REAL8
Double precision real floating-point number (8 bytes).
#define ABORT(statusptr, code, mesg)
A horizon coordinate system.
LAL status structure, see The LALStatus structure for more details.
Definition: LALDatatypes.h:947
#define LAL_REARTH_SI
Earth equatorial radius, m.
Definition: LALConstants.h:404
void LALHorizonToSystem(LALStatus *stat, SkyPosition *output, SkyPosition *input, const SkyPosition *zenith)
#define LAL_EARTHFLAT
static const INT4 q
Definition: Random.c:81
REAL8 longitude
The longitudinal coordinate (in radians), as defined above.
This structure stores the location of a point on (or near) the surface of the Earth in both geodetic ...
#define SKYCOORDINATESH_ENUL
Unexpected null pointer in arguments.
#define LAL_HSERIES
static const INT4 r
Definition: Random.c:82
SkyPosition geodetic
The geographic coordinates of the upward vertical direction from the point; that is, the point&#39;s geodetic latitude and longitude.
REAL8 z
The Earth-fixed geocentric Cartesian coordinates of the point, in metres.
void LALSystemToHorizon(LALStatus *stat, SkyPosition *output, SkyPosition *input, const SkyPosition *zenith)
#define ASSERT(assertion, statusptr, code, mesg)
#define LAL_TWOPI
2*pi is circumference of a circle divided by its radius
Definition: LALConstants.h:180
#define LAL_PI_2
pi/2
Definition: LALConstants.h:181