[gnuastro-commits] master 8b13edb: Library (units.h and arithmetic.h): new counts-to-mag function

[Mohammad Akhlaghi]

branch: master
commit 8b13edbf7a8f5de0b63aabe2554d84fb71951cee
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Library (units.h and arithmetic.h): new counts-to-mag function
    
    Until now, when we wanted to convert counts to magnitudes anywhere in
    Gnuastro we would use the raw equation. This could be buggy (mistakenly
    missing a negative, or etc!). On the other hand, for a Gnuastro user,
    converting counts to magnitudes for a certain count level required a
    relatively long and inconvenient AWK or Arithmetic command!
    
    With this commit, the 'gal_units_counts_to_mag' function has been added
    inside of Gnuastro's library. Thanks to that function, a new operator has
    been added to the Arithmetic library (and thus the Arithmetic program or
    Column Arithmetic in Table). Also, with the function in place, all the
    places where we had the raw magnitude estimations (using '-2.5*log10(B)+Z')
    are replaced with a call to this function and thus the code has been
    cleaned up.
---
 NEWS                      |  2 ++
 bin/mkcatalog/columns.c   |  3 ++-
 bin/mkcatalog/mkcatalog.c | 16 ++++++++--------
 bin/mknoise/mknoise.c     |  3 ++-
 bin/mkprof/mkprof.c       |  3 ++-
 doc/gnuastro.texi         | 21 ++++++++++++++++++++-
 lib/arithmetic.c          |  6 ++++++
 lib/gnuastro/arithmetic.h |  1 +
 lib/gnuastro/units.h      | 11 +++++++----
 lib/units.c               | 13 +++++++++++++
 10 files changed, 63 insertions(+), 16 deletions(-)

diff --git a/NEWS b/NEWS
index 6d05a81..d17b68e 100644
--- a/NEWS
+++ b/NEWS
@@ -36,6 +36,7 @@ See the end of the file for license conditions.
      - asinh: Inverse of hyperbolic sine.
      - acosh: Inverse of hyperbolic cosine.
      - atabh: Inverse of hyperbolic tangent.
+     - counts-to-mag: Convert counts to magnitudes with given zero point.
      - counts-to-jy: Convert counts to Janskys through a zero point based
           on AB magnitudes.
 
@@ -79,6 +80,7 @@ See the end of the file for license conditions.
          astquery ned --dataset=extinction --center=49.9507,41.5116
 
   Library:
+   - gal_units_counts_to_mag: Convert counts to magnitudes.
    - gal_units_counts_to_jy: Convert counts to Janskys.
    - New arithmetic operator macros (for the 'gal_arithmetic' function):
      - GAL_ARITHMETIC_OP_SIN: sine (input in deg).
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index e94bb79..69cc29d 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <pthread.h>
 
 #include <gnuastro/wcs.h>
+#include <gnuastro/units.h>
 #include <gnuastro/pointer.h>
 
 #include <gnuastro-internal/checkset.h>
@@ -1995,7 +1996,7 @@ columns_define_alloc(struct mkcatalogparams *p)
 /**********            Column calculation           ***************/
 /******************************************************************/
 #define MKC_RATIO(TOP,BOT) ( (BOT)!=0.0f ? (TOP)/(BOT) : NAN )
-#define MKC_MAG(B) ( ((B)>0) ? -2.5f * log10(B) + p->zeropoint : NAN )
+#define MKC_MAG(B) ( gal_units_counts_to_mag(B, p->zeropoint) )
 #define MKC_SB(B, A) ( ((B)>0 && (A)>0)                                 \
                        ? MKC_MAG(B) + 2.5f * log10((A) * p->pixelarcsecsq) \
                        : NAN )
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index e3e75b7..55fe6e7 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -35,6 +35,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/wcs.h>
 #include <gnuastro/data.h>
 #include <gnuastro/fits.h>
+#include <gnuastro/units.h>
 #include <gnuastro/threads.h>
 #include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
@@ -495,10 +496,9 @@ mkcatalog_outputs_same_start(struct mkcatalogparams *p, 
int o0c1,
         {
           /* Per pixel. */
           if( asprintf(&str, "%g sigma surface brightness (magnitude/pixel): "
-                       "%.3f", p->sfmagnsigma, ( -2.5f
-                                                 *log10( p->sfmagnsigma
-                                                         * p->medstd )
-                                                 + p->zeropoint ) )<0 )
+                       "%.3f", p->sfmagnsigma,
+                       gal_units_counts_to_mag(p->sfmagnsigma * p->medstd,
+                                               p->zeropoint) )<0 )
             error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
           gal_list_str_add(&comments, str, 0);
 
