[gnuastro-commits] master 4979e34: Table: angular-distance is a new column arithmetic operator

[Mohammad Akhlaghi]

branch: master
commit 4979e34ab427a00b4f3bcd2f2538b12871ddc217
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Table: angular-distance is a new column arithmetic operator
    
    `angular-distance' is a new operator to easily find the distance between
    points in various table columns, or the distances of all the rows with a
    fixed point.
---
 NEWS                   |  5 ++++
 bin/table/arithmetic.c | 76 +++++++++++++++++++++++++++++++++++++++++++++++---
 bin/table/arithmetic.h |  1 +
 doc/gnuastro.texi      | 29 +++++++++++++++++++
 4 files changed, 107 insertions(+), 4 deletions(-)

diff --git a/NEWS b/NEWS
index 14c5897..d777388 100644
--- a/NEWS
+++ b/NEWS
@@ -22,6 +22,11 @@ See the end of the file for license conditions.
      that have a value of 2, 4 and 5 in the `ID' column.
    --notequal: Output only rows that have a different value compared to the
      values given to this option in the given column.
+   - Column Arithmetic operators:
+     - `angular-distance': a new operator to easily find the angular
+       distance (along a great circle) between points in various table
+       columns, or the distances of all the points in the table rows with a
+       fixed point. See the book for examples and better explanation.
 
 ** Removed features
 
diff --git a/bin/table/arithmetic.c b/bin/table/arithmetic.c
index 5eb073f..279a9a5 100644
--- a/bin/table/arithmetic.c
+++ b/bin/table/arithmetic.c
@@ -106,6 +106,7 @@ arithmetic_operator_name(int operator)
       {
       case ARITHMETIC_TABLE_OP_WCSTOIMG: out="wcstoimg"; break;
       case ARITHMETIC_TABLE_OP_IMGTOWCS: out="imgtowcs"; break;
+      case ARITHMETIC_TABLE_OP_ANGULARDISTANCE: out="angular-distance"; break;
       default:
         error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
               "the problem. %d is not a recognized operator code", __func__,
@@ -153,9 +154,11 @@ arithmetic_set_operator(struct tableparams *p, char 
*string,
   if( op==GAL_ARITHMETIC_OP_INVALID )
     {
       if(      !strncmp(string, "wcstoimg", 8))
-        { op=ARITHMETIC_TABLE_OP_WCSTOIMG;   *num_operands=0; }
+        { op=ARITHMETIC_TABLE_OP_WCSTOIMG; *num_operands=0; }
       else if (!strncmp(string, "imgtowcs", 8))
-        { op=ARITHMETIC_TABLE_OP_IMGTOWCS;   *num_operands=0; }
+        { op=ARITHMETIC_TABLE_OP_IMGTOWCS; *num_operands=0; }
+      else if (!strncmp(string, "angular-distance", 8))
+        { op=ARITHMETIC_TABLE_OP_ANGULARDISTANCE; *num_operands=0; }
       else
         { op=GAL_ARITHMETIC_OP_INVALID; *num_operands=GAL_BLANK_INT; }
     }
@@ -430,6 +433,68 @@ arithmetic_wcs(struct tableparams *p, gal_data_t **stack, 
int operator)
 
 
 
+static void
+arithmetic_angular_dist(struct tableparams *p, gal_data_t **stack, int 
operator)
+{
+  size_t i, j;
+  double *o, *a1, *a2, *b1, *b2;
+  gal_data_t *a, *b, *tmp, *out;
+
+  /* Pop the columns for point `b'.*/
+  tmp=arithmetic_stack_pop(stack, operator);
+  tmp=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT64);
+  b=arithmetic_stack_pop(stack, operator);
+  b=gal_data_copy_to_new_type_free(b, GAL_TYPE_FLOAT64);
+  b->next=tmp;
+
+  /* Pop the columns for point `a'.*/
+  tmp=arithmetic_stack_pop(stack, operator);
+  tmp=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT64);
+  a=arithmetic_stack_pop(stack, operator);
+  a=gal_data_copy_to_new_type_free(a, GAL_TYPE_FLOAT64);
+  a->next=tmp;
+
+  /* Make sure the sizes are consistant: note that one point can have a
+     single coordinate, but we don't know which one. */
+  if(a->size!=a->next->size)
+    error(EXIT_FAILURE, 0, "the sizes of the third and fourth operands "
+          "of the `%s' operator (respectively containing %zu and %zu "
+          "numbers) must be equal", arithmetic_operator_name(operator),
+          a->next->size, a->size);
+  if(b->size!=b->next->size)
+    error(EXIT_FAILURE, 0, "the sizes of the third and fourth operands "
+          "of the `%s' operator (respectively containing %zu and %zu "
+          "numbers) must be equal", arithmetic_operator_name(operator),
+          b->next->size, b->size);
+
+  /* Make the output array based on the largest size. */
+  out=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1,
+                     (a->size>b->size ? &a->size : &b->size), NULL, 0,
+                     p->cp.minmapsize, p->cp.quietmmap, "angular-dist",
+                     NULL, "Angular distance between input points");
+
+  /* Measure the distances.  */
+  o=out->array;
+  a1=a->array; a2=a->next->array;
+  b1=b->array; b2=b->next->array;
+  if(a->size==1 || b->size==1) /* One of them is a single point. */
+    for(i=0;i<a->size;++i)
+      for(j=0;j<b->size;++j)
+        o[a->size>b->size?i:j] = gal_wcs_angular_distance_deg(a1[i], a2[i],
+                                                              b1[j], b2[j]);
+  else                         /* Both have the same length. */
+    for(i=0;i<a->size;++i)     /* (all were originally from the same table) */
+      o[i] = gal_wcs_angular_distance_deg(a1[i], a2[i], b1[i], b2[i]);
+
+  /* Clean up and put the output dataset onto the stack. */
+  gal_list_data_free(a);
+  gal_list_data_free(b);
+  gal_list_data_add(stack, out);
+}
+
+
+
+
 
 
 
