CR2RE Pipeline Reference Manual 1.7.0
hdrl_barycorr.c
1/*
2 * This file is part of the HDRL
3 * Copyright (C) 2022 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/*-----------------------------------------------------------------------------
21 Includes
22 -----------------------------------------------------------------------------*/
23
24#ifdef HAVE_CONFIG_H
25#include <config.h>
26#endif
27
28#include "hdrl_barycorr.h"
29#include <erfa.h>
30
31/*----------------------------------------------------------------------------*/
46/*----------------------------------------------------------------------------*/
47
49
50/* ---------------------------------------------------------------------------*/
68/* ---------------------------------------------------------------------------*/
69
70cpl_error_code
71hdrl_eop_interpolate(double mjd, const cpl_table * eop_table, hdrl_parameter *
72 resample_par, double *pmx, double *pmy, double *dut)
73{
74 cpl_ensure_code (eop_table, CPL_ERROR_NULL_INPUT);
75 cpl_ensure_code (resample_par, CPL_ERROR_NULL_INPUT);
76 cpl_ensure_code (pmx, CPL_ERROR_NULL_INPUT);
77 cpl_ensure_code (pmy, CPL_ERROR_NULL_INPUT);
78 cpl_ensure_code (dut, CPL_ERROR_NULL_INPUT);
79
80 cpl_table * eop_table_loc = cpl_table_duplicate(eop_table);
81
82 /* Doing some baisc checks on the local eop table */
83
84 if (!cpl_table_has_column(eop_table_loc, "MJD") ||
85 !cpl_table_has_column(eop_table_loc, "PMX") ||
86 !cpl_table_has_column(eop_table_loc, "PMY") ||
87 !cpl_table_has_column(eop_table_loc, "DUT")){
88 cpl_table_delete(eop_table_loc);
89 return cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
90 "The EOP table does not have all required "
91 "columns, i.e. MJD, PMX, PMY, DUT");
92 }
93
94 cpl_table_unselect_all (eop_table_loc);
95 cpl_table_or_selected_invalid (eop_table_loc, "MJD");
96 cpl_table_or_selected_invalid (eop_table_loc, "PMX");
97 cpl_table_or_selected_invalid (eop_table_loc, "PMY");
98 cpl_table_or_selected_invalid (eop_table_loc, "DUT");
99 cpl_table_erase_selected (eop_table_loc);
100
101 if (cpl_table_get_nrow(eop_table_loc) < 1){
102 cpl_table_delete(eop_table_loc);
103 return cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
104 "The EOP table does not have entries after "
105 "removing all invalid elements");
106 }
107
108 cpl_size mjd_max = cpl_table_get_column_max(eop_table_loc, "MJD");
109 cpl_size mjd_min = cpl_table_get_column_min(eop_table_loc, "MJD");
110
111 if (mjd < mjd_min || mjd > mjd_max) {
112 *pmx = cpl_table_get_column_median(eop_table_loc, "PMX");
113 *pmy = cpl_table_get_column_median(eop_table_loc, "PMY");
114 *dut = cpl_table_get_column_median(eop_table_loc, "DUT");
115 cpl_msg_warning(cpl_func,"The exposure MJD-OBS is outside the validity "
116 "range of the EOP calibration. Using median values instead of "
117 "interpolated values. Please provide a more up to date EOP file for "
118 "a higher accuracy.");
119 cpl_table_delete(eop_table_loc);
120 return cpl_error_get_code();
121 }
122
123 hdrl_spectrum1D * pmx_spec =
124 hdrl_spectrum1D_convert_from_table(eop_table_loc, "PMX", "MJD", NULL, NULL,
125 hdrl_spectrum1D_wave_scale_linear);
126 hdrl_spectrum1D * pmy_spec =
127 hdrl_spectrum1D_convert_from_table(eop_table_loc, "PMY", "MJD", NULL, NULL,
128 hdrl_spectrum1D_wave_scale_linear);
129 hdrl_spectrum1D * dut_spec =
130 hdrl_spectrum1D_convert_from_table(eop_table_loc, "DUT", "MJD", NULL, NULL,
131 hdrl_spectrum1D_wave_scale_linear);
132
133 cpl_array * mjd_out = cpl_array_new(1, CPL_TYPE_DOUBLE);
134 cpl_array_set_double(mjd_out, 0, mjd);
135
136 hdrl_spectrum1D * pmx_spec_resampled =
137 hdrl_spectrum1D_resample_on_array(pmx_spec, mjd_out, resample_par);
138 hdrl_spectrum1D * pmy_spec_resampled =
139 hdrl_spectrum1D_resample_on_array(pmy_spec, mjd_out, resample_par);
140 hdrl_spectrum1D * dut_spec_resampled =
141 hdrl_spectrum1D_resample_on_array(dut_spec, mjd_out, resample_par);
142
143 cpl_array_delete(mjd_out);
144 hdrl_spectrum1D_delete(&pmx_spec);
145 hdrl_spectrum1D_delete(&pmy_spec);
146 hdrl_spectrum1D_delete(&dut_spec);
147
148 int rej_pmx = 0, rej_pmy = 0, rej_dut = 0;
149
150 hdrl_value hvpmx = hdrl_spectrum1D_get_flux_value(pmx_spec_resampled, 0, &rej_pmx);
151 hdrl_value hvpmy = hdrl_spectrum1D_get_flux_value(pmy_spec_resampled, 0, &rej_pmy);
152 hdrl_value hvdut = hdrl_spectrum1D_get_flux_value(dut_spec_resampled, 0, &rej_dut);
153
154 hdrl_spectrum1D_delete(&pmx_spec_resampled);
155 hdrl_spectrum1D_delete(&pmy_spec_resampled);
156 hdrl_spectrum1D_delete(&dut_spec_resampled);
157
158 *pmx = hvpmx.data;
159 *pmy = hvpmy.data;
160 *dut = hvdut.data;
161
162 cpl_table_delete(eop_table_loc);
163 return cpl_error_get_code();
164}
165
166/* ---------------------------------------------------------------------------*/
196/* ---------------------------------------------------------------------------*/
197cpl_error_code
198hdrl_barycorr_compute(double ra, double dec, const cpl_table * eop_table,
199 double mjdobs, double time_to_mid_exposure,
200 double longitude, double latitude, double elevation,
201 double pressure, double temperature,
202 double humidity, double wavelength , double *barycorr){
203
204 cpl_ensure_code(ra >= 0. && ra < 360., CPL_ERROR_ILLEGAL_INPUT);
205 cpl_ensure_code(dec >= -90. && dec <= 90., CPL_ERROR_ILLEGAL_INPUT);
206 /* Todo check range for longitude and latidude !*/
207 cpl_ensure_code(longitude >= -180. && longitude <= 180., CPL_ERROR_ILLEGAL_INPUT);
208 cpl_ensure_code(latitude >= -90. && latitude <= 90., CPL_ERROR_ILLEGAL_INPUT);
209 cpl_ensure_code(eop_table != NULL, CPL_ERROR_NULL_INPUT);
210
211
212 cpl_error_code err = CPL_ERROR_NONE;
213 ra *= CPL_MATH_RAD_DEG;
214 dec *= CPL_MATH_RAD_DEG;
215 longitude *= CPL_MATH_RAD_DEG; // Lon in [rad], East positive
216 latitude *= CPL_MATH_RAD_DEG; // Lat in [rad]
217
218 /* Mean mjdobs from the middle of the exposure time */
219 double mean_mjd = mjdobs + (time_to_mid_exposure/3600./24.);
220 cpl_msg_debug(cpl_func, "Mean MJD-OBS used to derive barycorr: %g", mean_mjd);
221
222
223 /* Compute Earth Orientation Parameters for the mean MJD */
224 double dut1 = 0, pmx = 0, pmy = 0;
225
226 // ToDo Expose it in the API?
227 // hdrl_parameter * resample_par = hdrl_spectrum1D_resample_interpolate_parameter_create(hdrl_spectrum1D_interp_akima);
228 // hdrl_parameter * resample_par = hdrl_spectrum1D_resample_interpolate_parameter_create(hdrl_spectrum1D_interp_cspline);
229
230 hdrl_parameter * resample_par =
231 hdrl_spectrum1D_resample_interpolate_parameter_create(hdrl_spectrum1D_interp_linear);
232 err =
233 hdrl_eop_interpolate(mean_mjd, eop_table, resample_par, &pmx, &pmy, &dut1) ;
234 hdrl_parameter_delete(resample_par);
235
236 /* Check output result */
237 if (err) {
238 return cpl_error_set_message(cpl_func, err, "Could not interpolate the "
239 "Earth Orientation Parameter table");
240 }
241
242 cpl_msg_debug(cpl_func, "Using the following Earth Orientation Parameter for "
243 "MJD-OBS %g: pmx: %g, pmy: %g, dut1: %g", mean_mjd, pmx, pmy, dut1);
244
245 pmx = pmx / 3600.0 * CPL_MATH_RAD_DEG;
246 pmy = pmy / 3600.0 * CPL_MATH_RAD_DEG;
247
248 cpl_msg_debug(cpl_func, "Input to the erfa function eraApco13():");
249 cpl_msg_indent_more();
250 cpl_msg_debug (cpl_func, "%12s %50s: %20.20g", "utc1 :", "UTC as a 2-part... ", 2400000.5);
251 cpl_msg_debug (cpl_func, "%12s %50s: %20.20g", "utc2 :", "...quasi Julian Date (Notes 1,2) ", mean_mjd);
252 cpl_msg_debug (cpl_func, "%12s %50s: %20.20g", "dut1 :", "UT1-UTC (seconds, Note 3) ", dut1);
253 cpl_msg_debug (cpl_func, "%12s %50s: %20.20g", "elong :", "longitude (radians, east +ve, Note 4) ", longitude);
254 cpl_msg_debug (cpl_func, "%12s %50s: %20.20g", "phi :", "latitude (geodetic, radians, Note 4) ", latitude);
255 cpl_msg_debug (cpl_func, "%12s %50s: %20.20g", "hm :", "height above ellipsoid (m, geodetic, Notes 4,6) ", elevation);
256 cpl_msg_debug (cpl_func, "%12s %50s: %20.20g", "xp :", "polar motion coordinates (radians, Note 5) ", pmx);
257 cpl_msg_debug (cpl_func, "%12s %50s: %20.20g", "yp :", "polar motion coordinates (radians, Note 5) ", pmy);
258 cpl_msg_debug (cpl_func, "%12s %50s: %20.20g", "phpa :", "pressure at the observer (hPa = mB, Note 6) ", pressure);
259 cpl_msg_debug (cpl_func, "%12s %50s: %20.20g", "tc :", "ambient temperature at the observer (deg C) ", temperature);
260 cpl_msg_debug (cpl_func, "%12s %50s: %20.20g", "rh :", "relative humidity at the observer (range 0-1) ", humidity);
261 cpl_msg_debug (cpl_func, "%12s %50s: %20.20g", "wl :", "wavelength (micrometers, Note 7) ", wavelength);
262 cpl_msg_indent_less();
263 /* Allocate memory for the tmp computations */
264 eraASTROM astrom;
265 double eo;
266 int eraApco13_status = 0;
267 cpl_msg_debug(cpl_func, "Calling erfa function eraApco13() ...");
268 eraApco13_status = eraApco13 (2400000.5, mean_mjd, dut1,
269 longitude, latitude, elevation,
270 pmx, pmy,
271 pressure, temperature, humidity, wavelength,
272 &astrom, &eo);
273 if (eraApco13_status < 0) {
274 *barycorr = NAN;
275 return cpl_error_set_message(cpl_func, CPL_ERROR_UNSPECIFIED,
276 "Erfa function eraApco13() did not succeed in"
277 " computing the barycentric correction");
278 }
279
280 cpl_msg_debug(cpl_func, "Output of the erfa function eraApco13():");
281 cpl_msg_indent_more();
282
283 cpl_msg_debug (cpl_func, " pmt : /* PM time interval (SSB, Julian years) */ : %20.20g ", astrom.pmt);
284 cpl_msg_debug (cpl_func, " eb[0] : /* SSB to observer (vector, au) */ : %20.20g ", astrom.eb[0]);
285 cpl_msg_debug (cpl_func, " eb[1] : /* SSB to observer (vector, au) */ : %20.20g ", astrom.eb[1]);
286 cpl_msg_debug (cpl_func, " eb[2] : /* SSB to observer (vector, au) */ : %20.20g ", astrom.eb[2]);
287 cpl_msg_debug (cpl_func, " eh[0] : /* Sun to observer (unit vector) */ : %20.20g ", astrom.eh[0]);
288 cpl_msg_debug (cpl_func, " eh[1] : /* Sun to observer (unit vector) */ : %20.20g ", astrom.eh[1]);
289 cpl_msg_debug (cpl_func, " eh[2] : /* Sun to observer (unit vector) */ : %20.20g ", astrom.eh[2]);
290 cpl_msg_debug (cpl_func, " em : /* distance from Sun to observer (au) */ : %20.20g ", astrom.em);
291 cpl_msg_debug (cpl_func, " v[0] : /* barycentric observer velocity (vector, c) */ : %20.20g ", astrom.v[0]);
292 cpl_msg_debug (cpl_func, " v[1] : /* barycentric observer velocity (vector, c) */ : %20.20g ", astrom.v[1]);
293 cpl_msg_debug (cpl_func, " v[2] : /* barycentric observer velocity (vector, c) */ : %20.20g ", astrom.v[2]);
294 cpl_msg_debug (cpl_func, " bm1 : /* sqrt(1-|v|^2): reciprocal of Lorenz factor */ : %20.20g ", astrom.bm1);
295 cpl_msg_debug (cpl_func, " bpn[0][0] : /* bias-precession-nutation matrix */ : %20.20g ", astrom.bpn[0][0]);
296 cpl_msg_debug (cpl_func, " along: : /* longitude + s' + dERA(DUT) (radians) */ : %20.20g ", astrom.along);
297 // phi may be un-initialized according to valgrind - not printing it at this stage
298 // cpl_msg_debug (cpl_func, " phi: : /* geodetic latitude (radians) */ : %20.20g ", astrom.phi);
299 cpl_msg_debug (cpl_func, " xpl : /* polar motion xp wrt local meridian (radians) */ : %20.20g ", astrom.xpl);
300 cpl_msg_debug (cpl_func, " ypl : /* polar motion yp wrt local meridian (radians) */ : %20.20g ", astrom.ypl);
301 cpl_msg_debug (cpl_func, " sphi : /* sine of geodetic latitude */ : %20.20g ", astrom.sphi);
302 cpl_msg_debug (cpl_func, " cphi : /* cosine of geodetic latitude */ : %20.20g ", astrom.cphi);
303 cpl_msg_debug (cpl_func, " diurab : /* magnitude of diurnal aberration vector */ : %20.20g ", astrom.diurab);
304 cpl_msg_debug (cpl_func, " eral : /* local Earth rotation angle (radians) */ : %20.20g ", astrom.eral);
305 cpl_msg_debug (cpl_func, " refa : /* refraction constant A (radians) */ : %20.20g ", astrom.refa);
306 cpl_msg_debug (cpl_func, " refb : /* refraction constant B (radians) */ : %20.20g ", astrom.refb);
307 cpl_msg_debug (cpl_func, " eo : /* equation of the origins (ERA-GST) */ : %20.20g ", astrom.refb);
308 cpl_msg_indent_less();
309
310
311 /* From UVES Pipeline:
312 * ... REFERENCE: THE ASTRONOMICAL ALMANAC 1982 PAGE:B17 */
313
314 *barycorr =
315 astrom.v[0]*cos(ra)*cos(dec)+
316 astrom.v[1]*sin(ra)*cos(dec)+
317 astrom.v[2]*sin(dec);
318
319 *barycorr *= CPL_PHYS_C; // m/s
320
321 return cpl_error_get_code();
322
323}
324
hdrl_barycorr_compute
cpl_error_code hdrl_barycorr_compute(double ra, double dec, const cpl_table *eop_table, double mjdobs, double time_to_mid_exposure, double longitude, double latitude, double elevation, double pressure, double temperature, double humidity, double wavelength, double *barycorr)
Derives the barycentric correction using the erfa function eraApco13(). The latter For a terrestrial ...
Definition hdrl_barycorr.c:198
hdrl_parameter_delete
void hdrl_parameter_delete(hdrl_parameter *obj)
shallow delete of a parameter
Definition hdrl_parameter.c:148
hdrl_spectrum1D_resample_interpolate_parameter_create
hdrl_parameter * hdrl_spectrum1D_resample_interpolate_parameter_create(const hdrl_spectrum1D_interpolation_method method)
constructor for the hdrl_parameter in the case of interpolation
Definition hdrl_spectrum_resample.c:240
hdrl_spectrum1D_delete
void hdrl_spectrum1D_delete(hdrl_spectrum1D **p_self)
hdrl_spectrum1D destructor
Definition hdrl_spectrum.c:322
hdrl_spectrum1D_convert_from_table
hdrl_spectrum1D * hdrl_spectrum1D_convert_from_table(const cpl_table *self, const char *flux_col_name, const char *wavelength_col_name, const char *flux_e_col_name, const char *flux_bpm_col_name, hdrl_spectrum1D_wave_scale scale)
convert a table to a spectrum
Definition hdrl_spectrum.c:1200
hdrl_spectrum1D_resample_on_array
hdrl_spectrum1D * hdrl_spectrum1D_resample_on_array(const hdrl_spectrum1D *self, const cpl_array *waves, const hdrl_parameter *par)
resample a hdrl_spectrum1D on the wavelengths contained in waves
Definition hdrl_spectrum_resample.c:484
hdrl_spectrum1D_get_flux_value
hdrl_value hdrl_spectrum1D_get_flux_value(const hdrl_spectrum1D *self, int idx, int *rej)
hdrl_spectrum1D getter for a flux value
Definition hdrl_spectrum.c:416
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2