[gnuastro-commits] master b7d94d3: Table: new operator to calculate distance on a flat surface

[Mohammad Akhlaghi]

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

    Table: new operator to calculate distance on a flat surface
    
    Until now we had a column arithmetic operator to find the distance between
    two points on a sphere ('angular-distance'), but it is often necessary to
    find the distance on a flat surface also.
    
    To avoid having to write complex code in such cases, with this commit,
    Table now has the 'distance-flat' operator. In order to have a unified
    naming convention, the spherical distance operator's name was also changed
    to 'distance-on-sphere' (as mentioned above, originally it was
    'angular-distance').
---
 NEWS                   |  3 +++
 bin/table/arithmetic.c | 59 ++++++++++++++++++++++++++++++++++++++++----------
 bin/table/arithmetic.h |  3 ++-
 doc/gnuastro.texi      | 38 +++++++++++++++++++++++++-------
 4 files changed, 82 insertions(+), 21 deletions(-)

diff --git a/NEWS b/NEWS
index c820ea8..923de68 100644
--- a/NEWS
+++ b/NEWS
@@ -60,6 +60,7 @@ See the end of the file for license conditions.
      - 'dec-to-degree': Convert Declination (DD:MM:SS) to degrees.
      - 'degree-to-ra': Convert degrees to Right Ascension (HH:MM:SS).
      - 'degree-to-dec': Convert degrees to Declination (HH:MM:SS).
+     - 'distance-flat': Distance between two points, assuming flat space.
 
   Library:
    - GAL_SPECLINES_INVALID_MAX: Total number of spectral lines, plus 1.
@@ -125,6 +126,8 @@ See the end of the file for license conditions.
      the new identifying character is very similar to AWK, allowing easier
      adoption and is also more clear. It is just important to put the total
      'arith' string within single quotes, not double quotes.
