CR2RE Pipeline Reference Manual 1.6.11
hdrl_dar.c
1/*
2 * This file is part of the HDRL
3 * Copyright (C) 2013 European Southern Observatory
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20#ifdef HAVE_CONFIG_H
21#include <config.h>
22#endif
23
24#if !defined(_POSIX_C_SOURCE) || (_POSIX_C_SOURCE - 0) < 200112L
25#define _POSIX_C_SOURCE 200112L
26#endif
27
28
29/*---------------------------------------------------------------------------*
30 * Includes *
31 *---------------------------------------------------------------------------*/
32#include "hdrl_dar.h"
33
34#include <math.h>
35
36/*---------------------------------------------------------------------------*
37 * Debugging/feature Macros *
38 *---------------------------------------------------------------------------*/
39
40
41/*---------------------------------------------------------------------------*
42 * Static *
43 *---------------------------------------------------------------------------*/
44
45
46/*---------------------------------------------------------------------------*/
53/*----------------------------------------------------------------------------*/
54
56
57
59
60/*-----------------------------------------------------------------------------
61 DAR Parameters Definition
62 -----------------------------------------------------------------------------*/
63typedef struct {
64 HDRL_PARAMETER_HEAD;
65 hdrl_value airmass; /* Air mass */
66 hdrl_value parang; /* Parallactic angle during exposure */
67 hdrl_value posang; /* Position angle on the sky from the angles we have */
68 hdrl_value temp; /* Temperature [Celsius] */
69 hdrl_value rhum; /* Relative humidity [%] */
70 hdrl_value pres; /* Pressure [mbar] */
71 cpl_wcs *wcs; /* World Coordinate system (WCS) in degrees(CDi_j) */
72} hdrl_dar_parameter;
73
74/* Parameter type */
75static hdrl_parameter_typeobj hdrl_dar_parameter_type = {
76 HDRL_PARAMETER_DAR, /* type */
77 (hdrl_alloc *)&cpl_malloc, /* fp_alloc */
78 (hdrl_free *)&cpl_free, /* fp_free */
79 NULL, /* fp_destroy */
80 sizeof(hdrl_dar_parameter), /* obj_size */
81};
82
83/*----------------------------------------------------------------------------*/
92 /*----------------------------------------------------------------------------*/
93cpl_error_code hdrl_dar_parameter_verify(const hdrl_parameter *param)
94{
95 cpl_error_ensure(param != NULL, CPL_ERROR_NULL_INPUT,
96 return CPL_ERROR_NULL_INPUT, "NULL Input Parameters");
97
98 cpl_error_ensure(hdrl_parameter_check_type(param, &hdrl_dar_parameter_type), CPL_ERROR_ILLEGAL_INPUT,
99 return CPL_ERROR_ILLEGAL_INPUT, "Expected DAR parameter");
100
101 const hdrl_dar_parameter *p_loc = (const hdrl_dar_parameter *)param;
102 hdrl_value airmass = p_loc->airmass;
103 hdrl_value parang = p_loc->parang; /* Degrees */
104 hdrl_value posang = p_loc->posang; /* Degrees */
105 hdrl_value temp = p_loc->temp; /* T [Celsius] */
106 hdrl_value rhum = p_loc->rhum; /* Relative humidity [%] */
107 hdrl_value pres = p_loc->pres; /* Pressure [mbar] */
108
109 cpl_error_ensure(airmass.data >= 0. && airmass.error >= 0., CPL_ERROR_ILLEGAL_INPUT,
110 return CPL_ERROR_ILLEGAL_INPUT, "Airmass parameter not valid");
111
112 cpl_error_ensure(parang.data >= -180. && parang.data <= 180. && parang.error >= 0., CPL_ERROR_ILLEGAL_INPUT,
113 return CPL_ERROR_ILLEGAL_INPUT, "Paralactic angle not valid");
114
115 cpl_error_ensure(posang.data >= -360. && posang.data <= 360. && posang.error >= 0., CPL_ERROR_ILLEGAL_INPUT,
116 return CPL_ERROR_ILLEGAL_INPUT, "Position angle not valid");
117
118 cpl_error_ensure(temp.data >= -273.15 && temp.error >= 0., CPL_ERROR_ILLEGAL_INPUT,
119 return CPL_ERROR_ILLEGAL_INPUT, "Temperature not valid");
120
121 cpl_error_ensure(rhum.data >= 0. && rhum.data <=100 && rhum.error >= 0., CPL_ERROR_ILLEGAL_INPUT,
122 return CPL_ERROR_ILLEGAL_INPUT, "Humidity percent value not valid");
123
124 cpl_error_ensure( pres.data >= 0. && pres.error >= 0., CPL_ERROR_ILLEGAL_INPUT,
125 return CPL_ERROR_ILLEGAL_INPUT, "Pressure not valid");
126
127 cpl_wcs* wcs = p_loc->wcs; /* Degrees */
128
129 cpl_error_ensure(wcs != NULL, CPL_ERROR_NULL_INPUT,
130 return CPL_ERROR_NULL_INPUT, "NULL WCS Input");
131
132 return CPL_ERROR_NONE;
133}
134
136
137/* ---------------------------------------------------------------------------*/
158/* ---------------------------------------------------------------------------*/
159hdrl_parameter * hdrl_dar_parameter_create(hdrl_value airmass, hdrl_value parang,
160 hdrl_value posang, hdrl_value temp, hdrl_value rhum, hdrl_value pres, cpl_wcs *wcs)
161{
162 hdrl_dar_parameter *p = (hdrl_dar_parameter *)hdrl_parameter_new(&hdrl_dar_parameter_type);
163
164 p->airmass = airmass;
165 p->parang = parang;
166 p->posang = posang;
167 p->temp = temp;
168 p->rhum = rhum;
169 p->pres = pres;
170
171 p->wcs = wcs;
172
173 /* After add the parameters, verify if they are correct */
174 if (hdrl_dar_parameter_verify((hdrl_parameter *)p) != CPL_ERROR_NONE) {
175 hdrl_parameter_delete((hdrl_parameter*)p);
176 return NULL;
177 }
178
179 return (hdrl_parameter *)p;
180}
181
182
183/*----------------------------------------------------------------------------*/
217/*----------------------------------------------------------------------------*/
218cpl_error_code hdrl_dar_compute(const hdrl_parameter *params,
219 const hdrl_value lambdaRef, const cpl_vector *lambdaIn,
220 cpl_vector *xShift, cpl_vector *yShift, cpl_vector *xShiftErr, cpl_vector *yShiftErr)
221{
222
223 cpl_error_ensure(params && lambdaIn && xShift && yShift != NULL, CPL_ERROR_NULL_INPUT,
224 return CPL_ERROR_NULL_INPUT, "NULL Input Parameters");
225
226 /* Check entry values in hdrl_parameter */
227 if (hdrl_dar_parameter_verify(params) != CPL_ERROR_NONE) {
228 return CPL_ERROR_UNSPECIFIED;
229 }
230
231 cpl_error_ensure(lambdaRef.data >= 0, CPL_ERROR_ILLEGAL_INPUT,
232 return CPL_ERROR_ILLEGAL_INPUT, "Reference wavelength must be >=0");
233
234 /* Local Usage Parameters */
235 const hdrl_dar_parameter *p_loc = (const hdrl_dar_parameter *)params;
236 hdrl_value airmass = p_loc->airmass;
237 hdrl_value parang = p_loc->parang; /* Degrees */
238 hdrl_value posang = p_loc->posang; /* Degrees */
239 hdrl_value temp = p_loc->temp; /* T [Celsius] */
240 hdrl_value rhum = p_loc->rhum; /* Relative humidity [%] */
241 hdrl_value pres = p_loc->pres; /* Pressure [mbar] */
242 cpl_wcs *wcs = p_loc->wcs; /* Degrees */
243
244 /* Check if the airmass is at less 1. */
245 cpl_ensure_code(airmass.data >= 1., cpl_error_get_code());
246
247 /* simple zenith distance in radians */
248 hdrl_value z = {acos(1. / airmass.data),
249 airmass.error * fabs( (-1. / pow(airmass.data, 2)) / sqrt(1. - pow(1. / airmass.data, 2)) )};
250
251 /* ----------------------------------------------------------------- *
252 * Compute the refractive index at lambdaRef with FILIPPENKO method *
253 * in um and output properties in "natural" (for the formulae) units *
254 * ----------------------------------------------------------------- */
255
256 /* Calculate temperature and error in Kelvin */
257 double temp_kel_data = temp.data + 273.15;
258 double temp_kel_err = (temp.error / fabs(temp.data)) * fabs(temp_kel_data);
259 hdrl_value temp_kel = {temp_kel_data, temp_kel_err};
260
261 /* Use the Owens formula to derive saturation pressure. Needs T[K] */
262 hdrl_value sp = hdrl_dar_owens_saturation_pressure(temp_kel);
263
264 /* Conversion from hPa (or mbar) to mmHg, needed for Filippenko *
265 * using that, derive the water vapor pressure in mmHg */
266 double HDRL_PHYS_hPa_TO_mmHg = 0.75006158;
267
268 /* Convert relative humidity [%] to fraction */
269 rhum.data /= 100.;
270 rhum.error /= 100.;
271
272 /* water vapor pressure in mmHg */
273 hdrl_value fp = {rhum.data *sp.data *HDRL_PHYS_hPa_TO_mmHg,
274 fabs(HDRL_PHYS_hPa_TO_mmHg * sp.data ) * rhum.error
275 + fabs(HDRL_PHYS_hPa_TO_mmHg * rhum.data) * sp.error};
276
277 /* need the pressure in mmHg */
278 pres.data *= HDRL_PHYS_hPa_TO_mmHg;
279 pres.error *= HDRL_PHYS_hPa_TO_mmHg;
280
281 /* refractive index of air at reference wavelength. Needs lambda[um] */
282 hdrl_value lambdaRef_um = {lambdaRef.data * 1e-4,lambdaRef.error * 1e-4};
283 hdrl_value nr0 = hdrl_dar_filippenko_refractive_index(lambdaRef_um, pres, temp, fp);
284
285 /* Obtain shift with scale: Absolute Shift for a lambdaRef, xshift is in *
286 * E-W direction for posang = 0, yshift is N-S Shift units --> Degrees */
287 hdrl_value x_shift = {-sin( (parang.data + posang.data) * CPL_MATH_RAD_DEG),
288 parang.error * fabs(-CPL_MATH_RAD_DEG * cos(parang.data + posang.data))
289 + posang.error * fabs(-CPL_MATH_RAD_DEG * cos(parang.data + posang.data))};
290
291 hdrl_value y_shift = { cos( (parang.data + posang.data) * CPL_MATH_RAD_DEG),
292 parang.error * fabs(-CPL_MATH_RAD_DEG * sin(parang.data + posang.data))
293 + posang.error * fabs(-CPL_MATH_RAD_DEG * sin(parang.data + posang.data))};
294
295 /* Get scale in the world cordinate system (wcs) and apply them */
296 double xscale, yscale;
297 hdrl_dar_wcs_get_scales(wcs, &xscale, &yscale);
298
299 x_shift.data /= xscale;
300 x_shift.error /= xscale;
301
302 y_shift.data /= yscale;
303 y_shift.error /= yscale;
304
305 /* Diff. refr. base in arcsec converted from radians (Filippenko does *
306 * the conversion using x206265 which converts radians to arcsec) */
307 hdrl_value dr0 = {tan(z.data) * CPL_MATH_DEG_RAD,
308 z.error * fabs( (1. + pow(tan(z.data), 2)) * CPL_MATH_DEG_RAD)};
309
310 /* ------------------------------------------------------------------ *
311 * Calculate the relative lambda of in array (in) *
312 * apply the absolute shift (x_shift, y_shift) for lamdaRef *
313 * for obtain the out arrays (xShift, yShift ) *
314 * ------------------------------------------------------------------ */
315 cpl_size i;
316 cpl_size nmax = cpl_vector_get_size(lambdaIn);
317
318 HDRL_OMP(omp parallel for \
319 default(none) \
320 shared( nmax, lambdaIn, \
321 xShift, yShift, xShiftErr, yShiftErr, \
322 lambdaRef_um, pres, temp, fp, dr0, nr0, \
323 x_shift, y_shift ))
324 for (i = 0; i < nmax; i++) {
325
326 double lambda = cpl_vector_get(lambdaIn, i);
327 if (isfinite(lambda) != 0) {
328
329 hdrl_value lambda_um = {lambda * 1e-4, lambdaRef_um.error};
330 hdrl_value nr = hdrl_dar_filippenko_refractive_index(lambda_um, pres, temp, fp);
331
332 hdrl_value dr = {dr0.data * (nr0.data - nr.data),
333 dr0.error * fabs(nr0.data - nr.data)
334 + nr0.error * fabs( dr0.data)
335 + nr.error * fabs(-dr0.data)};
336
337 hdrl_value shiftPlaneX = {x_shift.data * dr.data,
338 x_shift.error * fabs(dr.data)
339 + dr.error * fabs(x_shift.data)};
340 cpl_vector_set(xShift, i, shiftPlaneX.data );
341 cpl_vector_set(xShiftErr, i, shiftPlaneX.error);
342
343 hdrl_value shiftPlaneY = {y_shift.data * dr.data,
344 y_shift.error * fabs(dr.data)
345 + dr.error * fabs(y_shift.data)};
346 cpl_vector_set(yShift, i, shiftPlaneY.data );
347 cpl_vector_set(yShiftErr, i, shiftPlaneY.error);
348
349 } else {
350
351 cpl_vector_set(xShift, i, NAN);
352 cpl_vector_set(xShiftErr, i, NAN);
353
354 cpl_vector_set(yShift, i, NAN);
355 cpl_vector_set(yShiftErr, i, NAN);
356 }
357 }
358
359 return CPL_ERROR_NONE;
360}
361
362
363/*----------------------------------------------------------------------------*/
379/*----------------------------------------------------------------------------*/
380hdrl_value hdrl_dar_owens_saturation_pressure(hdrl_value hvT)
381{
382 double T = hvT.data;
383 double errorT = hvT.error;
384
385 return (hdrl_value){-10474.0 + 116.43 * T - 0.43284 * T * T + 0.00053840 * pow(T, 3),
386 errorT * fabs(0.0016152 * T * T - 0.86568 * T + 116.43)};
387}
388
389/*----------------------------------------------------------------------------*/
425/*----------------------------------------------------------------------------*/
426hdrl_value hdrl_dar_filippenko_refractive_index(
427 hdrl_value hvL, hdrl_value hvP, hdrl_value hvT, hdrl_value hvF)
428{
429 double l = hvL.data,
430 P = hvP.data,
431 T = hvT.data,
432 f = hvF.data;
433
434 double errorL = hvL.error,
435 errorP = hvP.error,
436 errorT = hvT.error,
437 errorF = hvF.error;
438
439 /* inverse square of the wavelength */
440 double lisq = 1. / (l * l);
441 double errorLisq = errorL * fabs(-2. / pow(l, 3) );
442
443 /* 10^6 [n(lambda) - 1] at standard environmental conditions, Eq. (1) */
444 double nl1 = 64.328 + 29498.1 / (146. - lisq) + 255.4 / (41. - lisq);
445 double errorNl1 = errorLisq * fabs( 29498.1 / pow(146. - lisq, 2) + 255.4 / pow(41. - lisq, 2) );
446
447 /* correction for non-standard conditions, Eq. (2) */
448 double factor = 1.e-6;
449 double nl2A = nl1 * ( P / 720.883 * (1. + (1.049 -0.0157 * T) * 1e-6 * P) / (1. + 0.003661 * T) );
450 double errorNl2A1 = errorNl1 * fabs( factor *(P / 720.883 * (1. + (1.049 - 0.0157 * T) * 1e-6 * P) / (1. + 0.003661 * T) ) );
451 double errorNl2A2 = errorP * fabs( factor *(nl1 / (720.883 * (1. + 0.003661 * T)) *( (1. + (1.049 - 0.0157 * T) * 1e-6 * P) + P * (1.049 - 0.0157 * T) * 1e-6) ) );
452 double errorNl2A3 = errorT * fabs( factor *(nl1 * P / 720.883 * ( ( -0.0157 * 1e-6 * P * (1. + 0.003661 * T) - 0.003661 * (1. + (1.049 - 0.0157 * T) * 1e-6 * P) )/pow(1. + 0.003661 * T, 2) ) ) );
453 double errorNl2A = errorNl2A1 + errorNl2A2 + errorNl2A3;
454
455 /* Calcule correction for water vapor, Eq. (3) */
456 double nl2B = (0.0624 - 0.000680 * lisq) / (1. + 0.003661 * T) * f;
457 double errorNl2B1 = errorLisq * fabs( -0.000680 * f / (1. + 0.003661 * T) );
458 double errorNl2B2 = errorT * fabs( -0.003661 * (0.0624 - 0.000680 * lisq) * f / pow(1. + 0.003661 * T, 2) );
459 double errorNl2B3 = errorF * fabs( (0.0624 - 0.000680 * lisq) / (1. + 0.003661 * T) );
460 double errorNl2B = errorNl2B1 + errorNl2B2 + errorNl2B3;
461
462 /* Apply correction for water vapor, Eq. (3) */
463 double nl2 = nl2A - nl2B;
464 double errorNl2 = errorNl2A + errorNl2B;
465
466 /* convert to refractive index n(lambda) */
467 double nl3 = nl2 * 1e-6 + 1.;
468 double errorNl3 = fabs(errorNl2 * 1e-6);
469
470 return (hdrl_value){nl3,errorNl3};
471}
472
473
475
476/*----------------------------------------------------------------------------*/
494/*----------------------------------------------------------------------------*/
495cpl_error_code hdrl_dar_wcs_get_scales(
496 cpl_wcs *wcs,
497 double *aXScale,
498 double *aYScale)
499{
500 cpl_ensure_code(aXScale && aYScale, CPL_ERROR_NULL_INPUT);
501
502 cpl_errorstate prestate = cpl_errorstate_get();
503
504 const cpl_matrix *cd = cpl_wcs_get_cd(wcs);
505
506 /* take the absolute and scale by 3600 to get positive arcseconds */
507 double cd11 = cpl_matrix_get(cd, 0, 0),
508 cd12 = cpl_matrix_get(cd, 0, 1),
509 cd21 = cpl_matrix_get(cd, 1, 0),
510 cd22 = cpl_matrix_get(cd, 1, 1);
511
512 double det = cd11 * cd22 - cd12 * cd21;
513 cpl_ensure_code(cpl_errorstate_is_equal(prestate), cpl_error_get_code());
514
515 if (det < 0.) {
516 cd12 *= -1.;
517 cd11 *= -1.;
518 }
519
520 /* matrix without rotation */
521 if (cd12 == 0. && cd21 == 0.) {
522 *aXScale = cd11;
523 *aYScale = cd22;
524 return CPL_ERROR_NONE;
525 }
526
527 /* Only the absolute value */
528 *aXScale = sqrt(cd11 * cd11 + cd12 * cd12);
529 *aYScale = sqrt(cd22 * cd22 + cd21 * cd21);
530
531 return CPL_ERROR_NONE;
532}
533
535
hdrl_dar_compute
cpl_error_code hdrl_dar_compute(const hdrl_parameter *params, const hdrl_value lambdaRef, const cpl_vector *lambdaIn, cpl_vector *xShift, cpl_vector *yShift, cpl_vector *xShiftErr, cpl_vector *yShiftErr)
Correct the pixel coordinates of all pixels of a given pixel table for differential atmospheric refra...
Definition hdrl_dar.c:218
hdrl_dar_filippenko_refractive_index
hdrl_value hdrl_dar_filippenko_refractive_index(hdrl_value hvL, hdrl_value hvP, hdrl_value hvT, hdrl_value hvF)
Compute the refractive index for the given wavelength following Filippenko formulae....
Definition hdrl_dar.c:426
hdrl_dar_parameter_create
hdrl_parameter * hdrl_dar_parameter_create(hdrl_value airmass, hdrl_value parang, hdrl_value posang, hdrl_value temp, hdrl_value rhum, hdrl_value pres, cpl_wcs *wcs)
Creates DAR parameters object with the values in the header.
Definition hdrl_dar.c:159
hdrl_dar_owens_saturation_pressure
hdrl_value hdrl_dar_owens_saturation_pressure(hdrl_value hvT)
Compute the saturation pressure using the Owens calibration.
Definition hdrl_dar.c:380
hdrl_parameter_delete
void hdrl_parameter_delete(hdrl_parameter *obj)
shallow delete of a parameter
Definition hdrl_parameter.c:148
Generated on Tue Jun 16 2026 00:11:35 for CR2RE Pipeline Reference Manual by doxygen 1.13.2