Loading [MathJax]/extensions/TeX/AMSsymbols.js
LAL 7.7.0.1-3a66518
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
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 */
148void
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 */
196void
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 */
244void
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 */
289void
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_ALPHAGAL
#define LAL_DELTAGAL
#define LAL_LGAL
#define ABORT(statusptr, code, mesg)
#define ASSERT(assertion, statusptr, code, mesg)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
void LALEquatorialToGalactic(LALStatus *stat, SkyPosition *output, SkyPosition *input)
void LALEquatorialToEcliptic(LALStatus *stat, SkyPosition *output, SkyPosition *input)
void LALEclipticToEquatorial(LALStatus *stat, SkyPosition *output, SkyPosition *input)
void LALGalacticToEquatorial(LALStatus *stat, SkyPosition *output, SkyPosition *input)
#define LAL_IEARTH
Earth inclination (2000), radians.
Definition: LALConstants.h:445
#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 a
Definition: Random.c:79
#define SKYCOORDINATESH_ENUL
Unexpected null pointer in arguments.
#define SKYCOORDINATESH_ESYS
Wrong coordinate system in input.
@ COORDINATESYSTEM_GALACTIC
The galactic coordinate system.
@ COORDINATESYSTEM_ECLIPTIC
The ecliptic coordinate system.
@ COORDINATESYSTEM_EQUATORIAL
The sky-fixed equatorial coordinate system.
LAL status structure, see The LALStatus structure for more details.
Definition: LALDatatypes.h:947
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