[gnuastro-commits] master 59e2e77: Magnitude error column in MakeCatalog

[Mohammad Akhlaghi]

branch: master
commit 59e2e7728abe6a60c84cf189a6e2907afae8fe74
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Magnitude error column in MakeCatalog
    
    A new column is now added to MakeCatalog for estimating the magnitude error
    (called with the `--magnitudeerr' or `-e' option). This was done as part of
    the work for task #14124. That task is not complete with this commit,
    because we will also need to work on ways to incorporate correlated noise
    and aperture error into this error too (and the S/N calculation in
    general). The sky-region pseudo-detections and clumps can be very good
    reference to make those estimations (maybe as a correction factor).
---
 doc/gnuastro.texi         |   15 +++++++++++++-
 src/mkcatalog/args.h      |   14 ++++++++++++-
 src/mkcatalog/columns.c   |   48 ++++++++++++++++++++++++++++++++++++++-------
 src/mkcatalog/columns.h   |    2 +-
 src/mkcatalog/main.h      |    2 ++
 src/mkcatalog/mkcatalog.c |    6 +++++-
 src/mkcatalog/ui.c        |   16 +++++++++++++++
 7 files changed, 92 insertions(+), 11 deletions(-)

diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 2e16829..87743aa 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -11564,7 +11564,20 @@ with no threshold!
 
 @item -m
 @itemx --magnitude
-The magnitude of all clumps or objects, see @option{--brightness}.
+The magnitude of clumps or objects, see @option{--brightness}.
+
address@hidden -e
address@hidden --magnitudeerr
+The magnitude error of clumps or objects. The magnitude error is calculated
+from the signal-to-noise ratio (see @option{--sn}, abbreviated as
address@hidden): @mymath{\Delta{M}=2.5/(S\times\ln{10})}. Note that until now
+this error assumes un-correlated pixel values and also does not include the
+error in estimating the aperture (or error in generating the labeled
+image).
+
+For now these factors have to be found by other means.
address@hidden://savannah.gnu.org/task/index.php?14124, Task 14124} has been
+defined for work on adding these sources of error too.
 
 @item --clumpsmagnitude
 [Objects] The magnitude of all clumps in this object, see
diff --git a/src/mkcatalog/args.h b/src/mkcatalog/args.h
index 5b1c8ca..2883619 100644
--- a/src/mkcatalog/args.h
+++ b/src/mkcatalog/args.h
@@ -71,7 +71,7 @@ const char doc[] =
 
 /* Available letters for short options:
 
-   e f g k l u v w
+   f g k l u v w
    F G J L Q R U W X Y Z
 
    Number keys used: <=533
@@ -463,6 +463,14 @@ static struct argp_option options[] =
       3
     },
     {
+      "magnitudeerr",
+      'e',
+      0,
+      0,
+      "Total magnitude error.",
+      3
+    },
+    {
       "clumpsmagnitude",
       512,
       0,
@@ -790,6 +798,10 @@ parse_opt(int key, char *arg, struct argp_state *state)
       add_to_sll(&p->allcolsll, CATMAGNITUDE);
       p->up.magnitudeset=1;
       break;
+    case 'e':
+      add_to_sll(&p->allcolsll, CATMAGNITUDEERR);
+      p->up.magnitudeerrset=1;
+      break;
     case 512:
       add_to_sll(&p->allcolsll, CATCLUMPSMAGNITUDE);
       p->up.clumpsmagnitudeset=1;
diff --git a/src/mkcatalog/columns.c b/src/mkcatalog/columns.c
index 5ef2532..ed41bc5 100644
--- a/src/mkcatalog/columns.c
+++ b/src/mkcatalog/columns.c
@@ -919,10 +919,10 @@ skystd(struct mkcatalogparams *p, size_t col)
 
 
 void
-sncol(struct mkcatalogparams *p)
+sncol(struct mkcatalogparams *p, int sn0_magerr1, char *target)
 {
   size_t i;
-  double I, O, Ni, errpt, *row;
+  double sn, I, O, Ni, errpt, *row;
   size_t stdcol        = p->obj0clump1 ? CSTD        : OSTD;
   size_t areacol       = p->obj0clump1 ? CAREA       : OAREA;
   size_t brightnesscol = p->obj0clump1 ? CBrightness : OBrightness;
@@ -938,8 +938,13 @@ sncol(struct mkcatalogparams *p)
     setclumpbrightness(p);
 
   /* For the comments: */
-  p->unitp = CATUNITRATIO;
-  sprintf(p->description, "%lu: Signal to noise ratio.", p->curcol);
+  p->unitp = sn0_magerr1 ? CATUNITMAG : CATUNITRATIO;
+  if(sn0_magerr1)
+    sprintf(p->description, "%lu: %s Magnitude error.", p->curcol,
+            target);
+  else
+    sprintf(p->description, "%lu: %s signal to noise ratio.", p->curcol,
+            target);
 
   /* Calculate the signal to noise ratio. Recall that for the objects,
      the sky value was subtracted from oinfo, but for the clumps, it
@@ -995,9 +1000,38 @@ sncol(struct mkcatalogparams *p)
                     + errpt * (p->skysubtracted ? 1.0f : 2.0f) );
         }
 
-      /* Fill in the output column */
-      p->cat[i * p->numcols + p->curcol ] =
-        ( sqrt(Ni/p->cpscorr)*I / sqrt( errpt ) );
+      /* Fill in the output column. Note that magnitude error is directly
+         derivable from the S/N:
+
+         To derive the error in measuring the magnitude from the S/N, let's
+         take `F' as the flux, `Z' is the zeropoint, `M' is the magnitude,
+         `S' is the S/N, and `D' to stand for capital delta (or error in a
+         value) then from
+
+              `M = -2.5*log10(F) + Z'
+
+         we get the following equation after calculating the derivative
+         with respect to F.
+
+              `dM/df = -2.5 * ( 1 / ( F * ln(10) ) )'
+
+         From the Tailor series, `DM' can be written as:
+
+              `DM = dM/dF * DF'
+
+         So
+
+              `DM = |-2.5/ln(10)| * DF/F'
+
+         But `DF/F' is just the inverse of the Signal to noise ratio, or
+         `1/S'. So
+
+              `DM = 2.5 / ( S * ln(10) )'
+      */
+      sn = sqrt(Ni/p->cpscorr)*I / sqrt( errpt );
+      p->cat[i * p->numcols + p->curcol ] = ( sn0_magerr1
+                                              ? ( 2.5 / (sn*log(10)) )
+                                              : sn );
     }
 
 }
