LAL  7.1.7.1-ab2cc12
CelestialCoordinates.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 #include <math.h>
21 #include <lal/LALStdlib.h>
22 #include <lal/LALConstants.h>
23 #include <lal/SkyCoordinates.h>
24 
25 #define LAL_ALPHAGAL (3.366032942)
26 #define LAL_DELTAGAL (0.473477302)
27 #define LAL_LGAL (0.576)
28 
29 /**
30  * \author T.D. Creighton
31  * \addtogroup CelestialCoordinates_c
32  * \brief Converts among Galactic, ecliptic, and equatorial coordinates.
33  *
34  * These functions perform the specified coordinate transformation on the
35  * contents of \a input and store the result in \a *output. The
36  * two pointers may point to the same object, in which case the
37  * conversion is done in place. The functions will also check
38  * <tt>input->system</tt> and set <tt>output->system</tt> as appropriate.
39  *
40  * These routines are collected together because they involve fixed,
41  * absolute coordinate systems, so the transformations require no
42  * additional parameters such as the time or site of observation. We
43  * also note that there are no direct conversions between Galactic and
44  * ecliptic coordinates. At the risk of additional computational
45  * overhead, it is simple to use the equatorial coordinate system as an
46  * intermediate step.
47  *
48  * ### Description ###
49  *
50  * These functions perform the specified coordinate transformation on the
51  * contents of \a input and store the result in \a output. The
52  * two pointers may point to the same object, in which case the
53  * conversion is done in place. The functions will also check
54  * <tt>input->system</tt> and set <tt>output->system</tt> as appropriate.
55  *
56  * These routines are collected together because they involve fixed,
57  * absolute coordinate systems, so the transformations require no
58  * additional parameters such as the time or site of observation. We
59  * also note that there are no direct conversions between Galactic and
60  * ecliptic coordinates. At the risk of additional computational
61  * overhead, it is simple to use the equatorial coordinate system as an
62  * intermediate step.
63  *
64  * ### Algorithm ###
65  *
66  * These routines follow the spherical angle relations on p. 13
67  * of \cite Lang_K1999 . Note that the actual formulae for Galactic
68  * longitude and right ascension in this reference are wrong; we give
69  * corrected formulae below derived from the sine and cosine equations.
70  * (The Galactic to equatorial transformations can also be found in
71  * Sec. 12.3 of \cite GRASP2000 . All positions are assumed to
72  * be in the J2000 epoch.
73  *
74  * \par Galactic coordinates:
75  * The following formulae relate
76  * Galactic latitude \f$b\f$ and longitude \f$l\f$ to declination \f$\delta\f$ and
77  * right ascension \f$\alpha\f$:
78  * \f{eqnarray}{
79  * b & = & \arcsin[\cos\delta\cos\delta_\mathrm{NGP}
80  * \cos(\alpha-\alpha_\mathrm{NGP}) +
81  * \sin\delta\sin\delta_\mathrm{NGP}] \;,\\
82  * l & = & \arctan\!2[\sin\delta\cos\delta_\mathrm{NGP} -
83  * \cos\delta\cos(\alpha-\alpha_\mathrm{NGP})
84  * \sin\delta_\mathrm{NGP},
85  * \cos\delta\sin(\alpha-\alpha_\mathrm{NGP})]\\
86  * & & \quad + \; l_\mathrm{ascend} \;,
87  * \f}
88  * where \f$\arctan\!2(y,x)\f$ can be thought of as the argument of the
89  * complex number \f$x+iy\f$; unlike \f$\arctan(y/x)\f$, it ranges over the full
90  * range \f$[0,2\pi)\f$ instead of just half of it. The inverse
91  * transformations are:
92  * \f{eqnarray}{
93  * \delta & = & \arcsin[\cos b\cos\delta_\mathrm{NGP}\sin(l-l_\mathrm{ascend}) +
94  * \sin b\sin\delta_\mathrm{NGP}] \;,\\
95  * \alpha & = & \arctan\!2[\cos b\cos(l-l_\mathrm{ascend}),
96  * \sin b\cos\delta_\mathrm{NGP} -
97  * \cos b\sin(l-l_\mathrm{ascend})\sin\delta_\mathrm{NGP}] \\
98  * & & \quad + \; \alpha_\mathrm{NGP} \;.
99  * \f}
100  * In these equations we have defined the orientation of the Galaxy with
101  * the following parameters (which should eventually be placed in <tt>LALConstants.h</tt>:
102  * \f[
103  * \begin{array}{r@{\quad=\quad}l@{\quad=\quad}l}
104  * \alpha_\mathrm{NGP} & 192.8594813^\circ &
105  * \mbox{the right ascension (epoch J2000) of the north Galactic pole} \\
106  * \delta_\mathrm{NGP} & 27.1282511^\circ &
107  * \mbox{the declination (epoch J2000) of the north Galactic pole} \\
108  * l_\mathrm{ascend} & 33^\circ &
109  * \mbox{the longitude of the ascending node of the Galactic plane}
110  * \end{array}
111  * \f]
112  * The ascending node of the Galactic plane is defined as the direction
113  * along the intersection of the Galactic and equatorial planes where
114  * rotation in the positive sense about the Galactic \f$z\f$ axis carries a
115  * point from the southern to northern equatorial hemisphere. That is,
116  * if \f$\mathbf{u}\f$ points in the direction \f$\delta=90^\circ\f$
117  * (celestial north), and \f$\mathbf{v}\f$ points in the direction
118  * \f$b=90^\circ\f$ (Galactic north), then
119  * \f$\mathbf{u} \times \mathbf{v}\f$ points along the ascending node.
120  *
121  * \par Ecliptic coordinates:
122  * The following formulae relate
123  * Ecliptic latitude \f$\beta\f$ and longitude \f$\lambda\f$ to declination
124  * \f$\delta\f$ and right ascension \f$\alpha\f$:
125  * \f{eqnarray}{
126  * \beta & = & \arcsin(\sin\delta\cos\epsilon -
127  * \cos\delta\sin\alpha\sin\epsilon) \;, \\
128  * \lambda & = & \arctan\!2(\cos\delta\sin\alpha\cos\epsilon +
129  * \sin\delta\sin\epsilon, \cos\delta\cos\alpha) \;.
130  * \f}
131  * The inverse transformations are:
132  * \f{eqnarray}{
133  * \delta & = & \arcsin(\cos\beta\sin\lambda\sin\epsilon +
134  * \sin\beta\cos\epsilon) \;, \\
135  * \alpha & = & \arctan\!2(\cos\beta\sin\lambda\cos\epsilon -
136  * \sin\beta\sin\epsilon, \cos\beta\cos\lambda) \;.
137  * \f}
138  * Here \f$\epsilon\f$ is the obliquity (inclination) of the ecliptic plane,
139  * which varies over time; at epoch J200 it has a mean value of:
140  * \f[
141  * \epsilon = 23.4392911^\circ \; .
142  * \f]
143  *
144  */
145 /** @{ */
146 
147 /** \see See documentation in \ref CelestialCoordinates_c */
148 void
150  SkyPosition *output,
151  SkyPosition *input )
152 {
153  REAL8 sinDGal = sin( LAL_DELTAGAL ); /* sin(delta_NGP) */
154  REAL8 cosDGal = cos( LAL_DELTAGAL ); /* cos(delta_NGP) */
155  REAL8 l = -LAL_LGAL; /* will be l-l(ascend) */
156  REAL8 sinB, cosB, sinL, cosL; /* sin and cos of b and l */
157  REAL8 sinD, sinA, cosA; /* sin and cos of delta and alpha */
158 
159  INITSTATUS(stat);
160 
161  /* Make sure parameter structures exist. */
162  ASSERT( input, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
163  ASSERT( output, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
164 
165  /* Make sure we're given the right coordinate system. */
166  if ( input->system != COORDINATESYSTEM_GALACTIC ) {
167  ABORT( stat, SKYCOORDINATESH_ESYS, SKYCOORDINATESH_MSGESYS );
168  }
169 
170  /* Compute intermediates. */
171  l += input->longitude;
172  sinB = sin( input->latitude );
173  cosB = cos( input->latitude );
174  sinL = sin( l );
175  cosL = cos( l );
176 
177  /* Compute components. */
178  sinD = cosB*cosDGal*sinL + sinB*sinDGal;
179  sinA = cosB*cosL;
180  cosA = sinB*cosDGal - cosB*sinL*sinDGal;
181 
182  /* Compute final results. */
184  output->latitude = asin( sinD );
185  l = atan2( sinA, cosA ) + LAL_ALPHAGAL;
186 
187  /* Optional phase correction. */
188  if ( l < 0.0 )
189  l += LAL_TWOPI;
190  output->longitude = l;
191 
192  RETURN( stat );
193 }
194 
195 /** \see See documentation in \ref CelestialCoordinates_c */
196 void
198  SkyPosition *output,
199  SkyPosition *input )
200 {
201  REAL8 sinDGal = sin( LAL_DELTAGAL ); /* sin(delta_NGP) */
202  REAL8 cosDGal = cos( LAL_DELTAGAL ); /* cos(delta_NGP) */
203  REAL8 a = -LAL_ALPHAGAL; /* will be alpha-alpha_NGP */
204  REAL8 sinD, cosD, sinA, cosA; /* sin and cos of delta and alpha */
205  REAL8 sinB, sinL, cosL; /* sin and cos of b and l */
206 
207  INITSTATUS(stat);
208 
209  /* Make sure parameter structures exist. */
210  ASSERT( input, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
211  ASSERT( output, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
212 
213  /* Make sure we're given the right coordinate system. */
214  if ( input->system != COORDINATESYSTEM_EQUATORIAL ) {
215  ABORT( stat, SKYCOORDINATESH_ESYS, SKYCOORDINATESH_MSGESYS );
216  }
217 
218  /* Compute intermediates. */
219  a += input->longitude;
220  sinD = sin( input->latitude );
221  cosD = cos( input->latitude );
222  sinA = sin( a );
223  cosA = cos( a );
224 
225  /* Compute components. */
226  sinB = cosD*cosDGal*cosA + sinD*sinDGal;
227  sinL = sinD*cosDGal - cosD*cosA*sinDGal;
228  cosL = cosD*sinA;
229 
230  /* Compute final results. */
232  output->latitude = asin( sinB );
233  a = atan2( sinL, cosL ) + LAL_LGAL;
234 
235  /* Optional phase correction. */
236  if ( a < 0.0 )
237  a += LAL_TWOPI;
238  output->longitude = a;
239 
240  RETURN( stat );
241 }
242 
243 /** \see See documentation in \ref CelestialCoordinates_c */
244 void
246  SkyPosition *output,
247  SkyPosition *input )
248 {
249  REAL8 sinE = sin( LAL_IEARTH ); /* sin(epsilon) */
250  REAL8 cosE = cos( LAL_IEARTH ); /* cos(epsilon) */
251  REAL8 sinB, cosB, sinL, cosL; /* sin and cos of b and l */
252  REAL8 sinD, sinA, cosA; /* sin and cos of delta and alpha */
253 
254  INITSTATUS(stat);
255 
256  /* Make sure parameter structures exist. */
257  ASSERT( input, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
258  ASSERT( output, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
259 
260  /* Make sure we're given the right coordinate system. */
261  if ( input->system != COORDINATESYSTEM_ECLIPTIC ) {
262  ABORT( stat, SKYCOORDINATESH_ESYS, SKYCOORDINATESH_MSGESYS );
263  }
264 
265  /* Compute intermediates. */
266  sinB = sin( input->latitude );
267  cosB = cos( input->latitude );
268  sinL = sin( input->longitude );
269  cosL = cos( input->longitude );
270 
271  /* Compute components. */
272  sinD = cosB*sinL*sinE + sinB*cosE;
273  sinA = cosB*sinL*cosE - sinB*sinE;
274  cosA = cosB*cosL;
275 
276  /* Compute final results. */
278  output->latitude = asin( sinD );
279  output->longitude = atan2( sinA, cosA );
280 
281  /* Optional phase correction. */
282  if ( output->longitude < 0.0 )
283  output->longitude += LAL_TWOPI;
284 
285  RETURN( stat );
286 }
287 
288 /** \see See documentation in \ref CelestialCoordinates_c */
289 void
291  SkyPosition *output,
292  SkyPosition *input )
293 {
294  REAL8 sinE = sin( LAL_IEARTH ); /* sin(epsilon) */
295  REAL8 cosE = cos( LAL_IEARTH ); /* cos(epsilon) */
296  REAL8 sinD, cosD, sinA, cosA; /* sin and cos of delta and alpha */
297  REAL8 sinB, sinL, cosL; /* sin and cos of b and l */
298 
299  INITSTATUS(stat);
300 
301  /* Make sure parameter structures exist. */
302  ASSERT( input, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
303  ASSERT( output, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
304 
305  /* Make sure we're given the right coordinate system. */
306  if ( input->system != COORDINATESYSTEM_EQUATORIAL ) {
307  ABORT( stat, SKYCOORDINATESH_ESYS, SKYCOORDINATESH_MSGESYS );
308  }
309 
310  /* Compute intermediates. */
311  sinD = sin( input->latitude );
312  cosD = cos( input->latitude );
313  sinA = sin( input->longitude );
314  cosA = cos( input->longitude );
315 
316  /* Compute components. */
317  sinB = sinD*cosE - cosD*sinA*sinE;
318  sinL = cosD*sinA*cosE + sinD*sinE;
319  cosL = cosD*cosA;
320 
321  /* Compute final results. */
323  output->latitude = asin( sinB );
324  output->longitude = atan2( sinL, cosL );
325 
326  /* Optional phase correction. */
327  if ( output->longitude < 0.0 )
328  output->longitude += LAL_TWOPI;
329 
330  RETURN( stat );
331 }
332 /** @} */
#define LAL_IEARTH
Earth inclination (2000), radians.
Definition: LALConstants.h:433
void LALEclipticToEquatorial(LALStatus *stat, SkyPosition *output, SkyPosition *input)
The galactic coordinate system.
#define LAL_ALPHAGAL
#define LAL_LGAL
The sky-fixed equatorial coordinate system.
REAL8 latitude
The latitudinal coordinate (in radians), as defined above.
void LALEquatorialToEcliptic(LALStatus *stat, SkyPosition *output, SkyPosition *input)
void LALGalacticToEquatorial(LALStatus *stat, SkyPosition *output, SkyPosition *input)
#define SKYCOORDINATESH_ESYS
Wrong coordinate system in input.
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.
static const INT4 a
Definition: Random.c:79
#define LAL_DELTAGAL
#define RETURN(statusptr)
The ecliptic coordinate system.
#define INITSTATUS(statusptr)
void LALEquatorialToGalactic(LALStatus *stat, SkyPosition *output, SkyPosition *input)
double REAL8
Double precision real floating-point number (8 bytes).
#define ABORT(statusptr, code, mesg)
LAL status structure, see The LALStatus structure for more details.
Definition: LALDatatypes.h:947
REAL8 longitude
The longitudinal coordinate (in radians), as defined above.
#define SKYCOORDINATESH_ENUL
Unexpected null pointer in arguments.
#define ASSERT(assertion, statusptr, code, mesg)
#define LAL_TWOPI
2*pi is circumference of a circle divided by its radius
Definition: LALConstants.h:180