+   - Operators:
+     - 'distance-on-sphere': New name for old `angular-distance' operator.
 
   Library:
    - gal_polygon_is_inside_convex: new name for 'gal_polygon_pin'.
diff --git a/bin/table/arithmetic.c b/bin/table/arithmetic.c
index 3eae227..a4c6d69 100644
--- a/bin/table/arithmetic.c
+++ b/bin/table/arithmetic.c
@@ -106,7 +106,8 @@ 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;
+      case ARITHMETIC_TABLE_OP_DISTANCEFLAT: out="distance-flat"; break;
+      case ARITHMETIC_TABLE_OP_DISTANCEONSPHERE: out="distance-on-sphere"; 
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__,
@@ -157,8 +158,10 @@ arithmetic_set_operator(struct tableparams *p, char 
*string,
         { op=ARITHMETIC_TABLE_OP_WCSTOIMG; *num_operands=0; }
       else if (!strncmp(string, "imgtowcs", 8))
         { op=ARITHMETIC_TABLE_OP_IMGTOWCS; *num_operands=0; }
-      else if (!strncmp(string, "angular-distance", 8))
-        { op=ARITHMETIC_TABLE_OP_ANGULARDISTANCE; *num_operands=0; }
+      else if (!strncmp(string, "distance-flat", 13))
+        { op=ARITHMETIC_TABLE_OP_DISTANCEFLAT; *num_operands=0; }
+      else if (!strncmp(string, "distance-on-sphere", 18))
+        { op=ARITHMETIC_TABLE_OP_DISTANCEONSPHERE; *num_operands=0; }
       else
         { op=GAL_ARITHMETIC_OP_INVALID; *num_operands=GAL_BLANK_INT; }
     }
@@ -433,12 +436,25 @@ arithmetic_wcs(struct tableparams *p, gal_data_t **stack, 
int operator)
 
 
 
+static double
+arithmetic_distance_flat(double a1, double a2, double b1, double b2)
+{
+  double d1=a1-b1, d2=a2-b2;
+  return sqrt(d1*d1 + d2*d2);
+}
+
+
+
+
+
 static void
-arithmetic_angular_dist(struct tableparams *p, gal_data_t **stack, int 
operator)
+arithmetic_distance(struct tableparams *p, gal_data_t **stack, int operator)
 {
   size_t i, j;
+  char *colname, *colcomment;
   double *o, *a1, *a2, *b1, *b2;
   gal_data_t *a, *b, *tmp, *out;
+  double (*distance_func)(double, double, double, double);
 
   /* Pop the columns for point 'b'.*/
   tmp=arithmetic_stack_pop(stack, operator);
@@ -462,16 +478,35 @@ arithmetic_angular_dist(struct tableparams *p, gal_data_t 
**stack, int operator)
           "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 "
+    error(EXIT_FAILURE, 0, "the sizes of the first and second operands "
           "of the '%s' operator (respectively containing %zu and %zu "
           "numbers) must be equal", arithmetic_operator_name(operator),
           b->next->size, b->size);
 
+  /* Basic settings based on the operator. */
+  switch(operator)
+    {
+    case ARITHMETIC_TABLE_OP_DISTANCEFLAT:
+      colname="dist-flat";
+      distance_func=arithmetic_distance_flat;
+      colcomment="Distance measured on a flat surface.";
+      break;
+    case ARITHMETIC_TABLE_OP_DISTANCEONSPHERE:
+      colname="dist-spherical";
+      distance_func=gal_wcs_angular_distance_deg;
+      colcomment="Distance measured on a great circle.";
+      break;
+    default:
+      error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+            "fix the problem. The operator code %d isn't recognized",
+            __func__, PACKAGE_BUGREPORT, operator);
+    }
+
   /* 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");
+                     p->cp.minmapsize, p->cp.quietmmap, colname, NULL,
+                     colcomment);
 
   /* Measure the distances.  */
   o=out->array;
@@ -480,11 +515,10 @@ arithmetic_angular_dist(struct tableparams *p, gal_data_t 
**stack, int operator)
   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]);
+        o[a->size>b->size?i:j] = distance_func(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]);
+      o[i] = distance_func(a1[i], a2[i], b1[i], b2[i]);
 
   /* Clean up and put the output dataset onto the stack. */
   gal_list_data_free(a);
@@ -610,8 +644,9 @@ 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);
+        case ARITHMETIC_TABLE_OP_DISTANCEFLAT:
+        case ARITHMETIC_TABLE_OP_DISTANCEONSPHERE:
+          arithmetic_distance(p, stack, operator);
           break;
 
         default:
diff --git a/bin/table/arithmetic.h b/bin/table/arithmetic.h
index 2df0278..8128f92 100644
--- a/bin/table/arithmetic.h
+++ b/bin/table/arithmetic.h
@@ -37,7 +37,8 @@ enum arithmetic_operators
 {
  ARITHMETIC_TABLE_OP_WCSTOIMG = GAL_ARITHMETIC_OP_LAST_CODE,
  ARITHMETIC_TABLE_OP_IMGTOWCS,
- ARITHMETIC_TABLE_OP_ANGULARDISTANCE,
+ ARITHMETIC_TABLE_OP_DISTANCEFLAT,
+ ARITHMETIC_TABLE_OP_DISTANCEONSPHERE,
 };
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 29b38ac..e3748a4 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -9172,29 +9172,51 @@ $ 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
+@item distance-flat
+Return the distance between two points assuming they are on a flat surface.
+Note that each point needs two coordinates, so this operator needs four 
operands (currently it only works for 2D spaces).
+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 X1 Y1 X2 Y2 distance-flat'
+$ asttable table.fits \
+           -c'arith X Y 12.345 6.789 distance-flat'
+$ asttable table.fits \
+           -c'arith 12.345 6.789 X Y distance-flat'
+@end example
+
+In the first case we are assuming that @file{table.fits} has the following 
four columns @code{X1}, @code{Y1}, @code{X2}, @code{Y2}.
+The returned column by this operator will be the difference between two points 
in each row with coordinates like the following (@code{X1}, @code{Y1}) and 
(@code{X2}, @code{Y2}).
+In other words, for each row, the 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{X} and @code{Y}.
+The returned column by this operator will be the difference of each row with 
the fixed point at (12.345, 6.789).
+
+@item distance-on-sphere
 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.
+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'
+           -c'arith RA1 DEC1 RA2 DEC2 distance-on-sphere'
 $ asttable table.fits \
-           -c'arith RA DEC 12.345 6.789 angular-distance'
+           -c'arith RA DEC 9.876 5.432 distance-on-sphere'
 $ asttable table.fits \
-           -c'arith 12.345 6.789 RA DEC angular-distance'
+           -c'arith 9.876 5.432 RA DEC distance-on-sphere'
 @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}).
+The returned column by this operator will 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 returned column by this operator will be the difference of each row with 
the fixed point at (9.876, 5.432).
 
 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.