CR2RE Pipeline Reference Manual 1.7.0
irplib_wcs.c
1/* $Id: irplib_wcs.c,v 1.8 2010-10-07 14:10:55 llundin Exp $
2 *
3 * This file is part of the irplib package
4 * Copyright (C) 2002,2003 European Southern Observatory
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02111-1307 USA
19 */
20
21/*
22 * $Author: llundin $
23 * $Date: 2010-10-07 14:10:55 $
24 * $Revision: 1.8 $
25 * $Name: not supported by cvs2svn $
26 */
27
28/*-----------------------------------------------------------------------------
29 Includes
30 -----------------------------------------------------------------------------*/
31
32#ifdef HAVE_CONFIG_H
33#include <config.h>
34#endif
35
36#include "irplib_wcs.h"
37
38#include <math.h>
39
40/*----------------------------------------------------------------------------*/
44/*----------------------------------------------------------------------------*/
45
46static cpl_error_code irplib_wcs_is_iso8601(int, int, int, int, int, double);
47
48/* Standard year-2000 form: CCYY-MM-DD[Thh:mm:ss[.sss...]] */
49#define IRPLIB_ISO8601_FORMAT "%4d-%2d-%2dT%2d:%2d:%lf"
50
52
53/*----------------------------------------------------------------------------*/
64/*----------------------------------------------------------------------------*/
65cpl_error_code
66irplib_wcs_xytoradec(const cpl_wcs *wcs,
67 double x,
68 double y,
69 double *ra,
70 double *dec)
71{
72 cpl_matrix *xy;
73 cpl_matrix *radec = NULL;
74 cpl_array *status = NULL;
75 cpl_error_code error;
76
77 cpl_ensure_code(ra != NULL, CPL_ERROR_NULL_INPUT);
78 cpl_ensure_code(dec != NULL, CPL_ERROR_NULL_INPUT);
79
80 /* Load up the information */
81 xy = cpl_matrix_new(1, 2);
82 cpl_matrix_set(xy, 0, 0, x);
83 cpl_matrix_set(xy, 0, 1, y);
84
85 /* Call the conversion routine */
86 error = cpl_wcs_convert(wcs, xy, &radec, &status, CPL_WCS_PHYS2WORLD);
87
88 cpl_matrix_delete(xy);
89
90 if (!error) {
91 /* Pass it back now */
92 *ra = cpl_matrix_get(radec, 0, 0);
93 *dec = cpl_matrix_get(radec, 0, 1);
94 }
95
96 /* Tidy and propagate error, if any */
97 cpl_matrix_delete(radec);
98 cpl_array_delete(status);
99
100 return cpl_error_set_where(cpl_func);
101}
102
103/*----------------------------------------------------------------------------*/
114/*----------------------------------------------------------------------------*/
115cpl_error_code
116irplib_wcs_radectoxy(const cpl_wcs *wcs,
117 double ra,
118 double dec,
119 double *x,
120 double *y)
121{
122 cpl_matrix *radec;
123 cpl_matrix *xy = NULL;
124 cpl_array *status = NULL;
125 cpl_error_code error;
126
127 cpl_ensure_code(x != NULL, CPL_ERROR_NULL_INPUT);
128 cpl_ensure_code(y != NULL, CPL_ERROR_NULL_INPUT);
129
130 /* Feed the matrix with RA, DEC */
131 radec = cpl_matrix_new(1, 2);
132 cpl_matrix_set(radec, 0, 0, ra);
133 cpl_matrix_set(radec, 0, 1, dec);
134
135 error = cpl_wcs_convert(wcs, radec, &xy, &status, CPL_WCS_WORLD2PHYS);
136
137 cpl_matrix_delete(radec);
138
139 if (!error) {
140 *x = cpl_matrix_get(xy, 0, 0);
141 *y = cpl_matrix_get(xy, 0, 1);
142 }
143
144 /* Tidy and propagate error, if any */
145 cpl_array_delete(status);
146 cpl_matrix_delete(xy);
147
148 return cpl_error_set_where(cpl_func);
149}
150
151/*----------------------------------------------------------------------------*/
161/*----------------------------------------------------------------------------*/
162double
163irplib_wcs_great_circle_dist(double ra1, double dec1, double ra2, double dec2)
164{
165 /* Convert all input from degrees to radian - and back for the result */
166 const double dra = sin(CPL_MATH_RAD_DEG * (ra2 - ra1) / 2.0);
167 const double ddec = sin(CPL_MATH_RAD_DEG * (dec2 - dec1) / 2.0);
168
169 dec1 *= CPL_MATH_RAD_DEG;
170 dec2 *= CPL_MATH_RAD_DEG;
171
172 return 2.0 * asin(sqrt(ddec * ddec + cos(dec1) * cos(dec2) * dra * dra)) *
173 CPL_MATH_DEG_RAD;
174}
175
176
177/*----------------------------------------------------------------------------*/
191/*----------------------------------------------------------------------------*/
192cpl_error_code
193irplib_wcs_mjd_from_iso8601(double *pmjd,
194 int year,
195 int month,
196 int day,
197 int hour,
198 int minute,
199 double second)
200{
201 cpl_ensure_code(pmjd != NULL, CPL_ERROR_NULL_INPUT);
202 cpl_ensure_code(!irplib_wcs_is_iso8601(year, month, day, hour, minute,
203 second),
204 cpl_error_get_code());
205
206 /* Compute MJD. */
207 *pmjd = (double)((1461 * (year - (12 - month) / 10 + 4712)) / 4 +
208 (306 * ((month + 9) % 12) + 5) / 10 -
209 (3 * ((year - (12 - month) / 10 + 4900) / 100)) / 4 + day -
210 2399904) +
211 (hour + (minute + second / 60.0) / 60.0) / 24.0;
212
213 return CPL_ERROR_NONE;
214}
215
216
217/*----------------------------------------------------------------------------*/
233/*----------------------------------------------------------------------------*/
234cpl_error_code
235irplib_wcs_iso8601_from_string(int *pyear,
236 int *pmonth,
237 int *pday,
238 int *phour,
239 int *pminute,
240 double *psecond,
241 const char *iso8601)
242{
243 int nret;
244
245 cpl_ensure_code(pyear != NULL, CPL_ERROR_NULL_INPUT);
246 cpl_ensure_code(pmonth != NULL, CPL_ERROR_NULL_INPUT);
247 cpl_ensure_code(pday != NULL, CPL_ERROR_NULL_INPUT);
248 cpl_ensure_code(phour != NULL, CPL_ERROR_NULL_INPUT);
249 cpl_ensure_code(pminute != NULL, CPL_ERROR_NULL_INPUT);
250 cpl_ensure_code(psecond != NULL, CPL_ERROR_NULL_INPUT);
251 cpl_ensure_code(iso8601 != NULL, CPL_ERROR_NULL_INPUT);
252
253 nret = sscanf(iso8601, IRPLIB_ISO8601_FORMAT, pyear, pmonth, pday, phour,
254 pminute, psecond);
255
256 if (nret != 6) {
257 return cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
258 "Parsed"
259 " %d != 6: input %s is not in format %s",
260 nret, iso8601, IRPLIB_ISO8601_FORMAT);
261 }
262
263 return irplib_wcs_is_iso8601(*pyear, *pmonth, *pday, *phour, *pminute,
264 *psecond)
265 ? cpl_error_set_where(cpl_func)
266 : CPL_ERROR_NONE;
267}
268
269
270/*----------------------------------------------------------------------------*/
279/*----------------------------------------------------------------------------*/
280cpl_error_code
281irplib_wcs_mjd_from_string(double *pmjd, const char *iso8601)
282{
283 int year, day, month, hour, minute;
284 double second;
285
286 return irplib_wcs_iso8601_from_string(&year, &month, &day, &hour, &minute,
287 &second, iso8601) ||
288 irplib_wcs_mjd_from_iso8601(pmjd, year, month, day, hour,
289 minute, second)
290 ? cpl_error_set_where(cpl_func)
291 : CPL_ERROR_NONE;
292}
293
294
295/*----------------------------------------------------------------------------*/
309/*----------------------------------------------------------------------------*/
310cpl_error_code
311irplib_wcs_iso8601_from_mjd(int *pyear,
312 int *pmonth,
313 int *pday,
314 int *phour,
315 int *pminute,
316 double *psecond,
317 double mjd)
318{
319 int jd, n4, dd;
320 double t;
321
322 cpl_ensure_code(pyear != NULL, CPL_ERROR_NULL_INPUT);
323 cpl_ensure_code(pmonth != NULL, CPL_ERROR_NULL_INPUT);
324 cpl_ensure_code(pday != NULL, CPL_ERROR_NULL_INPUT);
325 cpl_ensure_code(phour != NULL, CPL_ERROR_NULL_INPUT);
326 cpl_ensure_code(pminute != NULL, CPL_ERROR_NULL_INPUT);
327 cpl_ensure_code(psecond != NULL, CPL_ERROR_NULL_INPUT);
328
329 /* Copied from datfix() in wcslib (v. 4.4.4) */
330
331 jd = 2400001 + (int)mjd;
332
333 n4 = 4 * (jd + ((2 * ((4 * jd - 17918) / 146097) * 3) / 4 + 1) / 2 - 37);
334 dd = 10 * (((n4 - 237) % 1461) / 4) + 5;
335
336 *pyear = n4 / 1461 - 4712;
337 *pmonth = (2 + dd / 306) % 12 + 1;
338 *pday = (dd % 306) / 10 + 1;
339
340 t = mjd - (int)mjd; /* t is now days */
341
342 t *= 24.0; /* t is now hours */
343 *phour = (int)t;
344 t = 60.0 * (t - *phour); /* t is now minutes */
345 *pminute = (int)t;
346 *psecond = 60.0 * (t - *pminute);
347
348 /* A failure here implies that this code has a bug */
349 cpl_ensure_code(!irplib_wcs_is_iso8601(*pyear, *pmonth, *pday, *phour,
350 *pminute, *psecond),
351 CPL_ERROR_UNSPECIFIED);
352
353 return CPL_ERROR_NONE;
354}
355
356
358
359/*----------------------------------------------------------------------------*/
371/*----------------------------------------------------------------------------*/
372static cpl_error_code
373irplib_wcs_is_iso8601(int year,
374 int month,
375 int day,
376 int hour,
377 int minute,
378 double second)
379{
380 const cpl_boolean is_leap = (year % 4) ? CPL_FALSE : CPL_TRUE;
381 const int mlen[] = {
382 0, 31, is_leap ? 29 : 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
383
384 cpl_ensure_code(month > 0, CPL_ERROR_ILLEGAL_INPUT);
385 cpl_ensure_code(month <= 12, CPL_ERROR_ILLEGAL_INPUT);
386
387 cpl_ensure_code(day > 0, CPL_ERROR_ILLEGAL_INPUT);
388 cpl_ensure_code(day <= mlen[month], CPL_ERROR_ILLEGAL_INPUT);
389
390 cpl_ensure_code(minute < 60, CPL_ERROR_ILLEGAL_INPUT);
391 cpl_ensure_code(minute >= 0, CPL_ERROR_ILLEGAL_INPUT);
392
393 cpl_ensure_code(second < 60.0, CPL_ERROR_ILLEGAL_INPUT);
394 cpl_ensure_code(second >= 0.0, CPL_ERROR_ILLEGAL_INPUT);
395
396 cpl_ensure_code(hour >= 0, CPL_ERROR_ILLEGAL_INPUT);
397 /* 24:00:00 is valid ISO-8601 */
398 cpl_ensure_code(hour <= (minute > 0 || second > 0.0 ? 23 : 24),
399 CPL_ERROR_ILLEGAL_INPUT);
400
401 return CPL_ERROR_NONE;
402}
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2