[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastrocommits] master b7d94d3: Table: new operator to calculate dist
From: 
Mohammad Akhlaghi 
Subject: 
[gnuastrocommits] master b7d94d3: Table: new operator to calculate distance on a flat surface 
Date: 
Fri, 15 May 2020 12:57:35 0400 (EDT) 
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 ('angulardistance'), 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 'distanceflat' operator. In order to have a unified
naming convention, the spherical distance operator's name was also changed
to 'distanceonsphere' (as mentioned above, originally it was
'angulardistance').

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.
 'dectodegree': Convert Declination (DD:MM:SS) to degrees.
 'degreetora': Convert degrees to Right Ascension (HH:MM:SS).
 'degreetodec': Convert degrees to Declination (HH:MM:SS).
+  'distanceflat': 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:
+  'distanceonsphere': New name for old `angulardistance' 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="angulardistance"; break;
+ case ARITHMETIC_TABLE_OP_DISTANCEFLAT: out="distanceflat"; break;
+ case ARITHMETIC_TABLE_OP_DISTANCEONSPHERE: out="distanceonsphere";
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, "angulardistance", 8))
 { op=ARITHMETIC_TABLE_OP_ANGULARDISTANCE; *num_operands=0; }
+ else if (!strncmp(string, "distanceflat", 13))
+ { op=ARITHMETIC_TABLE_OP_DISTANCEFLAT; *num_operands=0; }
+ else if (!strncmp(string, "distanceonsphere", 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=a1b1, d2=a2b2;
+ 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="distflat";
+ distance_func=arithmetic_distance_flat;
+ colcomment="Distance measured on a flat surface.";
+ break;
+ case ARITHMETIC_TABLE_OP_DISTANCEONSPHERE:
+ colname="distspherical";
+ 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, "angulardist",
 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 angulardistance
+@item distanceflat
+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 distanceflat'
+$ asttable table.fits \
+ c'arith X Y 12.345 6.789 distanceflat'
+$ asttable table.fits \
+ c'arith 12.345 6.789 X Y distanceflat'
+@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 distanceonsphere
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 angulardistance'
+ c'arith RA1 DEC1 RA2 DEC2 distanceonsphere'
$ asttable table.fits \
 c'arith RA DEC 12.345 6.789 angulardistance'
+ c'arith RA DEC 9.876 5.432 distanceonsphere'
$ asttable table.fits \
 c'arith 12.345 6.789 RA DEC angulardistance'
+ c'arith 9.876 5.432 RA DEC distanceonsphere'
@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.
[Prev in Thread] 
Current Thread 
[Next in Thread] 
 [gnuastrocommits] master b7d94d3: Table: new operator to calculate distance on a flat surface,
Mohammad Akhlaghi <=