[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master a9656bcd 1/2: Library (arithmetic.h): new box-
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master a9656bcd 1/2: Library (arithmetic.h): new box-vertices-on-sphere operator |
|
Date: |
Thu, 22 Dec 2022 12:43:38 -0500 (EST) |
branch: master
commit a9656bcdf3868a2e48bd0a4af99b8f945935943a
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Library (arithmetic.h): new box-vertices-on-sphere operator
Until now, in several places within Gnuastro we accounted for the fact that
a box's vertices on a sphere need special consideration. However, for an
outside user, doing the job manually would require a complex column
arithmetic operation (and to get to pen-and-paper to derive the
equations!).
With this commit, a simple arithmetic operator has been added for this job:
'box-vertices-on-sphere'. It will take the center and size of the box along
each dimension and will return the coordinates of the 8 corners, while
accounting for the spherical effects. Some real-world usage examples and a
description of the problem it is designed to solve have been given in the
extensive documentation of this new operator.
This was added after a discussion with Martin Kuemmel.
---
NEWS | 8 ++
bin/arithmetic/arithmetic.c | 11 ++-
bin/table/arithmetic.c | 11 ++-
doc/announce-acknowledge.txt | 1 +
doc/gnuastro.texi | 133 ++++++++++++++++++++++++---
lib/arithmetic.c | 208 +++++++++++++++++++++++++++++++++----------
lib/gnuastro/arithmetic.h | 1 +
lib/gnuastro/wcs.h | 5 ++
lib/wcs.c | 38 ++++++++
9 files changed, 350 insertions(+), 66 deletions(-)
diff --git a/NEWS b/NEWS
index 6fc6f1f3..52dc2fcc 100644
--- a/NEWS
+++ b/NEWS
@@ -30,6 +30,10 @@ See the end of the file for license conditions.
SDSS). This was suggested by Giulia Golini.
- nanomaggy-to-counts: convert nanomaggy to counts. This was suggested
by Giulia Golini.
+ - box-vertices-on-sphere: calculate the RA,Dec of vertices of the four
+ vertices of a rectangle from its center and width along RA and
+ Dec. This takes into account the curved nature of the coordinate
+ system. Added after discussion with Martin Kuemmel.
- New operators (only in the Arithmetic program):
- interpolate-meanngb: interpolate blank values with mean of the
requested number of nearest neighbors.
@@ -120,10 +124,14 @@ See the end of the file for license conditions.
- GAL_ARITHMETIC_OP_COUNTERONLY: Similar to 'GAL_ARITHMETIC_OP_COUNTER'.
- GAL_ARITHMETIC_OP_COUNTS_TO_NANOMAGGY: convert counts to nanomaggy.
- GAL_ARITHMETIC_OP_NANOMAGGY_TO_COUNTS: convert nanomaggy to counts.
+ - GAL_ARITHMETIC_OP_BOX_VERTICES_ON_SPHERE: calculate the coordinates of
+ vertices of a rectable on a sphere from its center and width/height.
- gal_units_counts_to_nanomaggy: Convert counts to nanomaggy.
- gal_units_nanomaggy_to_counts: Convert nanomaggy to counts.
- gal_data_alloc_empty: Allocate an empty dataset with a given number of
dimensions.
+ - gal_wcs_box_vertices_from_center: calculate the coordinates of
+ vertices of a rectable on a sphere from its center and width/height.
** Removed features
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 7c6cbdfe..c3f49023 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -1414,7 +1414,7 @@ arithmetic_operator_run(struct arithmeticparams *p, int
operator,
size_t i;
unsigned int numop;
int flags = GAL_ARITHMETIC_FLAGS_BASIC;
- gal_data_t *d1=NULL, *d2=NULL, *d3=NULL;
+ gal_data_t *d1=NULL, *d2=NULL, *d3=NULL, *d4=NULL;
/* Set the operating-mode flags if necessary. */
if(p->cp.quiet) flags |= GAL_ARITHMETIC_FLAG_QUIET;
@@ -1448,6 +1448,13 @@ arithmetic_operator_run(struct arithmeticparams *p, int
operator,
d1=operands_pop(p, operator_string);
break;
+ case 4:
+ d4=operands_pop(p, operator_string);
+ d3=operands_pop(p, operator_string);
+ d2=operands_pop(p, operator_string);
+ d1=operands_pop(p, operator_string);
+ break;
+
case -1:
/* This case is when the number of operands is itself an
operand. So except for operators that have high-level
@@ -1486,7 +1493,7 @@ arithmetic_operator_run(struct arithmeticparams *p, int
operator,
when the operator doesn't need three operands, the extra
arguments will be ignored. */
operands_add(p, NULL, gal_arithmetic(operator, p->cp.numthreads,
- flags, d1, d2, d3));
+ flags, d1, d2, d3, d4));
}
/* No need to call the arithmetic library, call the proper wrappers
diff --git a/bin/table/arithmetic.c b/bin/table/arithmetic.c
index ac9f0636..3246e64e 100644
--- a/bin/table/arithmetic.c
+++ b/bin/table/arithmetic.c
@@ -932,7 +932,7 @@ arithmetic_operator_run(struct tableparams *p,
gal_data_t **stack)
{
int flags=GAL_ARITHMETIC_FLAGS_BASIC;
- gal_data_t *d1=NULL, *d2=NULL, *d3=NULL;
+ gal_data_t *d1=NULL, *d2=NULL, *d3=NULL, *d4=NULL;
/* Set the operating-mode flags if necessary. */
if(p->cp.quiet) flags |= GAL_ARITHMETIC_FLAG_QUIET;
@@ -966,6 +966,13 @@ arithmetic_operator_run(struct tableparams *p,
d1=arithmetic_stack_pop(stack, token->operator, NULL);
break;
+ case 4:
+ d4=arithmetic_stack_pop(stack, token->operator, NULL);
+ d3=arithmetic_stack_pop(stack, token->operator, NULL);
+ d2=arithmetic_stack_pop(stack, token->operator, NULL);
+ d1=arithmetic_stack_pop(stack, token->operator, NULL);
+ break;
+
case -1:
error(EXIT_FAILURE, 0, "operators with a variable number of "
"operands are not yet implemented. Please contact us at "
@@ -987,7 +994,7 @@ arithmetic_operator_run(struct tableparams *p,
ignored. */
gal_list_data_add(stack, gal_arithmetic(token->operator,
p->cp.numthreads,
- flags, d1, d2, d3) );
+ flags, d1, d2, d3, d4) );
/* Reset the meta-data for the element that was just put on the
stack. */
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index fbe4e7db..867d2eed 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -5,6 +5,7 @@ Alejandro Serrano Borlaff
Sepideh Eskandarlou
Zohreh Ghaffari
Giulia Golini
+Martin Kuemmel
Samane Raji
Elham Saremi
Nafiseh Sedighi
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 26441f01..6d594f1e 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -234,7 +234,6 @@ To go to the sections or subsections, you have to click on
the menu entries that
* Makefile extensions:: Use Gnuastro's features as GNU Make functions.
* Library:: Gnuastro's library of useful functions.
* Developing:: The development environment.
-* Gnuastro programs list:: List and short summary of Gnuastro.
* Other useful software:: Installing other useful software.
* GNU Free Doc License:: Full text of the GNU Free documentation
license.
* GNU General Public License:: Full text of the GNU General public license.
@@ -246,6 +245,7 @@ To go to the sections or subsections, you have to click on
the menu entries that
Introduction
* Quick start:: A quick start to installation.
+* Gnuastro programs list::
* Science and its tools:: Some philosophy and history.
* Your rights:: User rights.
* Logo of Gnuastro:: Meaning behind Gnuastro's logo.
@@ -521,7 +521,7 @@ Arithmetic operators
* Bitwise operators:: Work on bits within one pixel.
* Numerical type conversion operators:: Convert the numeric datatype of a
dataset.
* Random number generators:: Random numbers can be used to add noise for
example.
-* Elliptical shape operators:: Operations that are focused on an ellipse.
+* Box shape operators:: Dealing with box shapes and coordinates of
vertices.
* Loading external columns:: Read a column from a table into the stack.
* Size and position operators:: Extracting image size and pixel positions.
* Building new dataset and stack management:: How to construct an empty
dataset from scratch.
@@ -921,6 +921,7 @@ As discussed in @ref{Science and its tools} this is a
founding principle of the
@menu
* Quick start:: A quick start to installation.
+* Gnuastro programs list::
* Science and its tools:: Some philosophy and history.
* Your rights:: User rights.
* Logo of Gnuastro:: Meaning behind Gnuastro's logo.
@@ -934,7 +935,7 @@ As discussed in @ref{Science and its tools} this is a
founding principle of the
* Acknowledgments:: People who helped in the production.
@end menu
-@node Quick start, Science and its tools, Introduction, Introduction
+@node Quick start, Gnuastro programs list, Introduction, Introduction
@section Quick start
@cindex Test
@@ -1000,7 +1001,7 @@ We therefore recommend to read (an run the commands in)
the tutorials before sta
-@node Gnuastro programs list, Other useful software, Developing, Top
+@node Gnuastro programs list, Science and its tools, Quick start, Introduction
@section Gnuastro programs list
GNU Astronomy Utilities @value{VERSION}, contains the following programs.
@@ -1125,7 +1126,7 @@ Because of saturation and non-linearity, to get a good
estimate of the extended
-@node Science and its tools, Your rights, Quick start, Introduction
+@node Science and its tools, Your rights, Gnuastro programs list, Introduction
@section Gnuastro manifesto: Science and its tools
History of science indicates that there are always inevitably unseen faults,
hidden assumptions, simplifications and approximations in all our theoretical
models, data acquisition and analysis techniques.
@@ -16665,7 +16666,7 @@ Reading NaN as a floating point number in Gnuastro is
not case-sensitive.
* Bitwise operators:: Work on bits within one pixel.
* Numerical type conversion operators:: Convert the numeric datatype of a
dataset.
* Random number generators:: Random numbers can be used to add noise for
example.
-* Elliptical shape operators:: Operations that are focused on an ellipse.
+* Box shape operators:: Dealing with box shapes and coordinates of
vertices.
* Loading external columns:: Read a column from a table into the stack.
* Size and position operators:: Extracting image size and pixel positions.
* Building new dataset and stack management:: How to construct an empty
dataset from scratch.
@@ -18036,7 +18037,7 @@ Convert the type of the popped operand to 64-bit
(double precision) floating poi
The internal conversion of C will be used.
@end table
-@node Random number generators, Elliptical shape operators, Numerical type
conversion operators, Arithmetic operators
+@node Random number generators, Box shape operators, Numerical type conversion
operators, Arithmetic operators
@subsubsection Random number generators
When you simulate data (for example, see @ref{Sufi simulates a detection}),
everything is ideal and there is no noise!
@@ -18322,10 +18323,10 @@ $ echo 1000 \
These columns can easily be placed in the format for @ref{MakeProfiles} to be
inserted into an image automatically.
@end table
-@node Elliptical shape operators, Loading external columns, Random number
generators, Arithmetic operators
-@subsubsection Elliptical shape operators
+@node Box shape operators, Loading external columns, Random number generators,
Arithmetic operators
+@subsubsection Box shape operators
-The operators here describe certain functions that will be necessary when
dealing with objects that have a certain elliptical shape.
+The operators here help you in defining or using coordinates that form a
``box'' (a rectangular region).
@table @command
@item box-around-ellipse
@@ -18368,9 +18369,83 @@ $ asttable catalog.fits \
set-height set-width \
width float32 height float32'
@end example
+
+@item box-vertices-on-sphere
+@cindex Polygon
+@cindex Vertices on sphere (sky)
+Convert a box center and width to the coordinates of the vertices of the box
on a left-hand spherical coordinate system.
+In a left-handed spherical coordinate system, the longitude increases towards
the left while north is up (as in the RA and Dec direction of the equatorial
coordinate system used in astronomy).
+This operator therefore takes four input operands (the RA and Dec of the box's
center, as well as the width of the box in each direction).
+
+After it is complete, this operator places 8 operands on the stack which
contain the RA and Dec of the four vertices of the box in the following
anti-clockwise order:
+@enumerate
+@item
+Bottom-left vertice Logitude (RA)
+@item
+Bottom-left vertice Latitude (Dec)
+@item
+Bottom-right vertice Logitude (RA)
+@item
+Bottom-right vertice Latitude (Dec)
+@item
+Top-right vertice Logitude (RA)
+@item
+Top-right vertice Latitude (Dec)
+@item
+Top-left vertice Logitude (RA)
+@item
+Top-left vertice Latitude (Dec)
+@end enumerate
+
+For example, with the command below, we will retrieve the vertice coordinates
of a rectangle around a point with RA=20 and Dec=0 (on the equator).
+The rectangle will have a 1 degree edge along the RA direction and a 2 degree
edge along the declination.
+In this example, we are using the @option{-Afixed -B2} only for demonstration
purposes here due to the round numbers!
+In general, it is best to write your outputs to a binary FITS table to
preserve the full precision (see @ref{Printing floating point numbers}).
+
+@example
+$ echo "20 0 1 2" \
+ | asttable -Afixed -B2 \
+ -c'arith $1 $2 $3 $4 box-vertices-on-sphere'
+20.50 -1.00 19.50 -1.00 19.50 1.00 20.50 1.00
+@end example
+
+We see that the bottom-left vertice is at (RA,Dec) of @mymath{(20.50,-1.0)}
and the top-right vertice is at @mymath{(19.50,1.00)}.
+These could have easily been done by manually adding and subtracting!
+But you will see that the complextiy arises at higher/lower declinations.
+For example, with the command below, let's see how vertice coordinates of the
same box, but after moving its center to (RA,Dec) of (20,85):
+
+@example
+$ echo "20 85 1 2" \
+ | asttable -Afixed -B2 \
+ -c'arith $1 $2 $3 $4 box-vertices-on-sphere'
+24.78 84.00 15.22 84.00 12.83 86.00 27.17 86.00
+@end example
+
+Even though, we didn't change the central RA (20) or the size of the box
along the RA (1 degree), the RA of the bottom-left vertice is now at 24.78;
almost 5 degrees away!
+This occurs because of the spherical coordinate system, we measure the
longitude (e.g., RA) with the following way:
+@enumerate
+@item
+@cindex Meridian
+@cindex Great circle
+@cindex Circle (great)
+Draw a meridian that passes your point.
+The meridian is half of a @url{https://en.wikipedia.org/wiki/Great_circle,
great-circle} (which has a diameter that is equal to the sphere's diameter)
passes both poles.
+@item
+Find the intersection of that meridian with the equator.
+@item
+The distance of the intersection and the reference point (along the equator)
defines the longitude angle.
+@end enumerate
+
+@cindex Small circle
+@cindex Circle (small)
+As you get more distant from the equator (declination becomes non-zero), any
change along the RA (towards the east; 1 degree in the example above) will on
longer be on a great circle, but along a
``@url{https://en.wikipedia.org/wiki/Circle_of_a_sphere, small circle}''.
+On a small circle that is defined by the fixed declination @mymath{\delta},
the distance of two points is closer than the distances of their projection on
the equator (as described in the definition of longitude above).
+It is smaller by a factor of @mymath{\cos({\delta})}.
+
+Therefore, an angular change (let's call it @mymath{\Delta_{lon}}) along the
small circle defined by the fixed declination of @mymath{\delta} corresponds to
@mymath{\Delta_{lon}/\cos(\delta)} on the equator.
@end table
-@node Loading external columns, Size and position operators, Elliptical shape
operators, Arithmetic operators
+@node Loading external columns, Size and position operators, Box shape
operators, Arithmetic operators
@subsubsection Loading external columns
In the Arithmetic program, you can always load new dataset by simply giving
their name.
@@ -34953,6 +35028,32 @@ It will be use the
@url{https://en.wikipedia.org/wiki/Haversine_formula, Haversi
@dispmath{{\sin^2(d)\over 2}=\sin^2\left( {d_1-d_2\over 2}
\right)+\cos(d_1)\cos(d_2)\sin^2\left( {r_1-r_2\over 2} \right)}
@end deftypefun
+@deftypefun void gal_wcs_box_vertices_from_center (double @code{ra_center},
double @code{dec_center}, double @code{ra_delta}, double @code{dec_delta},
double @code{*out})
+Calculate the vertices of a rectangular box given the centeral RA and Dec and
delta of each.
+The vertice coordinates are written in the space that @code{out} points to
(assuming it has space for eight @code{double}s).
+
+Given the spherical nature of the coordinate system, the vertice lengths can't
be calculated with a simple addition/subtraction.
+For the declination, a simple addition/subtraction is enough.
+Also, on the equator (where the RA is defined), a simple addition/subtraction
along the RA is fine.
+However, at other declinations, the new RA after a shift needs special
treatment, such that close to the poles, a shift of 1 degree can correspond to
a new RA that is much more distant than the original RA.
+Assuming a point at Right Ascension (RA) and Declination of @mymath{\alpha}
and @mymath{\delta}, a shift of @mymath{R} degrees along the positive RA
direction corresponds to a right ascension of
@mymath{\alpha+\frac{R}{\cos(\delta)}}.
+For more, see the description of @code{box-vertices-on-sphere} in @ref{Box
shape operators}.
+
+The 8 coordinates of the 4 vertices of the box are written in the order below.
+Where ``bottom'' corresponds to a lower declination and ``top'' to higher
declination, ``left'' corresponds to a larger RA and ``right'' corresponds to
lower RA.
+
+@example
+out[0]: bottom-left RA
+out[1]: bottom-left Dec
+out[2]: bottom-right RA
+out[3]: bottom-right Dec
+out[4]: top-right RA
+out[5]: top-right Dec
+out[6]: top-left RA
+out[7]: top-left Dec
+@end example
+@end deftypefun
+
@deftypefun {double *} gal_wcs_pixel_scale (struct wcsprm @code{*wcs})
Return the pixel scale for each dimension of @code{wcs} in degrees.
The output is an allocated array of double precision floating point type with
one element for each dimension.
@@ -35324,6 +35425,12 @@ This function returns two datasets as a
@code{gal_data_t} linked list.
The top element of the list is the height and its next element is the width.
@end deffn
+@deffn Macro GAL_ARITHMETIC_OP_BOX_VERTICES_ON_SPHERE
+Return the vertices of a (possibly rectangluar) box on a sphere, given its
center RA, Dec and the width of the box along the two dimensions.
+It will take the spherical nature of the coordinate system into account (for
more, see the description of @code{gal_wcs_box_vertices_from_center} in
@ref{World Coordinate System}).
+This function returns 8 datasets as a @code{gal_data_t} linked list in the
following order: bottom-left RA, bottom-left Dec, bottom-right RA, bottom-right
Dec, top-right RA, top-right Dec, top-left RA, top-left Dec.
+@end deffn
+
@deffn Macro GAL_ARITHMETIC_OP_MAKENEW
Create a new, zero-valued dataset with an unsigned 8-bit data type.
The length along each dimension of the dataset should be given as a single
list of @code{gal_data_t}s.
@@ -39564,7 +39671,7 @@ main(void)
-@node Developing, Gnuastro programs list, Library, Top
+@node Developing, Other useful software, Library, Top
@chapter Developing
The basic idea of GNU Astronomy Utilities is for an interested astronomer
@@ -41448,7 +41555,7 @@ This will enable you to start your branch (work) from
the most recent commit and
-@node Other useful software, GNU Free Doc License, Gnuastro programs list, Top
+@node Other useful software, GNU Free Doc License, Developing, Top
@appendix Other useful software
In this appendix the installation of programs and libraries that are
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 2842f745..4a9ee667 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -37,6 +37,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gsl/gsl_const_mks.h>
#include <gnuastro/box.h>
+#include <gnuastro/wcs.h>
#include <gnuastro/list.h>
#include <gnuastro/blank.h>
#include <gnuastro/units.h>
@@ -2540,49 +2541,13 @@ arithmetic_counts_to_from_sb(int operator, int flags,
gal_data_t *d1,
static gal_data_t *
-arithmetic_box_around_ellipse(gal_data_t *d1, gal_data_t *d2,
- gal_data_t *d3, int operator, int flags)
+arithmetic_box_around_ellipse(gal_data_t *a_data, gal_data_t *b_data,
+ gal_data_t *pa_data, int flags)
{
size_t i;
- double *a, *b, *pa, out[2];
- gal_data_t *a_data, *b_data, *pa_data;
+ double *a=a_data->array, *b=b_data->array, *pa=pa_data->array, out[2];
- /* Basic sanity check. */
- if( d1->size != d2->size || d1->size != d3->size )
- error(EXIT_FAILURE, 0, "%s: the three operands to this "
- "function don't have the same number of elements",
- __func__);
-
- /* The datasets may be empty. In this case the output should also be
- empty (we can have tables and images with 0 rows or pixels!). */
- if( d1->size==0 || d1->array==NULL
- || d2->size==0 || d2->array==NULL
- || d3->size==0 || d3->array==NULL )
- {
- if(flags & GAL_ARITHMETIC_FLAG_FREE)
- { gal_data_free(d2); gal_data_free(d3); }
- if(d1->array) {free(d1->array); d1->array=NULL;}
- if(d1->dsize) for(i=0;i<d1->ndim;++i) d1->dsize[i]=0;
- d1->size=0; return d1;
- }
-
- /* Convert the inputs into double. Note that if the user doesn't want to
- free the inputs, we should make a copy of 'a_data' and 'b_data'
- because the output will also be written in them. */
- a_data=( d1->type==GAL_TYPE_FLOAT64
- ? d1
- : gal_data_copy_to_new_type(d1, GAL_TYPE_FLOAT64) );
- b_data=( d2->type==GAL_TYPE_FLOAT64
- ? d2
- : gal_data_copy_to_new_type(d2, GAL_TYPE_FLOAT64) );
- pa_data=( d3->type==GAL_TYPE_FLOAT64
- ? d3
- : gal_data_copy_to_new_type(d3, GAL_TYPE_FLOAT64) );
-
- /* Set the double pointers and do the calculation for each. */
- a=a_data->array;
- b=b_data->array;
- pa=pa_data->array;
+ /* Loop over all the elements. */
for(i=0;i<a_data->size;++i)
{
/* If the minor axis has a larger value, print a warning and set the
@@ -2603,18 +2568,153 @@ arithmetic_box_around_ellipse(gal_data_t *d1,
gal_data_t *d2,
b[i]=out[1];
}
+ /* Make the output as a list and return it. */
+ b_data->next=a_data;
+ return b_data;
+}
+
+
+
+
+
+/* Calculate the vertices of a box from its center and width. The longitude
+ coordinate/width is specified with a 'lon' and the latitude
+ coordinate/width is specified with a 'lat'. */
+static gal_data_t *
+arithmetic_box_vertices_on_sphere(gal_data_t *d1, gal_data_t *d2,
+ gal_data_t *d3, gal_data_t *d4,
+ int flags)
+{
+ size_t i;
+ double vertices[8];
+ double *clon=d1->array, *clat=d2->array;
+ double *dlon=d3->array, *dlat=d4->array;
+
+ /* Output datasets. */
+ double *blr, *bld, *brr, *brd, *tlr, *tld, *trr, *trd;
+ gal_data_t *blrd, *bldd, *brrd, *brdd, *tlrd, *tldd, *trrd, *trdd;
+
+ /* Allocate the extra output datasets (four of them will be the same as
+ the input (which will be over-written after usage) and set the array
+ pointers. */
+ blrd=d1; bldd=d2; brrd=d3; brdd=d4;
+ tlrd=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, d1->ndim, d1->dsize, NULL,
+ 0, d1->minmapsize, d1->quietmmap, NULL, NULL, NULL);
+ tldd=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, d1->ndim, d1->dsize, NULL,
+ 0, d1->minmapsize, d1->quietmmap, NULL, NULL, NULL);
+ trrd=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, d1->ndim, d1->dsize, NULL,
+ 0, d1->minmapsize, d1->quietmmap, NULL, NULL, NULL);
+ trdd=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, d1->ndim, d1->dsize, NULL,
+ 0, d1->minmapsize, d1->quietmmap, NULL, NULL, NULL);
+ blr=blrd->array; bld=bldd->array; brr=brrd->array; brd=brdd->array;
+ tlr=tlrd->array; tld=tldd->array; trr=trrd->array; trd=trdd->array;
+
+ /* Loop over the elements. */
+ for(i=0;i<d1->size;++i)
+ {
+ /* Compute the positions of the vertices. */
+ gal_wcs_box_vertices_from_center(clon[i], clat[i], dlon[i],
+ dlat[i], vertices);
+
+ /* Write the results into the the output arrays. */
+ blr[i]=vertices[0]; bld[i]=vertices[1];
+ brr[i]=vertices[2]; brd[i]=vertices[3];
+ tlr[i]=vertices[4]; tld[i]=vertices[5];
+ trr[i]=vertices[6]; trd[i]=vertices[7];
+ }
+
+ /* Create the output list, note that it has to in the inverse order of
+ the final output order. In the end, we want the order of the corners
+ to be bottom_left_ra -> bottom_left_dec -> bottom_right_ra ->
+ bottom_right_dec -> top_right_ra -> top_right_dec -> top_left_ra ->
+ top_left_dec. */
+ blrd->next=bldd;
+ bldd->next=brrd;
+ brrd->next=brdd;
+ brdd->next=trrd;
+ trrd->next=trdd;
+ trdd->next=tlrd;
+ tlrd->next=tldd;
+ gal_list_data_reverse(&blrd);
+ return blrd;
+}
+
+
+
+
+static gal_data_t *
+arithmetic_box(gal_data_t *d1, gal_data_t *d2, gal_data_t *d3,
+ gal_data_t *d4, int operator, int flags)
+{
+ size_t i;
+ gal_data_t *out=NULL;
+ gal_data_t *d1d=NULL, *d2d=NULL, *d3d=NULL, *d4d=NULL;
+
+ /* Basic sanity check. */
+ if( d1->size != d2->size || d1->size != d3->size
+ || (d4 && d1->size!=d4->size) )
+ error(EXIT_FAILURE, 0, "%s: the operands to this function "
+ "don't have the same number of elements", __func__);
+
+ /* The datasets may be empty. In this case the output should also be
+ empty (we can have tables and images with 0 rows or pixels!). */
+ if( d1->size==0 || d1->array==NULL
+ || d2->size==0 || d2->array==NULL
+ || d3->size==0 || d3->array==NULL
+ || (d4 && (d4->size==0 || d4->array==NULL) ) )
+ {
+ if(flags & GAL_ARITHMETIC_FLAG_FREE)
+ { gal_data_free(d2); gal_data_free(d3); gal_data_free(d4); }
+ if(d1->array) {free(d1->array); d1->array=NULL;}
+ if(d1->dsize) for(i=0;i<d1->ndim;++i) d1->dsize[i]=0;
+ d1->size=0; return d1;
+ }
+
+ /* Convert the inputs into double. Note that if the user doesn't want to
+ free the inputs, we should make a copy of 'a_data' and 'b_data' because
+ the output will also be written in them. */
+ d1d=( d1->type==GAL_TYPE_FLOAT64
+ ? d1
+ : gal_data_copy_to_new_type(d1, GAL_TYPE_FLOAT64) );
+ d2d=( d2->type==GAL_TYPE_FLOAT64
+ ? d2
+ : gal_data_copy_to_new_type(d2, GAL_TYPE_FLOAT64) );
+ d3d=( d3->type==GAL_TYPE_FLOAT64
+ ? d3
+ : gal_data_copy_to_new_type(d3, GAL_TYPE_FLOAT64) );
+ if(d4)
+ d4d=( d4->type==GAL_TYPE_FLOAT64
+ ? d4
+ : gal_data_copy_to_new_type(d4, GAL_TYPE_FLOAT64) );
+
+ /* Call the appropriate function. */
+ switch(operator)
+ {
+ case GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE:
+ out=arithmetic_box_around_ellipse(d1d, d2d, d3d, flags);
+ break;
+ case GAL_ARITHMETIC_OP_BOX_VERTICES_ON_SPHERE:
+ out=arithmetic_box_vertices_on_sphere(d1d, d2d, d3d, d4d, flags);
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+ "fix the problem. The code '%d' is not a recognized "
+ "operator for this function", __func__, PACKAGE_BUGREPORT,
+ operator);
+ }
+
/* Clean up. */
- gal_data_free(pa_data);
if(flags & GAL_ARITHMETIC_FLAG_FREE)
{
- if(a_data !=d1) gal_data_free(d1);
- if(b_data !=d2) gal_data_free(d2);
- if(pa_data!=d3) gal_data_free(d3);
+ if(d1d!=d1) gal_data_free(d1);
+ if(d2d!=d2) gal_data_free(d2);
+ if(d3d!=d3) gal_data_free(d3);
}
+ if(operator==GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE)
+ { gal_data_free(d3d); gal_data_free(d4d); }
- /* Make the output as a list and return it. */
- b_data->next=a_data;
- return b_data;
+ /* Return the output. */
+ return out;
}
@@ -3221,6 +3321,8 @@ gal_arithmetic_set_operator(char *string, size_t
*num_operands)
/* Surrounding box. */
else if (!strcmp(string, "box-around-ellipse"))
{ op=GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE;*num_operands=3; }
+ else if (!strcmp(string, "box-vertices-on-sphere"))
+ { op=GAL_ARITHMETIC_OP_BOX_VERTICES_ON_SPHERE; *num_operands=4; }
/* Size and position operators. */
else if (!strcmp(string, "swap"))
@@ -3378,6 +3480,7 @@ gal_arithmetic_operator_string(int operator)
case GAL_ARITHMETIC_OP_FINESTRUCTURE: return "fine-structure";
case GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE: return "box-around-ellipse";
+ case GAL_ARITHMETIC_OP_BOX_VERTICES_ON_SPHERE: return "vertices-on-sphere";
case GAL_ARITHMETIC_OP_SWAP: return "swap";
case GAL_ARITHMETIC_OP_INDEX: return "index";
@@ -3400,7 +3503,7 @@ gal_data_t *
gal_arithmetic(int operator, size_t numthreads, int flags, ...)
{
va_list va;
- gal_data_t *d1, *d2, *d3, *out=NULL;
+ gal_data_t *d1, *d2, *d3, *d4, *out=NULL;
/* Prepare the variable arguments (starting after the flags argument). */
va_start(va, flags);
@@ -3618,10 +3721,17 @@ gal_arithmetic(int operator, size_t numthreads, int
flags, ...)
/* Calculate the width and height of a box surrounding an ellipse with
a certain major axis, minor axis and position angle. */
case GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE:
+ case GAL_ARITHMETIC_OP_BOX_VERTICES_ON_SPHERE:
d1 = va_arg(va, gal_data_t *);
d2 = va_arg(va, gal_data_t *);
d3 = va_arg(va, gal_data_t *);
- out=arithmetic_box_around_ellipse(d1, d2, d3, operator, flags);
+ if(operator==GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE)
+ out=arithmetic_box(d1, d2, d3, NULL, operator, flags);
+ else
+ {
+ d4=va_arg(va, gal_data_t *);
+ out=arithmetic_box(d1, d2, d3, d4, operator, flags);
+ }
break;
/* Size and position operators. */
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index 7f5212c3..7c2b5b7b 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -207,6 +207,7 @@ enum gal_arithmetic_operators
GAL_ARITHMETIC_OP_TO_FLOAT64, /* Convert to float64. */
GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE, /* Width/Height of box over ellipse*/
+ GAL_ARITHMETIC_OP_BOX_VERTICES_ON_SPHERE, /* Vert. from center and width*/
/* Meta operators */
GAL_ARITHMETIC_OP_MAKENEW, /* Build a new dataset, containing zeros.*/
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 6523f83c..6c0f2a75 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -197,6 +197,11 @@ gal_wcs_to_cd(struct wcsprm *wcs);
double
gal_wcs_angular_distance_deg(double r1, double d1, double r2, double d2);
+void
+gal_wcs_box_vertices_from_center(double ra_center, double dec_center,
+ double ra_delta, double dec_delta,
+ double *out);
+
double *
gal_wcs_pixel_scale(struct wcsprm *wcs);
diff --git a/lib/wcs.c b/lib/wcs.c
index 0aef1ab7..60659b9d 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -1857,6 +1857,44 @@ gal_wcs_angular_distance_deg(double r1, double d1,
double r2, double d2)
+/* Calculate the vertices of a box in the given center. The output should
+ have 8 allocated spaces ready to be written into:
+ out[0]=bottom-left-ra
+ out[1]=bottom-left-dec
+ out[2]=bottom-right-ra
+ out[3]=bottom-right-dec
+ out[4]=top-right-ra
+ out[5]=top-right-dec
+ out[6]=top-left-ra
+ out[7]=top-left-dec
+
+ The 'ra_delta' and 'dec_delta' should be the full length of the box along
+ the respective axis (not half of it!). */
+void
+gal_wcs_box_vertices_from_center(double ra_center, double dec_center,
+ double ra_delta, double dec_delta,
+ double *out)
+{
+ double corr;
+
+ /* The bottom vertices, note that the positive right ascension is on
+ the left side. */
+ out[1] = out[3] = dec_center - dec_delta/2;
+ corr=1/(2*cos(out[1]*M_PI/180));
+ out[0] = ra_center + ra_delta*corr;
+ out[2] = ra_center - ra_delta*corr;
+
+ /* The top vertices. */
+ out[5] = out[7] = dec_center + dec_delta/2;
+ corr=1/(2*cos(out[5]*M_PI/180));
+ out[4] = ra_center + ra_delta*corr;
+ out[6] = ra_center - ra_delta*corr;
+}
+
+
+
+
+
/* Return the pixel scale of the dataset in units of the WCS. */
double *
gal_wcs_pixel_scale(struct wcsprm *wcs)