@@ -520,10 +520,10 @@ mkcatalog_outputs_same_start(struct mkcatalogparams *p, 
int o0c1,
               if( asprintf(&str, "%g sigma surface brightness "
                            "(magnitude/%sarcsec^2): %.3f", p->sfmagnsigma,
                            tstr ? tstr : "",
-                           ( -2.5f * log10( p->sfmagnsigma
-                                            * p->medstd
-                                            * sqrt( p->sfmagarea / pixarea) )
-                             + p->zeropoint ) )<0 )
+                           gal_units_counts_to_mag(p->sfmagnsigma
+                                                   * p->medstd
+                                                   * sqrt( p->sfmagarea / 
pixarea),
+                                                   p->zeropoint) ) <0 )
                 error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
 
               /* Add the final string/line to the catalog comments. */
diff --git a/bin/mknoise/mknoise.c b/bin/mknoise/mknoise.c
index 6731ed7..5320d5c 100644
--- a/bin/mknoise/mknoise.c
+++ b/bin/mknoise/mknoise.c
@@ -32,6 +32,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gsl/gsl_rng.h>         /* Used in setrandoms.   */
 #include <gnuastro/fits.h>
+#include <gnuastro/units.h>
 #include <gsl/gsl_randist.h>     /* To make noise.        */
 
 #include <gnuastro-internal/timing.h>