diff --git a/src/mkcatalog/columns.h b/src/mkcatalog/columns.h
index 886d4f5..91cce4b 100644
--- a/src/mkcatalog/columns.h
+++ b/src/mkcatalog/columns.h
@@ -76,6 +76,6 @@ void
 skystd(struct mkcatalogparams *p, size_t col);
 
 void
-sncol(struct mkcatalogparams *p);
+sncol(struct mkcatalogparams *p, int sn0_magerr1, char *target);
 
 #endif
diff --git a/src/mkcatalog/main.h b/src/mkcatalog/main.h
index 3160b1f..483e207 100644
--- a/src/mkcatalog/main.h
+++ b/src/mkcatalog/main.h
@@ -190,6 +190,7 @@ enum outcols
     CATCLUMPSBRIGHTNESS,
     CATNORIVERBRIGHTNESS,
     CATMAGNITUDE,
+    CATMAGNITUDEERR,
     CATCLUMPSMAGNITUDE,
     CATRIVERAVE,
     CATRIVERNUM,
@@ -269,6 +270,7 @@ struct uiparams
   int     clumpsbrightnessset;
   int    noriverbrightnessset;
   int            magnitudeset;
+  int         magnitudeerrset;
   int      clumpsmagnitudeset;
   int             riveraveset;
   int             rivernumset;
diff --git a/src/mkcatalog/mkcatalog.c b/src/mkcatalog/mkcatalog.c
index de86d9f..c9897d3 100644
--- a/src/mkcatalog/mkcatalog.c
+++ b/src/mkcatalog/mkcatalog.c
@@ -652,6 +652,10 @@ makeoutput(struct mkcatalogparams *p)
               brightnessmag(p, tmpcol, target, MKCATMAG);
               break;
 
+            case CATMAGNITUDEERR:
+              sncol(p, 1, target);
+              break;
+
             case CATCLUMPSMAGNITUDE:
               brightnessmag(p, OBrightnessC, MKCATCINO, MKCATMAG);
               break;
@@ -665,7 +669,7 @@ makeoutput(struct mkcatalogparams *p)
               break;
 
             case CATSN:
-              sncol(p);
+              sncol(p, 0, target);
               break;
 
             case CATSKY:
diff --git a/src/mkcatalog/ui.c b/src/mkcatalog/ui.c
index 19daaa2..f6df921 100644
--- a/src/mkcatalog/ui.c
+++ b/src/mkcatalog/ui.c
@@ -440,6 +440,15 @@ readconfig(char *filename, struct mkcatalogparams *p)
           add_to_sll(&p->allcolsll, CATMAGNITUDE);
           up->magnitudeset=1;
         }
+      else if(strcmp(name, "magnitudeerr")==0)
+        {
+          if(up->magnitudeerrset) continue;
+          gal_checkset_int_zero_or_one(value, &yes, name, key, SPACK,
+                                       filename, lineno);
+          if(!yes) continue;
+          add_to_sll(&p->allcolsll, CATMAGNITUDEERR);
+          up->magnitudeerrset=1;
+        }
       else if(strcmp(name, "clumpsmagnitude")==0)
         {
           if(up->clumpsmagnitudeset) continue;
@@ -716,6 +725,9 @@ printvalues(FILE *fp, struct mkcatalogparams *p)
       case CATMAGNITUDE:
         fprintf(fp, CONF_SHOWFMT"%d\n", "magnitude", 1);
         break;
+      case CATMAGNITUDEERR:
+        fprintf(fp, CONF_SHOWFMT"%d\n", "magnitudeerr", 1);
+        break;
       case CATCLUMPSMAGNITUDE:
         fprintf(fp, CONF_SHOWFMT"%d\n", "clumpsmagnitude", 1);
         break;
@@ -1193,6 +1205,10 @@ preparearrays(struct mkcatalogparams *p)
           p->objcols[p->objncols++] = p->allcols[i];
           p->clumpcols[p->clumpncols++] = p->allcols[i];
           break;
+        case CATMAGNITUDEERR:
+          p->objcols[p->objncols++] = p->allcols[i];
+          p->clumpcols[p->clumpncols++] = p->allcols[i];
+          break;
         case CATCLUMPSMAGNITUDE:
           p->objcols[p->objncols++] = p->allcols[i];
           break;