@@ -545,6 +610,10 @@ arithmetic_operator_run(struct tableparams *p, gal_data_t 
**stack,
           arithmetic_wcs(p, stack, operator);
           break;
 
+        case ARITHMETIC_TABLE_OP_ANGULARDISTANCE:
+          arithmetic_angular_dist(p, stack, operator);
+          break;
+
         default:
           error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
                 "fix the problem. The operator code %d is not recognized",
@@ -605,7 +674,7 @@ arithmetic_reverse_polish(struct tableparams *p, struct 
column_pack *outpack)
       /* A small sanity check. */
       if(single->size==1 && p->table && single->size!=p->table->size)
         error(EXIT_FAILURE, 0, "the arithmetic operation resulted in a "
-              "a single value, but other columns have also been requested "
+              "single value, but other columns have also been requested "
               "which have more elements/rows");
 
       /* Set `single->next' to NULL so it isn't treated as a list and
@@ -643,7 +712,6 @@ arithmetic_operate(struct tableparams *p)
   size_t i;
   struct column_pack *outpack;
 
-
   /* From now on, we will be looking for columns from the index in
      `colarray', so to keep things clean, we'll set all the `next' elements
      to NULL. */
diff --git a/bin/table/arithmetic.h b/bin/table/arithmetic.h
index dc75d61..7695484 100644
--- a/bin/table/arithmetic.h
+++ b/bin/table/arithmetic.h
@@ -37,6 +37,7 @@ enum arithmetic_operators
 {
  ARITHMETIC_TABLE_OP_WCSTOIMG = GAL_ARITHMETIC_OP_LAST_CODE,
  ARITHMETIC_TABLE_OP_IMGTOWCS,
+ ARITHMETIC_TABLE_OP_ANGULARDISTANCE,
 };
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 48e4e7a..af7d08a 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -8801,6 +8801,35 @@ $ asttable table.fits -cID,RA -cDEC \
 
 @item imgtowcs
 Similar to @code{wcstoimg}, except that image/dataset coordinates are 
converted to WCS coordinates.
+
+@item angular-distance
+Return the spherical angular distance (along a great circle, in degrees) 
between the given two points.
+Note that each point needs two coordinates (in degrees), so this operator 
needs four operands.
+The first and second popped operands are considered to belong to one point and 
the third and fourth popped operands to the second point.
+
+Each of the input points can be a single coordinate or a full table column 
(containing many points.
+In other words, the following commands are all valid:
+
+@example
+$ asttable table.fits \
+           -c"arith RA1 DEC1 RA2 DEC2 angular-distance"
+$ asttable table.fits \
+           -c"arith RA DEC 12.345 6.789 angular-distance"
+$ asttable table.fits \
+           -c"arith 12.345 6.789 RA DEC angular-distance"
+@end example
+
+In the first case we are assuming that @file{table.fits} has the following 
four columns @code{RA1}, @code{DEC1}, @code{RA2}, @code{DEC2}.
+The returned column by this operator would be the difference between two 
points in each row with coordinates like the following (@code{RA1}, 
@code{DEC1}) and (@code{RA2}, @code{DEC2}).
+In other words, for each row, the angular distance between different points is 
calculated.
+
+In the second and third cases (which are identical), it is assumed that 
@file{table.fits} has the two columns @code{RA} and @code{DEC}.
+The returned column by this operator will be the difference of each row with 
the fixed point of (12.345, 6.789).
+
+The distance (along a great circle) on a sphere between two points is 
calculated with the equation below, where @mymath{r_1}, @mymath{r_2}, 
@mymath{d_1} and @mymath{d_2} are the right ascensions and declinations of 
points 1 and 2.
+
+@dispmath {\cos(d)=\sin(d_1)\sin(d_2)+\cos(d_1)\cos(d_2)\cos(r_1-r_2)}
+
 @end table