LAL  7.1.7.1-2d066e5
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 */
72 void
74  SkyPosition *output,
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 ) &&
122  ( params->zenith->system == COORDINATESYSTEM_EQUATORIAL ) )
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  */
165 void
166 LALNormalizeSkyPosition (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  */
194 void
195 XLALNormalizeSkyPosition ( 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 /** @} */
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
void LALEclipticToEquatorial(LALStatus *stat, SkyPosition *output, SkyPosition *input)
The galactic coordinate system.
This structure stores parameters for the function LALConvertSkyPosition().
SkyPosition * zenith
The position of the zenith of the horizon coordinate system; may be NULL if one is neither converting...
void LALEquatorialToGeographic(LALStatus *, SkyPosition *output, SkyPosition *input, LIGOTimeGPS *gpsTime)
The sky-fixed equatorial coordinate system.
REAL8 latitude
The latitudinal coordinate (in radians), as defined above.
void LALEquatorialToEcliptic(LALStatus *stat, SkyPosition *output, SkyPosition *input)
CoordinateSystem system
The coordinate system to which one is transforming.
void LALConvertSkyCoordinates(LALStatus *stat, SkyPosition *output, SkyPosition *input, ConvertSkyParams *params)
void LALGalacticToEquatorial(LALStatus *stat, SkyPosition *output, SkyPosition *input)
#define SKYCOORDINATESH_ESYS
Wrong coordinate system in input.
void LALNormalizeSkyPosition(LALStatus *stat, SkyPosition *posOut, const SkyPosition *posIn)
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 *, SkyPosition *output, SkyPosition *input, LIGOTimeGPS *gpsTime)
#define ATTATCHSTATUSPTR(statusptr)
#define RETURN(statusptr)
The ecliptic coordinate system.
#define INITSTATUS(statusptr)
The Earth-fixed geographic coordinate system.
void LALEquatorialToGalactic(LALStatus *stat, SkyPosition *output, SkyPosition *input)
#define DETATCHSTATUSPTR(statusptr)
#define ABORT(statusptr, code, mesg)
A horizon coordinate system.
LAL status structure, see The LALStatus structure for more details.
Definition: LALDatatypes.h:947
void LALHorizonToSystem(LALStatus *, SkyPosition *output, SkyPosition *input, const SkyPosition *zenith)
REAL8 longitude
The longitudinal coordinate (in radians), as defined above.
#define SKYCOORDINATESH_ENUL
Unexpected null pointer in arguments.
LIGOTimeGPS * gpsTime
The GPS time for conversions between Earth-fixed and sky-fixed coordinates; may be NULL if no such co...
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...
#define TRY(func, statusptr)
#define LAL_PI
Archimedes&#39;s constant, pi.
Definition: LALConstants.h:179
void LALSystemToHorizon(LALStatus *, 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