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
SkyCoordinates.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 Creighton, T. D.
31 * \addtogroup SkyCoordinates_c
32 * \brief Automatically converts among sky coordinate systems.
33 *
34 * The function <tt>LALConvertSkyCoordinates()</tt> transforms the contents
35 * of <tt>*input</tt> to the system
36 * specified in <tt>*params</tt>, storing the result in <tt>*output</tt>
37 * (which may point to the same object as <tt>*input</tt> for an in-place
38 * transformation). The routine makes calls to the functions in
39 * \ref CelestialCoordinates.c and \ref TerrestrialCoordinates.c as
40 * required; the <tt>*params</tt> object must store any data fields
41 * required by these functions, or an error will occur.
42 *
43 * The function <tt>LALNormalizeSkyPosition()</tt> "normalizes" any given
44 * (spherical) sky-position (in radians), which means it projects the
45 * angles into \f$[0, 2\pi) \times [-\pi/2, \pi/2]\f$ if they lie outside.
46 *
47 * ### Algorithm ###
48 *
49 * <tt>LALConvertSkyCoordinates()</tt> is structured as a simple loop over
50 * transformations, each
51 * of which moves the output sky position one step closer to the desired
52 * final coordinates system. The usual "flow" of the algorithm is:
53 *
54 * \anchor SkyCoordinates_conversions
55 * \image html SkyCoordinates_conversions.png ""
56 *
57 * although one can also convert directly between equatorial and horizon
58 * coordinate systems if <tt>params->zenith</tt> is given in equatorial
59 * coordinates (i.e.\ if its longitudinal coordinate is the local mean
60 * sidereal time rather than the geographic longitude of the observer).
61 * This leads to the only error checking done within this function: when
62 * transforming to horizon coordinates, it checks that
63 * <tt>params->zenith</tt> is either in sky-fixed equatorial or Earth-fixed
64 * geographic coordinates. Other than this, error checking is left to
65 * the secondary function call; if a parameter is absent or poorly
66 * formatted, the called function will return an error.
67 *
68 */
69/** @{ */
70
71/** \see See documentation in \ref SkyCoordinates_c */
72void
75 SkyPosition *input,
76 ConvertSkyParams *params )
77{
78 SkyPosition temp; /* temporary sky position (duh) */
79
80 INITSTATUS(stat);
81 ATTATCHSTATUSPTR( stat );
82
83 /* Make sure parameter structures exist. */
84 ASSERT( input, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
85 ASSERT( output, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
86 ASSERT( params, stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
87
88 /* Start looping! */
89 temp = *input;
90 while ( temp.system != params->system ) {
91
92 if ( temp.system == COORDINATESYSTEM_HORIZON )
93 /* Only one possible direction. */
94 TRY( LALHorizonToSystem( stat->statusPtr, &temp, &temp,
95 params->zenith ), stat );
96
97 else if ( temp.system == COORDINATESYSTEM_GEOGRAPHIC ) {
98 /* Two possible directions. But, if headed towards the horizon
99 coordinate system, make sure we can get there! */
100 if ( params->system == COORDINATESYSTEM_HORIZON ) {
101 if ( params->zenith ) {
102 if ( params->zenith->system == COORDINATESYSTEM_GEOGRAPHIC )
103 TRY( LALSystemToHorizon( stat->statusPtr, &temp, &temp,
104 params->zenith ), stat );
105 else if ( params->zenith->system == COORDINATESYSTEM_EQUATORIAL )
106 TRY( LALGeographicToEquatorial( stat->statusPtr, &temp, &temp,
107 params->gpsTime ), stat );
108 else
109 ABORT( stat, SKYCOORDINATESH_ESYS, SKYCOORDINATESH_MSGESYS );
110 } else
111 ABORT( stat, SKYCOORDINATESH_ENUL, SKYCOORDINATESH_MSGENUL );
112 } else
113 TRY( LALGeographicToEquatorial( stat->statusPtr, &temp, &temp,
114 params->gpsTime ), stat );
115 }
116
117 else if ( temp.system == COORDINATESYSTEM_EQUATORIAL ) {
118 /* Up to four possible directions, depending on
119 params->zenith->system. */
120 if ( ( params->system == COORDINATESYSTEM_HORIZON ) &&
121 ( params->zenith ) &&
123 TRY( LALSystemToHorizon( stat->statusPtr, &temp, &temp,
124 params->zenith ), stat );
125 else if ( params->system == COORDINATESYSTEM_ECLIPTIC )
126 TRY( LALEquatorialToEcliptic( stat->statusPtr, &temp, &temp ),
127 stat );
128 else if ( params->system == COORDINATESYSTEM_GALACTIC )
129 TRY( LALEquatorialToGalactic( stat->statusPtr, &temp, &temp ),
130 stat );
131 else
132 TRY( LALEquatorialToGeographic( stat->statusPtr, &temp, &temp,
133 params->gpsTime ), stat );
134 }
135
136 else if ( temp.system == COORDINATESYSTEM_ECLIPTIC ) {
137 /* Only one possible direction. */
138 TRY( LALEclipticToEquatorial( stat->statusPtr, &temp, &temp ),
139 stat );
140 }
141
142 else if ( temp.system == COORDINATESYSTEM_GALACTIC ) {
143 /* Only one possible direction. */
144 TRY( LALGalacticToEquatorial( stat->statusPtr, &temp, &temp ),
145 stat );
146 }
147
148 else
149 ABORT( stat, SKYCOORDINATESH_ESYS, SKYCOORDINATESH_MSGESYS );
150 }
151
152 /* We've gotten to the correct coordinate system. Set the output
153 and return. */
154 *output = temp;
155 DETATCHSTATUSPTR( stat );
156 RETURN( stat );
157
158} /* LALConvertSkyCoordinates */
159
160
161
162/**
163 * \deprecated Use XLALNormalizeSkyPosition() instead.
164 */
165void
166LALNormalizeSkyPosition (LALStatus *stat, /**< pointer to LALStatus structure */
167 SkyPosition *posOut, /**< [out] normalized sky-position */
168 const SkyPosition *posIn) /**< [in] general sky-position */
169{
170 SkyPosition tmp; /* allow posOut == posIn */
171
172 INITSTATUS(stat);
173
174 ASSERT (posIn, stat, SKYCOORDINATESH_ENUL , SKYCOORDINATESH_MSGENUL );
175 ASSERT (posOut, stat, SKYCOORDINATESH_ENUL , SKYCOORDINATESH_MSGENUL );
176
177 tmp = *posIn;
178
180
181 *posOut = tmp;
182
183 RETURN (stat);
184
185} /* LALNormalizeSkyPosition() */
186
187
188/**
189 * If sky-position is not in the canonical range
190 * \f$(\alpha,\delta)\in [0,2\pi]\times[-\pi/2, \pi/2]\f$, normalize it
191 * by mapping it into this coordinate-interval.
192 * Based on Alicia's function with some additional "unwinding" added.
193 */
194void
195XLALNormalizeSkyPosition ( double *restrict longitude, /**< [in,out] sky-position longitude to normalize*/
196 double *restrict latitude /**< [in,out] sky-position latitude to normalize*/
197 )
198{
199
200 /* FIRST STEP: completely "unwind" positions, i.e. make sure that
201 * [0 <= alpha < 2pi] and [-pi < delta <= pi] */
202 /* normalize longitude */
203 *longitude -= floor(*longitude / LAL_TWOPI) * LAL_TWOPI;
204
205 /* pre-normalize (unwind) latitude */
206 *latitude += LAL_PI;
207 *latitude -= floor(*latitude / LAL_TWOPI) * LAL_TWOPI;
208 *latitude -= LAL_PI;
209
210 /* SECOND STEP: get latitude into canonical interval [-pi/2 <= delta <= pi/2 ] */
211 /* this requires also a change in longitude by adding/subtracting PI */
212 if (*latitude > LAL_PI_2)
213 {
214 *latitude = LAL_PI - *latitude;
215 if (*longitude < LAL_PI)
216 {
217 *longitude += LAL_PI;
218 }
219 else
220 {
221 *longitude -= LAL_PI;
222 }
223 }
224
225 if (*latitude < -LAL_PI_2)
226 {
227 *latitude = -LAL_PI - *latitude;
228 if (*longitude < LAL_PI)
229 {
230 *longitude += LAL_PI;
231 }
232 else
233 {
234 *longitude -= LAL_PI;
235 }
236 }
237
238 return;
239
240} /* XLALNormalizeSkyPosition() */
241
242/** @} */
#define ABORT(statusptr, code, mesg)
#define TRY(func, statusptr)
#define ATTATCHSTATUSPTR(statusptr)
#define ASSERT(assertion, statusptr, code, mesg)
#define DETATCHSTATUSPTR(statusptr)
#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_PI_2
pi/2
Definition: LALConstants.h:181
#define LAL_PI
Archimedes's constant, pi.
Definition: LALConstants.h:179
#define LAL_TWOPI
2*pi is circumference of a circle divided by its radius
Definition: LALConstants.h:180
void LALConvertSkyCoordinates(LALStatus *stat, SkyPosition *output, SkyPosition *input, ConvertSkyParams *params)
void XLALNormalizeSkyPosition(double *restrict longitude, double *restrict latitude)
If sky-position is not in the canonical range , normalize it by mapping it into this coordinate-inter...
void LALNormalizeSkyPosition(LALStatus *stat, SkyPosition *posOut, const SkyPosition *posIn)
#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_GEOGRAPHIC
The Earth-fixed geographic coordinate system.
@ COORDINATESYSTEM_HORIZON
A horizon coordinate system.
@ COORDINATESYSTEM_EQUATORIAL
The sky-fixed equatorial coordinate system.
void LALEquatorialToGeographic(LALStatus *, SkyPosition *output, SkyPosition *input, LIGOTimeGPS *gpsTime)
void LALHorizonToSystem(LALStatus *, SkyPosition *output, SkyPosition *input, const SkyPosition *zenith)
void LALGeographicToEquatorial(LALStatus *, SkyPosition *output, SkyPosition *input, LIGOTimeGPS *gpsTime)
void LALSystemToHorizon(LALStatus *, SkyPosition *output, SkyPosition *input, const SkyPosition *zenith)
This structure stores parameters for the function LALConvertSkyPosition().
CoordinateSystem system
The coordinate system to which one is transforming.
SkyPosition * zenith
The position of the zenith of the horizon coordinate system; may be NULL if one is neither converting...
LIGOTimeGPS * gpsTime
The GPS time for conversions between Earth-fixed and sky-fixed coordinates; may be NULL if no such co...
LAL status structure, see The LALStatus structure for more details.
Definition: LALDatatypes.h:947
struct tagLALStatus * statusPtr
Pointer to the next node in the list; NULL if this function is not reporting a subroutine error.
Definition: LALDatatypes.h:954
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