@@ -70,7 +71,7 @@ convertsaveoutput(struct mknoiseparams *p)
                                 0, NULL, 0);
       if( !isnan(p->zeropoint) )
         {
-          tmp=-2.5 * log10(p->background) + p->zeropoint;
+          tmp=gal_units_counts_to_mag(p->background, p->zeropoint);
           gal_checkset_allocate_copy("BCKGMAG", &keyname);
           gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname, 1,
                                     &tmp, 0,
diff --git a/bin/mkprof/mkprof.c b/bin/mkprof/mkprof.c
index de7144b..13b8cde 100644
--- a/bin/mkprof/mkprof.c
+++ b/bin/mkprof/mkprof.c
@@ -32,6 +32,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/box.h>
 #include <gnuastro/git.h>
 #include <gnuastro/fits.h>
+#include <gnuastro/units.h>
 #include <gnuastro/threads.h>
 #include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
@@ -655,7 +656,7 @@ mkprof_write(struct mkprofparams *p)
                 break;
               case 2:
                 ((float *)(log->array))[ibq->id] =
-                  sum>0.0f ? -2.5f*log10(sum)+p->zeropoint : NAN;
+                  gal_units_counts_to_mag(sum, p->zeropoint);
                 break;
               case 1:
                 ((unsigned long *)(log->array))[ibq->id]=ibq->id+1;
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 4122f32..3aba1d1 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -11688,10 +11688,23 @@ These operators take a single operand.
 Inverse Hyperbolic sine, cosine, and tangent.
 These operators take a single operand.
 
+@item counts-to-mag
+Convert counts (usually CCD outputs) to magnitudes using the given zeropoint.
+The zero point is the first popped operand and the count value is the second.
+For example assume you have measured the standard deviation of the noise in an 
image to be @code{0.1}, and the image's zero point is @code{22.5}.
+You can therefore measure the @emph{per-pixel} surface brightness limit of the 
dataset (which is the magnitude of the noise standrard deviation) with this 
simple command below.
+Note that because the output is a simple number, we are using @option{--quiet} 
to avoid printing extra information.
+
+@example
+astarithmetic 0.1 22.5 counts-to-mag --quiet
+@end example
+
+Of course, you can also convert every pixel in an image (or table column in 
Table's @ref{Column arithmetic}) with this operator if you replace the second 
popped operand with an image/column.
+
 @item counts-to-jy
 @cindex AB magnitude
 @cindex Magnitude, AB
-Convert counts (CCD outputs) to Janskys through an AB-magnitude based 
zeropoint.
+Convert counts (usually CCD outputs) to Janskys through an AB-magnitude based 
zeropoint.
 The top-popped operand is assumed to be the AB-magnitude zero point and the 
second-popped operand is assumed to be a dataset in units of counts (an image 
in Arithmetic, and a column in Table's @ref{Column arithmetic}).
 For the full equation and basic definitions, see @ref{Brightness flux 
magnitude}.
 
@@ -28387,6 +28400,12 @@ Convert the input Declination (Dec) degree (a single 
floating point number) to o
 If @code{usecolon!=0}, then the delimiters between the components will be 
colons: @code{_:_:_}.
 @end deftypefun
 
+@deftypefun double gal_units_counts_to_mag (double @code{counts}, double 
@code{zeropoint})
+@cindex Magnitude
+Convert counts to magnitudes through the given zero point.
+For more on the equation, see @ref{Brightness flux magnitude}.
+@end deftypefun
+
 @deftypefun double gal_units_counts_to_jy (double @code{counts}, double 
@code{zeropoint_ab})
 @cindex Jansky (Jy)
 @cindex AB Magnitude
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 5f427cc..8f98575 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -1777,6 +1777,8 @@ arithmetic_function_binary_flt(int operator, int flags, 
gal_data_t *il,
       BINFUNC_F_OPERATOR_SET( pow,   +0 );         break;
     case GAL_ARITHMETIC_OP_ATAN2:
       BINFUNC_F_OPERATOR_SET( atan2, *180.0f/pi ); break;
+    case GAL_ARITHMETIC_OP_COUNTS_TO_MAG:
+      BINFUNC_F_OPERATOR_SET( gal_units_counts_to_mag, +0 ); break;
     case GAL_ARITHMETIC_OP_COUNTS_TO_JY:
       BINFUNC_F_OPERATOR_SET( gal_units_counts_to_jy, +0 ); break;
     default:
@@ -1962,6 +1964,8 @@ gal_arithmetic_set_operator(char *string, size_t 
*num_operands)
     { op=GAL_ARITHMETIC_OP_DEGREE_TO_RA;      *num_operands=1;  }
   else if (!strcmp(string, "degree-to-dec"))
     { op=GAL_ARITHMETIC_OP_DEGREE_TO_DEC;     *num_operands=1;  }
+  else if (!strcmp(string, "counts-to-mag"))
+    { op=GAL_ARITHMETIC_OP_COUNTS_TO_MAG;     *num_operands=2;  }
   else if (!strcmp(string, "counts-to-jy"))
     { op=GAL_ARITHMETIC_OP_COUNTS_TO_JY;      *num_operands=2;  }
 
@@ -2137,6 +2141,7 @@ gal_arithmetic_operator_string(int operator)
     case GAL_ARITHMETIC_OP_DEC_TO_DEGREE:   return "dec-to-degree";
     case GAL_ARITHMETIC_OP_DEGREE_TO_RA:    return "degree-to-ra";
     case GAL_ARITHMETIC_OP_DEGREE_TO_DEC:   return "degree-to-dec";
+    case GAL_ARITHMETIC_OP_COUNTS_TO_MAG:   return "counts-to-mag";
     case GAL_ARITHMETIC_OP_COUNTS_TO_JY:    return "counts-to-jy";
 
     case GAL_ARITHMETIC_OP_MINVAL:          return "minvalue";
@@ -2260,6 +2265,7 @@ gal_arithmetic(int operator, size_t numthreads, int 
flags, ...)
     /* Binary function operators. */
     case GAL_ARITHMETIC_OP_POW:
     case GAL_ARITHMETIC_OP_ATAN2:
+    case GAL_ARITHMETIC_OP_COUNTS_TO_MAG:
     case GAL_ARITHMETIC_OP_COUNTS_TO_JY:
       d1 = va_arg(va, gal_data_t *);
       d2 = va_arg(va, gal_data_t *);
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index 6ad36b0..49d6930 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -124,6 +124,7 @@ enum gal_arithmetic_operators
   GAL_ARITHMETIC_OP_DEC_TO_DEGREE,/* declination to decimal. */
   GAL_ARITHMETIC_OP_DEGREE_TO_RA, /* right ascension to decimal. */
   GAL_ARITHMETIC_OP_DEGREE_TO_DEC,/* declination to decimal. */
+  GAL_ARITHMETIC_OP_COUNTS_TO_MAG,/* Counts to magnitude. */
   GAL_ARITHMETIC_OP_COUNTS_TO_JY, /* Counts to Janskys with AB-mag zeropoint. 
*/
 
   GAL_ARITHMETIC_OP_MINVAL,       /* Minimum value of array.               */
diff --git a/lib/gnuastro/units.h b/lib/gnuastro/units.h
index 52ee3b8..3243f1f 100644
--- a/lib/gnuastro/units.h
+++ b/lib/gnuastro/units.h
@@ -62,16 +62,19 @@ gal_units_extract_decimal(char *convert, const char 
*delimiter,
                           double *args, size_t n);
 
 double
-gal_units_ra_to_degree (char *convert);
+gal_units_ra_to_degree(char *convert);
 
 double
-gal_units_dec_to_degree (char *convert);
+gal_units_dec_to_degree(char *convert);
 
 char *
-gal_units_degree_to_ra (double decimal, int usecolon);
+gal_units_degree_to_ra(double decimal, int usecolon);
 
 char *
-gal_units_degree_to_dec (double decimal, int usecolon);
+gal_units_degree_to_dec(double decimal, int usecolon);
+
+double
+gal_units_counts_to_mag(double counts, double zeropoint);
 
 double
 gal_units_counts_to_jy(double counts, double zeropoint_ab);
diff --git a/lib/units.c b/lib/units.c
index b20cf65..ce9e873 100644
--- a/lib/units.c
+++ b/lib/units.c
@@ -355,6 +355,19 @@ gal_units_degree_to_dec(double decimal, int usecolon)
 /****************          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 Pixel values to Janskys with an AB-magnitude based
    zero-point. See the "Brightness, Flux, Magnitude and Surface
    brightness". */