1 /*********************************************************************
2 Units -- Convert data from one unit to other.
3 This is part of GNU Astronomy Utilities (Gnuastro) package.
4
5 Original author:
6 Kartik Ohri <kartikohri13@gmail.com>
7 Contributing author(s):
8 Mohammad Akhlaghi <mohammad@akhlaghi.org>
9 Copyright (C) 2020-2021, Free Software Foundation, Inc.
10
11 Gnuastro is free software: you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation, either version 3 of the License, or (at your
14 option) any later version.
15
16 Gnuastro is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 General Public License for more details.
20
21 You should have received a copy of the GNU General Public License
22 along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
23 **********************************************************************/
24 #include <config.h>
25
26 #include <math.h>
27 #include <errno.h>
28 #include <error.h>
29 #include <stdio.h>
30 #include <stdlib.h>
31 #include <string.h>
32
33 #include <gnuastro/type.h>
34 #include <gnuastro/pointer.h>
35
36
37
38 /**********************************************************************/
39 /**************** Functions to parse strings *****************/
40 /**********************************************************************/
41 /* Parse the input string consisting of numbers separated by given
42 delimiter into an array. */
43 int
gal_units_extract_decimal(char * convert,const char * delimiter,double * args,size_t n)44 gal_units_extract_decimal(char *convert, const char *delimiter,
45 double *args, size_t n)
46 {
47 size_t i = 0;
48 char *copy, *token, *end;
49
50 /* Create a copy of the string to be parsed and parse it. This is because
51 it will be modified during the parsing. */
52 copy=strdup(convert);
53 do
54 {
55 /* Check if the required number of arguments are passed */
56 if(i==n+1)
57 {
58 free(copy);
59 error(0, 0, "%s: input '%s' exceeds maximum number of arguments "
60 "(%zu)", __func__, convert, n);
61 return 0;
62 }
63
64 /* Extract the substring till the next delimiter */
65 token=strtok(i==0?copy:NULL, delimiter);
66 if(token)
67 {
68 /* Parse extracted string as a number, and check if it worked. */
69 args[i++] = strtod (token, &end);
70 if (*end && *end != *delimiter)
71 {
72 /* In case a warning is necessary
73 error(0, 0, "%s: unable to parse element %zu in '%s'\n",
74 __func__, i, convert);
75 */
76 free(copy);
77 return 0;
78 }
79 }
80 }
81 while(token && *token);
82 free (copy);
83
84 /* Check if the number of elements parsed. */
85 if (i != n)
86 {
87 /* In case a warning is necessary
88 error(0, 0, "%s: input '%s' must contain %lu numbers, but has "
89 "%lu numbers\n", __func__, convert, n, i);
90 */
91 return 0;
92 }
93
94 /* Numbers are written, return successfully. */
95 return 1;
96 }
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115 /**********************************************************************/
116 /**************** Convert string to decimal *****************/
117 /**********************************************************************/
118
119 /* Parse the right ascension input as a string in form of hh:mm:ss to a
120 * single decimal value calculated by (hh + mm / 60 + ss / 3600 ) * 15. */
121 double
gal_units_ra_to_degree(char * convert)122 gal_units_ra_to_degree(char *convert)
123 {
124 double val[3];
125 double decimal=0.0;
126
127 /* Check whether the string is successfully parsed */
128 if(gal_units_extract_decimal(convert, ":hms", val, 3))
129 {
130 /* Check whether the first value is in within limits, and add it. */
131 if(val[0]<0.0 || val[0]>24.0) return NAN;
132 decimal += val[0];
133
134 /* Check whether value of minutes is within limits, and add it. */
135 if(val[1]<0.0 || val[1]>60.0) return NAN;
136 decimal += val[1] / 60;
137
138 /* Check whether value of seconds is in within limits, and add it. */
139 if(val[2]<0.0 || val[2]>60.0) return NAN;
140 decimal += val[2] / 3600;
141
142 /* Convert value to degrees and return. */
143 decimal *= 15.0;
144 return decimal;
145 }
146 else return NAN;
147
148 /* Control shouldn't reach this point. If it does, its a bug! */
149 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
150 "problem. Control should not reach the end of this function",
151 __func__, PACKAGE_BUGREPORT);
152 return NAN;
153 }
154
155
156
157
158
159 /* Parse the declination input as a string in form of dd:mm:ss to a decimal
160 * calculated by (dd + mm / 60 + ss / 3600 ). */
161 double
gal_units_dec_to_degree(char * convert)162 gal_units_dec_to_degree(char *convert)
163 {
164 int sign;
165 double val[3], decimal=0.0;
166
167 /* Parse the values in the input string. */
168 if(gal_units_extract_decimal(convert, ":dms", val, 3))
169 {
170 /* Check whether the first value is in within limits. */
171 if(val[0]<-90.0 || val[0]>90.0) return NAN;
172
173 /* If declination is negative, the first value in the array will be
174 negative and all other values will be positive. In that case, we
175 set sign equal to -1. Therefore, we multiply the first value by
176 sign to make it positive. The final answer is again multiplied by
177 sign to make its sign same as original. */
178 sign = val[0]<0.0 ? -1 : 1;
179 decimal += val[0] * sign;
180
181 /* Check whether value of arc-minutes is in within limits. */
182 if(val[1]<0.0 || val[1]>60.0) return NAN;
183 decimal += val[1] / 60;
184
185 /* Check whether value of arc-seconds is in within limits */
186 if (val[2] < 0.0 || val[2] > 60.0) return NAN;
187 decimal += val[2] / 3600;
188
189 /* Make the sign of the decimal value same as input and return. */
190 decimal *= sign;
191 return decimal;
192 }
193 else return NAN;
194
195 /* Control shouldn't reach this point. If it does, its a bug! */
196 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
197 "problem. Control should not reach the end of this function",
198 __func__, PACKAGE_BUGREPORT);
199 return NAN;
200 }
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221 /**********************************************************************/
222 /**************** Convert decimal to string *****************/
223 /**********************************************************************/
224
225 /* Max-length of output string. */
226 #define UNITS_RADECSTR_MAXLENGTH 50
227
228 /* Parse the right ascension input as a decimal to a string in form of
229 hh:mm:ss.ss . */
230 char *
gal_units_degree_to_ra(double decimal,int usecolon)231 gal_units_degree_to_ra(double decimal, int usecolon)
232 {
233 size_t nchars;
234 int hours=0, minutes=0;
235 float seconds=0.0; /* For sub-second accuracy */
236
237 /* Allocate a long string which is large enough for string of format
238 hh:mm:ss.ss and sign */
239 char *ra=gal_pointer_allocate(GAL_TYPE_UINT8, UNITS_RADECSTR_MAXLENGTH,
240 0, __func__, "ra");
241
242 /* Check if decimal value is within bounds otherwise return error */
243 if (decimal<0 || decimal>360)
244 {
245 error (0, 0, "%s: value of decimal should be between be 0 and 360, "
246 "but is %g\n", __func__, decimal);
247 return NULL;
248 }
249
250 /* Divide decimal value by 15 and extract integer part of decimal value
251 to obtain hours */
252 decimal /= 15.0;
253 hours = (int)decimal;
254
255 /* Subtract hours from decimal and multiply remaining value by 60 to
256 obtain minutes. */
257 minutes = (int)((decimal - hours) * 60);
258
259 /* Subtract hours and minutes from decimal and multiply remaining value
260 by 3600 to obtain seconds. */
261 seconds = (decimal - hours - minutes / 60.0) * 3600;
262
263 /* Format the extracted hours, minutes and seconds as a string with
264 leading zeros if required, in hh:mm:ss format */
265 nchars = snprintf(ra, UNITS_RADECSTR_MAXLENGTH-1,
266 usecolon ? "%02d:%02d:%g" : "%02dh%02dm%gs",
267 hours, minutes, seconds);
268 if(nchars>UNITS_RADECSTR_MAXLENGTH)
269 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to address "
270 "the problem. The output string has an unreasonable length of "
271 "%zu characters", __func__, PACKAGE_BUGREPORT, nchars);
272
273 /* Return the final string. */
274 return ra;
275 }
276
277
278
279
280
281 /* Parse the declination input as a decimal to a string in form of dd:mm:ss*/
282 char *
gal_units_degree_to_dec(double decimal,int usecolon)283 gal_units_degree_to_dec(double decimal, int usecolon)
284 {
285 size_t nchars;
286 float arc_seconds=0.0;
287 int sign, degrees=0, arc_minutes=0;
288
289 /* Allocate string of fixed length which is large enough for string of
290 * format hh:mm:ss.ss and sign */
291 char *dec=gal_pointer_allocate(GAL_TYPE_UINT8, UNITS_RADECSTR_MAXLENGTH,
292 0, __func__, "ra");
293
294 /* Check if decimal value is within bounds otherwise return error */
295 if(decimal<-90 || decimal>90)
296 {
297 error (0, 0, "%s: value of decimal should be between be -90 and 90, "
298 "but is %g\n", __func__, decimal);
299 return NULL;
300 }
301
302 /* If declination is negative, we set 'sign' equal to -1. We multiply the
303 decimal by to make sure it is positive. We then extract degrees,
304 arc-minutes and arc-seconds from the decimal. Finally, we add a minus
305 sign in beginning of string if input was negative. */
306 sign = decimal<0.0 ? -1 : 1;
307 decimal *= sign;
308
309 /* Extract integer part of decimal value to obtain degrees. */
310 degrees=(int)decimal;
311
312 /* Subtract degrees from decimal and multiply remaining value by 60 to
313 obtain arc-minutes. */
314 arc_minutes=(int)( (decimal - degrees) * 60 );
315
316 /* Subtract degrees and arc-minutes from decimal and multiply remaining
317 value by 3600 to obtain arc-seconds. */
318 arc_seconds = (decimal - degrees - arc_minutes / 60.0) * 3600;
319
320 /* Format the extracted degrees, arc-minutes and arc-seconds as a string
321 with leading zeros if required, in hh:mm:ss format with correct
322 sign. */
323 nchars = snprintf(dec, UNITS_RADECSTR_MAXLENGTH-1,
324 usecolon ? "%s%02d:%02d:%g" : "%s%02dd%02dm%gs",
325 sign<0?"-":"+", degrees, arc_minutes, arc_seconds);
326 if(nchars>UNITS_RADECSTR_MAXLENGTH)
327 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to address "
328 "the problem. The output string has an unreasonable length of "
329 "%zu characters", __func__, PACKAGE_BUGREPORT, nchars);
330
331 /* Return the final string. */
332 return dec;
333 }
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354 /**********************************************************************/
355 /**************** Flux conversions *****************/
356 /**********************************************************************/
357
358 /* Convert counts to magnitude using the given zeropoint. */
359 double
gal_units_counts_to_mag(double counts,double zeropoint)360 gal_units_counts_to_mag(double counts, double zeropoint)
361 {
362 return ( counts > 0.0f
363 ? ( -2.5f * log10(counts) + zeropoint )
364 : NAN );
365 }
366
367
368
369
370
371 /* Convert magnitude to counts using the given zeropoint. */
372 double
gal_units_mag_to_counts(double mag,double zeropoint)373 gal_units_mag_to_counts(double mag, double zeropoint)
374 {
375 return pow(10, (mag - zeropoint)/(-2.5f));
376 }
377
378
379
380
381 /* Convert Pixel values to Janskys with an AB-magnitude based
382 zero-point. See the "Brightness, Flux, Magnitude and Surface
383 brightness". */
384 double
gal_units_counts_to_jy(double counts,double zeropoint_ab)385 gal_units_counts_to_jy(double counts, double zeropoint_ab)
386 {
387 return counts * 3631 * pow(10, -1 * zeropoint_ab / 2.5);
388 }
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409 /**********************************************************************/
410 /**************** Distance conversions *****************/
411 /**********************************************************************/
412 /* Convert Astronomical Units (AU) to Parsecs (PC). From the definition
413 of Parsecs, 648000/pi AU = 1 PC. So */
414 double
gal_units_au_to_pc(double au)415 gal_units_au_to_pc(double au)
416 {
417 return au / 648000.0f * 3.14159265358979323846264338327f;
418 }
419
420
421
422
423
424 /* Convert Parsecs (PC) to Astronomical units (AU), see comment of
425 'gal_units_au_to_pc'. */
426 double
gal_units_pc_to_au(double au)427 gal_units_pc_to_au(double au)
428 {
429 return au * 648000.0f / 3.14159265358979323846264338327f;
430 }
431
432
433
434
435
436 /* Convert Light-years to Parsecs, according to
437 https://en.wikipedia.org/wiki/Light-year#Definitions:
438
439 1 light-year = 9460730472580800 metres (exactly)
440 ~ 9.461 petametres
441 ~ 9.461 trillion kilometres (5.879 trillion miles)
442 ~ 63241.077 astronomical units
443 ~ 0.306601 parsecs */
444 double
gal_units_ly_to_pc(double ly)445 gal_units_ly_to_pc(double ly)
446 {
447 return ly * 0.306601f;
448 }
449
450
451
452
453
454 /* Convert Parsecs to Light-years (see comment of gal_units_ly_to_pc). */
455 double
gal_units_pc_to_ly(double pc)456 gal_units_pc_to_ly(double pc)
457 {
458 return pc / 0.306601f;
459 }
460
461
462
463
464
465 /* Convert Astronomical Units to Light-years (see comment of
466 gal_units_ly_to_pc). */
467 double
gal_units_au_to_ly(double au)468 gal_units_au_to_ly(double au)
469 {
470 return au / 63241.077f;
471 }
472
473
474
475
476
477 /* Convert Light-years to Astronomical Units (see comment of
478 gal_units_ly_to_pc). */
479 double
gal_units_ly_to_au(double ly)480 gal_units_ly_to_au(double ly)
481 {
482 return ly * 63241.077f;
483 }
484