/********************************************************************* Units -- Convert data from one unit to other. This is part of GNU Astronomy Utilities (Gnuastro) package. Original author: Kartik Ohri Contributing author(s): Mohammad Akhlaghi Copyright (C) 2020-2021, Free Software Foundation, Inc. Gnuastro is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. Gnuastro is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with Gnuastro. If not, see . **********************************************************************/ #include #include #include #include #include #include #include #include #include /**********************************************************************/ /**************** Functions to parse strings *****************/ /**********************************************************************/ /* Parse the input string consisting of numbers separated by given delimiter into an array. */ int gal_units_extract_decimal(char *convert, const char *delimiter, double *args, size_t n) { size_t i = 0; char *copy, *token, *end; /* Create a copy of the string to be parsed and parse it. This is because it will be modified during the parsing. */ copy=strdup(convert); do { /* Check if the required number of arguments are passed */ if(i==n+1) { free(copy); error(0, 0, "%s: input '%s' exceeds maximum number of arguments " "(%zu)", __func__, convert, n); return 0; } /* Extract the substring till the next delimiter */ token=strtok(i==0?copy:NULL, delimiter); if(token) { /* Parse extracted string as a number, and check if it worked. */ args[i++] = strtod (token, &end); if (*end && *end != *delimiter) { /* In case a warning is necessary error(0, 0, "%s: unable to parse element %zu in '%s'\n", __func__, i, convert); */ free(copy); return 0; } } } while(token && *token); free (copy); /* Check if the number of elements parsed. */ if (i != n) { /* In case a warning is necessary error(0, 0, "%s: input '%s' must contain %lu numbers, but has " "%lu numbers\n", __func__, convert, n, i); */ return 0; } /* Numbers are written, return successfully. */ return 1; } /**********************************************************************/ /**************** Convert string to decimal *****************/ /**********************************************************************/ /* Parse the right ascension input as a string in form of hh:mm:ss to a * single decimal value calculated by (hh + mm / 60 + ss / 3600 ) * 15. */ double gal_units_ra_to_degree(char *convert) { double val[3]; double decimal=0.0; /* Check whether the string is successfully parsed */ if(gal_units_extract_decimal(convert, ":hms", val, 3)) { /* Check whether the first value is in within limits, and add it. */ if(val[0]<0.0 || val[0]>24.0) return NAN; decimal += val[0]; /* Check whether value of minutes is within limits, and add it. */ if(val[1]<0.0 || val[1]>60.0) return NAN; decimal += val[1] / 60; /* Check whether value of seconds is in within limits, and add it. */ if(val[2]<0.0 || val[2]>60.0) return NAN; decimal += val[2] / 3600; /* Convert value to degrees and return. */ decimal *= 15.0; return decimal; } else return NAN; /* Control shouldn't reach this point. If it does, its a bug! */ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the " "problem. Control should not reach the end of this function", __func__, PACKAGE_BUGREPORT); return NAN; } /* Parse the declination input as a string in form of dd:mm:ss to a decimal * calculated by (dd + mm / 60 + ss / 3600 ). */ double gal_units_dec_to_degree(char *convert) { int sign; double val[3], decimal=0.0; /* Parse the values in the input string. */ if(gal_units_extract_decimal(convert, ":dms", val, 3)) { /* Check whether the first value is in within limits. */ if(val[0]<-90.0 || val[0]>90.0) return NAN; /* If declination is negative, the first value in the array will be negative and all other values will be positive. In that case, we set sign equal to -1. Therefore, we multiply the first value by sign to make it positive. The final answer is again multiplied by sign to make its sign same as original. */ sign = val[0]<0.0 ? -1 : 1; decimal += val[0] * sign; /* Check whether value of arc-minutes is in within limits. */ if(val[1]<0.0 || val[1]>60.0) return NAN; decimal += val[1] / 60; /* Check whether value of arc-seconds is in within limits */ if (val[2] < 0.0 || val[2] > 60.0) return NAN; decimal += val[2] / 3600; /* Make the sign of the decimal value same as input and return. */ decimal *= sign; return decimal; } else return NAN; /* Control shouldn't reach this point. If it does, its a bug! */ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the " "problem. Control should not reach the end of this function", __func__, PACKAGE_BUGREPORT); return NAN; } /**********************************************************************/ /**************** Convert decimal to string *****************/ /**********************************************************************/ /* Max-length of output string. */ #define UNITS_RADECSTR_MAXLENGTH 50 /* Parse the right ascension input as a decimal to a string in form of hh:mm:ss.ss . */ char * gal_units_degree_to_ra(double decimal, int usecolon) { size_t nchars; int hours=0, minutes=0; float seconds=0.0; /* For sub-second accuracy */ /* Allocate a long string which is large enough for string of format hh:mm:ss.ss and sign */ char *ra=gal_pointer_allocate(GAL_TYPE_UINT8, UNITS_RADECSTR_MAXLENGTH, 0, __func__, "ra"); /* Check if decimal value is within bounds otherwise return error */ if (decimal<0 || decimal>360) { error (0, 0, "%s: value of decimal should be between be 0 and 360, " "but is %g\n", __func__, decimal); return NULL; } /* Divide decimal value by 15 and extract integer part of decimal value to obtain hours */ decimal /= 15.0; hours = (int)decimal; /* Subtract hours from decimal and multiply remaining value by 60 to obtain minutes. */ minutes = (int)((decimal - hours) * 60); /* Subtract hours and minutes from decimal and multiply remaining value by 3600 to obtain seconds. */ seconds = (decimal - hours - minutes / 60.0) * 3600; /* Format the extracted hours, minutes and seconds as a string with leading zeros if required, in hh:mm:ss format */ nchars = snprintf(ra, UNITS_RADECSTR_MAXLENGTH-1, usecolon ? "%02d:%02d:%g" : "%02dh%02dm%gs", hours, minutes, seconds); if(nchars>UNITS_RADECSTR_MAXLENGTH) error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to address " "the problem. The output string has an unreasonable length of " "%zu characters", __func__, PACKAGE_BUGREPORT, nchars); /* Return the final string. */ return ra; } /* Parse the declination input as a decimal to a string in form of dd:mm:ss*/ char * gal_units_degree_to_dec(double decimal, int usecolon) { size_t nchars; float arc_seconds=0.0; int sign, degrees=0, arc_minutes=0; /* Allocate string of fixed length which is large enough for string of * format hh:mm:ss.ss and sign */ char *dec=gal_pointer_allocate(GAL_TYPE_UINT8, UNITS_RADECSTR_MAXLENGTH, 0, __func__, "ra"); /* Check if decimal value is within bounds otherwise return error */ if(decimal<-90 || decimal>90) { error (0, 0, "%s: value of decimal should be between be -90 and 90, " "but is %g\n", __func__, decimal); return NULL; } /* If declination is negative, we set 'sign' equal to -1. We multiply the decimal by to make sure it is positive. We then extract degrees, arc-minutes and arc-seconds from the decimal. Finally, we add a minus sign in beginning of string if input was negative. */ sign = decimal<0.0 ? -1 : 1; decimal *= sign; /* Extract integer part of decimal value to obtain degrees. */ degrees=(int)decimal; /* Subtract degrees from decimal and multiply remaining value by 60 to obtain arc-minutes. */ arc_minutes=(int)( (decimal - degrees) * 60 ); /* Subtract degrees and arc-minutes from decimal and multiply remaining value by 3600 to obtain arc-seconds. */ arc_seconds = (decimal - degrees - arc_minutes / 60.0) * 3600; /* Format the extracted degrees, arc-minutes and arc-seconds as a string with leading zeros if required, in hh:mm:ss format with correct sign. */ nchars = snprintf(dec, UNITS_RADECSTR_MAXLENGTH-1, usecolon ? "%s%02d:%02d:%g" : "%s%02dd%02dm%gs", sign<0?"-":"+", degrees, arc_minutes, arc_seconds); if(nchars>UNITS_RADECSTR_MAXLENGTH) error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to address " "the problem. The output string has an unreasonable length of " "%zu characters", __func__, PACKAGE_BUGREPORT, nchars); /* Return the final string. */ return dec; } /**********************************************************************/ /**************** Flux conversions *****************/ /**********************************************************************/ /* Convert counts to magnitude using the given zeropoint. */ double gal_units_counts_to_mag(double counts, double zeropoint) { return ( counts > 0.0f ? ( -2.5f * log10(counts) + zeropoint ) : NAN ); } /* Convert magnitude to counts using the given zeropoint. */ double gal_units_mag_to_counts(double mag, double zeropoint) { return pow(10, (mag - zeropoint)/(-2.5f)); } /* Convert Pixel values to Janskys with an AB-magnitude based zero-point. See the "Brightness, Flux, Magnitude and Surface brightness". */ double gal_units_counts_to_jy(double counts, double zeropoint_ab) { return counts * 3631 * pow(10, -1 * zeropoint_ab / 2.5); } /**********************************************************************/ /**************** Distance conversions *****************/ /**********************************************************************/ /* Convert Astronomical Units (AU) to Parsecs (PC). From the definition of Parsecs, 648000/pi AU = 1 PC. So */ double gal_units_au_to_pc(double au) { return au / 648000.0f * 3.14159265358979323846264338327f; } /* Convert Parsecs (PC) to Astronomical units (AU), see comment of 'gal_units_au_to_pc'. */ double gal_units_pc_to_au(double au) { return au * 648000.0f / 3.14159265358979323846264338327f; } /* Convert Light-years to Parsecs, according to https://en.wikipedia.org/wiki/Light-year#Definitions: 1 light-year = 9460730472580800 metres (exactly) ~ 9.461 petametres ~ 9.461 trillion kilometres (5.879 trillion miles) ~ 63241.077 astronomical units ~ 0.306601 parsecs */ double gal_units_ly_to_pc(double ly) { return ly * 0.306601f; } /* Convert Parsecs to Light-years (see comment of gal_units_ly_to_pc). */ double gal_units_pc_to_ly(double pc) { return pc / 0.306601f; } /* Convert Astronomical Units to Light-years (see comment of gal_units_ly_to_pc). */ double gal_units_au_to_ly(double au) { return au / 63241.077f; } /* Convert Light-years to Astronomical Units (see comment of gal_units_ly_to_pc). */ double gal_units_ly_to_au(double ly) { return ly * 63241.077f; }