[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master de51acf8: Library (wcs.h): convert points from
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master de51acf8: Library (wcs.h): convert points from one coordsys to another |
|
Date: |
Sun, 17 Sep 2023 18:32:32 -0400 (EDT) |
branch: master
commit de51acf8cf14368a6582eac6eb97fa8a09d1f058
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Library (wcs.h): convert points from one coordsys to another
Until now, we had no function to convert a set of coordinates in one
celestial coordinate system to another (for example to convert the RA and
Dec columns in a table to galactic coordinates).
With this commit, the necessary functions have been added within the WCS
library; enabling 30 new operators for Arithmetic (for conversion to and
from all the 6 exisiting coordinate systems).
---
NEWS | 53 ++++++-
doc/gnuastro.texi | 168 ++++++++++++++++++++-
lib/arithmetic.c | 337 ++++++++++++++++++++++++++++++++++++++++++
lib/gnuastro/arithmetic.h | 36 ++++-
lib/gnuastro/wcs.h | 7 +
lib/wcs.c | 362 ++++++++++++++++++++++++++++++++++++++++------
6 files changed, 907 insertions(+), 56 deletions(-)
diff --git a/NEWS b/NEWS
index 3c5222e0..b6b7048c 100644
--- a/NEWS
+++ b/NEWS
@@ -43,17 +43,55 @@ See the end of the file for license conditions.
Arithmetic:
- New operators (see book for a full description):
- pool-min: Min-pooling to reduce the size of the input by calculating
- the minimum of a the pixels within the pooling window (accounting for
- a stride). See the new "Pooling opeators" section of the book for
- more. The pooling operators were all implemented by Faezeh
- Bidjarchian.
+ the minimum of a the pixels within the pooling window
+ (accounting for a stride). See the new "Pooling opeators"
+ section of the book for more. The pooling operators were
+ all implemented by Faezeh Bidjarchian.
- pool-max: Similar to 'pool-min' but using maximum.
- pool-sum: Similar to 'pool-min' but using sum.
- pool-mean: Similar to 'pool-min' but using mean.
- pool-median: Similar to 'pool-min' but using median.
- to-1d: convert the input operand into a 1D array, no matter how many
- dimensions it has. Suggested by Faezeh Bidjarchian.
+ dimensions it has. Suggested by Faezeh Bidjarchian.
- trim: remove all fully-blank outer regions of the input dataset.
+ - eq-b1950-to-eq-j2000: Convert input equatorial coordinates (RA and
+ Dec) in the B1950 equinox to equatorial J2000
+ equinox.
+ - eq-b1950-to-ec-b1950: Same input as above; to ecliptic B1950.
+ - eq-b1950-to-ec-j2000: Same input as above; to ecliptic J2000.
+ - eq-b1950-to-galactic: Same input as above; to Galactic.
+ - eq-b1950-to-supergalactic: Same input as above; to Supergalactic.
+ - eq-j2000-to-eq-b1950: Convert input equatorial coordinates (RA and
+ Dec) in the J2000 equinox to equatorial B1950
+ equinox.
+ - eq-j2000-to-ec-b1950: Same input as above; to ecliptic B1950.
+ - eq-j2000-to-ec-j2000: Same input as above; to ecliptic J2000.
+ - eq-j2000-to-galactic: Same input as above; to Galactic.
+ - eq-j2000-to-supergalactic: Same input as above; to Supergalactic.
+ - ec-b1950-to-eq-b1950: Convert input ecliptic coordinates in the
+ B1950 equinox to equatorial B1950 equinox.
+ - ec-b1950-to-eq-j2000: Same input as above; to equatorial B1950.
+ - ec-b1950-to-ec-j2000: Same input as above; to ecliptic J2000.
+ - ec-b1950-to-galactic: Same input as above; to Galactic.
+ - ec-b1950-to-supergalactic: Same input as above; to Supergalactic
+ - ec-j2000-to-eq-b1950: Convert input ecliptic coordinates in the
+ J2000 equinox to equatorial B1950 equinox.
+ - ec-j2000-to-eq-j2000: Same input as above; to equatorial J2000.
+ - ec-j2000-to-ec-b1950: Same input as above; to ecliptic B1950.
+ - ec-j2000-to-galactic: Same input as above; to Galactic.
+ - ec-j2000-to-supergalactic: Same input as above; to Supergalactic.
+ - galactic-to-eq-b1950: Convert input Galactic coordinates to
+ equatorial B1950 equinox.
+ - galactic-to-eq-j2000: Same input as above; to equatorial J2000.
+ - galactic-to-ec-b1950: Same input as above; to ecliptic B1950.
+ - galactic-to-ec-j2000: Same input as above; to ecliptic J2000.
+ - galactic-to-supergalactic: Same input as above; to Supergalactic.
+ - supergalactic-to-eq-b1950: Convert input Supergalactic coordinates to
+ equatorial B1950 equinox.
+ - supergalactic-to-eq-j2000: Same input as above; to equatorial J2000.
+ - supergalactic-to-ec-b1950: Same input as above; to ecliptic B1950.
+ - supergalactic-to-ec-j2000: Same input as above; to ecliptic J2000.
+ - supergalactic-to-galactic: Same input as above; to Galactic.
ConvertType:
- It is now possible to write TIFF files in the output (until now,
@@ -103,6 +141,11 @@ See the end of the file for license conditions.
Fathma Mehnoor.
- gal_wcs_projection_name_to_id: convert projection names to IDs.
- gal_wcs_projection_name_from_id: convert projection IDs to their name.
+ - gal_wcs_coordsys_convert_points: convert the input set of points from
+ one celestial coordinate system to another.
+ - gal_wcs_coordsys_sys1_ref_in_sys2: return the longitude reference point
+ of the first celestial coordinate system in the second coordinate
+ system's longitude and latitude.
** Removed features
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 60903615..9afe4844 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -550,6 +550,7 @@ Arithmetic operators
* Basic mathematical operators:: For example, +, -, /, log, and pow.
* Trigonometric and hyperbolic operators:: sin, cos, atan, asinh, etc.
* Constants:: Physical and Mathematical constants.
+* Coordinate conversion operators:: For example equatorial J2000 to Galactic.
* Unit conversion operators:: Various unit conversions necessary.
* Statistical operators:: Statistics of a single dataset (for example,
mean).
* Stacking operators:: Coadding or combining multiple datasets into
one.
@@ -562,7 +563,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.
-* Box shape operators::
+* Box shape operators:: Return edges of 2D boxes.
* 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.
@@ -15266,6 +15267,7 @@ This is a very useful option for operations on the FITS
date values, for example
The advantage of working with the Unix epoch time is that you do not have to
worry about calendar details (for example, the number of days in different
months, or leap years).
@item --wcscoordsys=STR
+@cindex FK5
@cindex Galactic coordinate system
@cindex Ecliptic coordinate system
@cindex Equatorial coordinate system
@@ -15290,6 +15292,17 @@ The currently recognized coordinate systems are listed
below (the most common on
@table @code
@item eq-j2000
2000.0 (Julian-year) equatorial coordinates.
+This is also known as FK5 (short for ``Fundamental Katalog No 5'' which was
the source of the star coordinates used to define it).
+
+This coordinate system is based on the motion of the Sun and has epochs when
the mean equator was used (for example @code{eq-b1950} below).
+Furthermore, the definition of year is different: either the Besselian year in
1950.0, or the Julian year in 2000.
+For more on their difference and links for further reading about epochs in
astronomy, please see the description in
@url{https://en.wikipedia.org/wiki/Epoch_(astronomy), Wikipedia}.
+
+Because of these difficulties, the equatorial J2000.0 coordinate system has
been deprecated by the IAU in favor of International Celestial Refernece System
(ICRS) but is still used extensively.
+ICRS is defined based on extra-galactic quasars, so it does not depend on the
dynamics of the solar system any more.
+But to enable historical continuity, ICRS has been defined to be equivalent to
the equatorial J2000.0 within its accepted error bars of the latter (tens of
milli-arcseconds).
+This justifies the reason that moving to ICRS has been relatively slow.
+
@item eq-b1950
1950.0 (Besselian-year) equatorial coordinates.
@item ec-j2000
@@ -15302,9 +15315,6 @@ Galactic coordinates.
Supergalactic coordinates.
@end table
-The Equatorial and Ecliptic coordinate systems are defined by the mean equator
and equinox epoch: either the Besselian year 1950.0, or the Julian year 2000.
-For more on their difference and links for further reading about epochs in
astronomy, please see the description in
@url{https://en.wikipedia.org/wiki/Epoch_(astronomy), Wikipedia}.
-
@item --wcsdistortion=STR
@cindex WCS distortion
@cindex Distortion, WCS
@@ -19883,6 +19893,7 @@ Reading NaN as a floating point number in Gnuastro is
not case-sensitive.
* Basic mathematical operators:: For example, +, -, /, log, and pow.
* Trigonometric and hyperbolic operators:: sin, cos, atan, asinh, etc.
* Constants:: Physical and Mathematical constants.
+* Coordinate conversion operators:: For example equatorial J2000 to Galactic.
* Unit conversion operators:: Various unit conversions necessary.
* Statistical operators:: Statistics of a single dataset (for example,
mean).
* Stacking operators:: Coadding or combining multiple datasets into
one.
@@ -19895,7 +19906,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.
-* Box shape operators::
+* Box shape operators:: Return edges of 2D boxes.
* 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.
@@ -20102,7 +20113,7 @@ Inverse Hyperbolic sine, cosine, and tangent.
These operators take a single operand.
@end table
-@node Constants, Unit conversion operators, Trigonometric and hyperbolic
operators, Arithmetic operators
+@node Constants, Coordinate conversion operators, Trigonometric and hyperbolic
operators, Arithmetic operators
@subsubsection Constants
@cindex Pi
During your analysis it is often necessary to have certain constants like the
number @mymath{\pi}.
@@ -20162,7 +20173,104 @@ The fine-structure constant (no units).
See @url{https://en.wikipedia.org/wiki/Fine-structure_constant, Wikipedia}.
@end table
-@node Unit conversion operators, Statistical operators, Constants, Arithmetic
operators
+@node Coordinate conversion operators, Unit conversion operators, Constants,
Arithmetic operators
+@subsubsection Coordinate conversion operators
+
+@cindex Galactic coordinate system
+@cindex Ecliptic coordinate system
+@cindex Equatorial coordinate system
+Different celestial coordinate systems are useful for different scenarios.
+For example, assume you have the RA and Dec of large sample of galaxies that
you plan to study the halos of galaxies from.
+For such studies, you prefer to stay as far away as possible from the Galactic
plane, because the density of stars and interstellar filaments (cirrus)
significantly increases as you get close to the Milky way's disk.
+But the @url{https://en.wikipedia.org/wiki/Equatorial_coordinate_system,
Equatorial coordinate system} which defines the RA and Dec and is based on
Earth's equator; and does not show the position of your objects in relationt to
the galactic disk.
+
+The best way forward in the example above is to convert your RA and Dec table
into the @url{https://en.wikipedia.org/wiki/Galactic_coordinate_system,
Galactic coordinate system}; and select those with a large (positive or
negative) Galactic latitude.
+Alternatively, if you observe a bright point on a galaxy and want to confirm
if it was actually a super-nova and not a moving asteriod, a first step is to
convert your RA and Dec to the
@url{https://en.wikipedia.org/wiki/Ecliptic_coordinate_system, Ecliptic
coordinate system} and confirm if you are sufficiently distant from the
ecliptic (plane of the Solar System; where fast moving objects are most common).
+
+The operators described in this section are precisely for the purpose above:
to convert various celestial coordinate systems that are supported within
Gnuastro into each other.
+For example, if you want to convert the RA and Dec equatorial (at the Julian
year 2000 equinox) coordinates (within the @code{RA} and @code{DEC} columns) of
@file{points.fits} into Galactic longitude and latitutde, you can use the
command below (the column metadata are not mandatory, but to avoid later
confusion, it is always good to have them in your output.
+
+@example
+$ asttable points.fits -c'arith RA DEC eq-j2000-to-galactic' \
+ --colmetadata=1,GLON,deg,"Galactic longitude" \
+ --colmetadata=2,GLAT,deg,"Galactic latitude" \
+ --output=points-gal.fits
+@end example
+
+One important thing to consider is that the equatorial and ecliptic
coordinates are not static: they include the dynamics of Earth in the solar
system: in particular, the reference point on the equator moves over decades.
+Therefore these two (equatorial and ecliptic) coordinate systems are defined
within epochs: the 1950 epoch is defined by
@url{https://en.wikipedia.org/wiki/Epoch_(astronomy)#Besselian_years, Besselian
years}, while the 2000 epoch is defined in
@url{https://en.wikipedia.org/wiki/Epoch_(astronomy)#Julian_years_and_J2000,
Julian years}.
+So when dealing with these coordinates one of the `@code{-b1950}' or
`@code{-j2000}' suffixes are necessary (for example @code{eq-j2000} or
@code{ec-b1950}).
+
+@cindex ICRS
+The Galactic or Supergalactic coordiantes are not defined based on the Earth's
dynamics; therefore they do not have any epoch associated with them.
+Extra-galactic studies do not depend on the dynamics of the earth, but the
equatorial coordinate system is the most dominant in that field.
+Therefore in its 23rd General Assembly, the International Astronomical Union
approved the
@url{https://en.wikipedia.org/wiki/International_Celestial_Reference_System_and_its_realizations,
International Celestial Reference System} or ICRS based on quasars (which are
static within our observational limitations)viewed through long baseline radio
interferometry (the most accurate method of observation that we currently have).
+ICRS is designed to be within the errors of the Equatorial J2000 coordinate
system, so they are currently very similar; but ICRS has much better accuracy.
+We will be adding ICRS in the operators below soon.
+
+@strong{Floating point errors:} The operation to convert between the
coordinate systems involves many sines, cosines (and their inverse).
+Therefore, floating point errors (due to the limited precision of the
definition of floating points in bits) can cause small offests.
+For example see the code below were we convert equatorial to galactic and
back, then compare the input and output (which is in the 5th and 6th decimal of
a degree; or about 0.2 or 0.01 arcseconds).
+
+@example
+$ sys1=eq-j2000
+$ sys2=galactic
+$ echo "10.2345689 45.6789012" \
+ | asttable -Afixed -B8 \
+ -c'arith $1 $2 '$sys1'-to-'$sys2' \
+ '$sys2'-to-'$sys1' set-lat set-lng \
+ lng $1 - lat $2 -'
+0.00000363 -0.00007725
+@end example
+
+If you set @code{sys2=ec-j2000} or @code{sys2=supergalactic}, it will be zero
until the full set of 8 decimals that are printed here (the displayed precision
can be changed with the value of the @code{-B} option above).
+It is therefore useful to have your original coordinates (in the same table
for example) and not do too many conversions on conversions (to propagate this
problem).
+
+@table @code
+@item eq-b1950-to-eq-j2000
+@itemx eq-b1950-to-ec-b1950
+@itemx eq-b1950-to-ec-j2000
+@itemx eq-b1950-to-galactic
+@itemx eq-b1950-to-supergalactic
+Convert Equatorial (B1950 equinox) coordinates into the respective coordinate
system within each operator's name.
+
+@item eq-j2000-to-eq-b1950
+@itemx eq-j2000-to-ec-b1950
+@itemx eq-j2000-to-ec-j2000
+@itemx eq-j2000-to-galactic
+@itemx eq-j2000-to-supergalactic
+Convert Equatorial (J2000 equinox) coordinates into the respective coordinate
system within each operator's name.
+
+@item ec-b1950-to-eq-b1950
+@itemx ec-b1950-to-eq-j2000
+@itemx ec-b1950-to-ec-j2000
+@itemx ec-b1950-to-galactic
+@itemx ec-b1950-to-supergalactic
+Convert Ecliptic (B1950 equinox) coordinates into the respective coordinate
system within each operator's name.
+
+@item ec-j2000-to-eq-b1950
+@itemx ec-j2000-to-eq-j2000
+@itemx ec-j2000-to-ec-b1950
+@itemx ec-j2000-to-galactic
+@itemx ec-j2000-to-supergalactic
+Convert Ecliptic (J2000 equinox) coordinates into the respective coordinate
system within each operator's name.
+
+@item galactic-to-eq-b1950
+@itemx galactic-to-eq-j2000
+@itemx galactic-to-ec-b1950
+@itemx galactic-to-ec-j2000
+@itemx galactic-to-supergalactic
+Convert Galactic coordinates into the respective coordinate system within each
operator's name.
+
+@item supergalactic-to-eq-b1950
+@itemx supergalactic-to-eq-j2000
+@itemx supergalactic-to-ec-b1950
+@itemx supergalactic-to-ec-j2000
+@itemx supergalactic-to-galactic
+Convert Supergalactic coordinates into the respective coordinate system within
each operator's name.
+@end table
+
+@node Unit conversion operators, Statistical operators, Coordinate conversion
operators, Arithmetic operators
@subsubsection Unit conversion operators
It often happens that you have data in one unit (for example, counts on your
CCD), but would like to convert it into another (for example, magnitudes, to
measure the brightness of a galaxy).
@@ -38566,6 +38674,18 @@ The Gnuastro WCS distortion identifiers are defined in
the @code{GAL_WCS_COORDSY
Since the returned dataset is newly allocated, if you do not need the original
dataset after this, use the WCSLIB library function @code{wcsfree} to free the
input, for example, @code{wcsfree(inwcs)}.
@end deftypefun
+@deftypefun void gal_wcs_coordsys_convert_points (int @code{sys1}, double
@code{*lng1_d}, double @code{*lat1_d}, int @code{sys2}, double @code{*lng2_d},
double @code{*lat2_d}, size_t @code{number})
+Convert the input set of longitudes (@code{lng1_d}, in degrees) and latitudes
(@code{lat1_d}, in degrees) within a recognized coordinate system (@code{sys1};
one of the @code{GAL_WCS_COORDSYS_*} macros above) into an output coordinate
system (@code{sys2}).
+The output values are written in @code{lng2_d} and @code{lng2_d}.
+The total number of points should be given in @code{number}.
+If you want the operation to be done in place (without allocating a new
dataset), give the same pointers to the coordinate arguments.
+@end deftypefun
+
+@deftypefun void gal_wcs_coordsys_sys1_ref_in_sys2 (int @code{sys1}, int
@code{sys2}, double @code{*lng2}, double @code{*lat2})
+Return the longitude and latitude of the reference point (on the equator) of
the first coordinate system (@code{sys1}) within the second system
(@code{sys2}).
+Coordinate systems are identifed by the @code{GAL_WCS_COORDSYS_*} macros above.
+@end deftypefun
+
@cindex WCS distortion
@cindex Distortion, WCS
@deftypefun {int} gal_wcs_distortion_identify (struct wcsprm @code{*wcs})
@@ -39067,6 +39187,40 @@ Return the first dataset, but with the second dataset
being placed in the @code{
This is useful to swap the operators on the stacks of the higher-level
programs that call the arithmetic library.
@end deffn
+@deffn Macro GAL_ARITHMETIC_OP_EQB1950_TO_EQJ2000
+@deffnx Macro GAL_ARITHMETIC_OP_EQB1950_TO_ECB1950
+@deffnx Macro GAL_ARITHMETIC_OP_EQB1950_TO_ECJ2000
+@deffnx Macro GAL_ARITHMETIC_OP_EQB1950_TO_GALACTIC
+@deffnx Macro GAL_ARITHMETIC_OP_EQB1950_TO_SUPERGALACTIC
+@deffnx Macro GAL_ARITHMETIC_OP_EQJ2000_TO_EQB1950
+@deffnx Macro GAL_ARITHMETIC_OP_EQJ2000_TO_ECB1950
+@deffnx Macro GAL_ARITHMETIC_OP_EQJ2000_TO_ECJ2000
+@deffnx Macro GAL_ARITHMETIC_OP_EQJ2000_TO_GALACTIC
+@deffnx Macro GAL_ARITHMETIC_OP_EQJ2000_TO_SUPERGALACTIC
+@deffnx Macro GAL_ARITHMETIC_OP_ECB1950_TO_EQB1950
+@deffnx Macro GAL_ARITHMETIC_OP_ECB1950_TO_EQJ2000
+@deffnx Macro GAL_ARITHMETIC_OP_ECB1950_TO_ECJ2000
+@deffnx Macro GAL_ARITHMETIC_OP_ECB1950_TO_GALACTIC
+@deffnx Macro GAL_ARITHMETIC_OP_ECB1950_TO_SUPERGALACTIC
+@deffnx Macro GAL_ARITHMETIC_OP_ECJ2000_TO_EQB1950
+@deffnx Macro GAL_ARITHMETIC_OP_ECJ2000_TO_EQJ2000
+@deffnx Macro GAL_ARITHMETIC_OP_ECJ2000_TO_ECB1950
+@deffnx Macro GAL_ARITHMETIC_OP_ECJ2000_TO_GALACTIC
+@deffnx Macro GAL_ARITHMETIC_OP_ECJ2000_TO_SUPERGALACTIC
+@deffnx Macro GAL_ARITHMETIC_OP_GALACTIC_TO_EQB1950
+@deffnx Macro GAL_ARITHMETIC_OP_GALACTIC_TO_EQJ2000
+@deffnx Macro GAL_ARITHMETIC_OP_GALACTIC_TO_ECB1950
+@deffnx Macro GAL_ARITHMETIC_OP_GALACTIC_TO_ECJ2000
+@deffnx Macro GAL_ARITHMETIC_OP_GALACTIC_TO_SUPERGALACTIC
+@deffnx Macro GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_EQB1950
+@deffnx Macro GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_EQJ2000
+@deffnx Macro GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_ECB1950
+@deffnx Macro GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_ECJ2000
+@deffnx Macro GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_GALACTIC
+Operators that convert recognized celestial coordinates to and from each other.
+They all take two operands and return two @code{gal_data_t}s (as a list).
+For more on celestial coordinate conversion, see @ref{Coordinate conversion
operators}.
+@end deffn
@deftypefun {gal_data_t *} gal_arithmetic (int @code{operator}, size_t
@code{numthreads}, int @code{flags}, ...)
Apply the requested arithmetic operator on the operand(s).
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 4f0f6438..193a6b9a 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -678,6 +678,183 @@ arithmetic_function_unary(int operator, int flags,
gal_data_t *in)
+/* These are unit-conversion operators that need two components (like
+ 2D coordinate system conversion). */
+static gal_data_t *
+arithmetic_unit_binary(int operator, int flags, gal_data_t *a_in,
+ gal_data_t *b_in)
+{
+ int sys1, sys2;
+ gal_data_t *a, *b;
+
+ /* If the input datasets are empty, abort. */
+ if(a_in==NULL || b_in==NULL)
+ error(EXIT_FAILURE, 0, "%s: at least one of the input "
+ "datasets ('a_in' and 'b_in') is NULL", __func__);
+
+ /* Make sure the two inputs have the same size. */
+ if(a_in->size != b_in->size)
+ error(EXIT_FAILURE, 0, "%s: the two inputs have different "
+ "sizes: %zu for 'a_in' and %zu for 'b_in'", __func__,
+ a_in->size, b_in->size);
+
+ /* Set any potentially 'next' element of the inputs to NULL. */
+ if(a_in->next && a_in->next!=b_in) gal_data_free(a_in->next);
+ if(b_in->next && b_in->next!=a_in) gal_data_free(b_in->next);
+ a_in->next = b_in->next = NULL;
+
+ /* Convert the inputs to double-precision floating points.*/
+ a = ( flags & GAL_ARITHMETIC_FLAG_FREE
+ ? gal_data_copy_to_new_type_free(a_in, GAL_TYPE_FLOAT64)
+ : gal_data_copy_to_new_type(a_in, GAL_TYPE_FLOAT64) );
+ b = ( flags & GAL_ARITHMETIC_FLAG_FREE
+ ? gal_data_copy_to_new_type_free(b_in, GAL_TYPE_FLOAT64)
+ : gal_data_copy_to_new_type(b_in, GAL_TYPE_FLOAT64) );
+
+ /* Set the two systems. */
+ switch(operator)
+ {
+ case GAL_ARITHMETIC_OP_EQB1950_TO_EQJ2000:
+ sys1=GAL_WCS_COORDSYS_EQB1950;
+ sys2=GAL_WCS_COORDSYS_EQJ2000;
+ break;
+ case GAL_ARITHMETIC_OP_EQB1950_TO_ECB1950:
+ sys1=GAL_WCS_COORDSYS_EQB1950;
+ sys2=GAL_WCS_COORDSYS_ECB1950;
+ break;
+ case GAL_ARITHMETIC_OP_EQB1950_TO_ECJ2000:
+ sys1=GAL_WCS_COORDSYS_EQB1950;
+ sys2=GAL_WCS_COORDSYS_ECJ2000;
+ break;
+ case GAL_ARITHMETIC_OP_EQB1950_TO_GALACTIC:
+ sys1=GAL_WCS_COORDSYS_EQB1950;
+ sys2=GAL_WCS_COORDSYS_GALACTIC;
+ break;
+ case GAL_ARITHMETIC_OP_EQB1950_TO_SUPERGALACTIC:
+ sys1=GAL_WCS_COORDSYS_EQB1950;
+ sys2=GAL_WCS_COORDSYS_SUPERGALACTIC;
+ break;
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_EQB1950:
+ sys1=GAL_WCS_COORDSYS_EQJ2000;
+ sys2=GAL_WCS_COORDSYS_EQB1950;
+ break;
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_ECB1950:
+ sys1=GAL_WCS_COORDSYS_EQJ2000;
+ sys2=GAL_WCS_COORDSYS_ECB1950;
+ break;
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_ECJ2000:
+ sys1=GAL_WCS_COORDSYS_EQJ2000;
+ sys2=GAL_WCS_COORDSYS_ECJ2000;
+ break;
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_GALACTIC:
+ sys1=GAL_WCS_COORDSYS_EQJ2000;
+ sys2=GAL_WCS_COORDSYS_GALACTIC;
+ break;
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_SUPERGALACTIC:
+ sys1=GAL_WCS_COORDSYS_EQJ2000;
+ sys2=GAL_WCS_COORDSYS_SUPERGALACTIC;
+ break;
+ case GAL_ARITHMETIC_OP_ECB1950_TO_EQB1950:
+ sys1=GAL_WCS_COORDSYS_ECB1950;
+ sys2=GAL_WCS_COORDSYS_EQB1950;
+ break;
+ case GAL_ARITHMETIC_OP_ECB1950_TO_EQJ2000:
+ sys1=GAL_WCS_COORDSYS_ECB1950;
+ sys2=GAL_WCS_COORDSYS_EQJ2000;
+ break;
+ case GAL_ARITHMETIC_OP_ECB1950_TO_ECJ2000:
+ sys1=GAL_WCS_COORDSYS_ECB1950;
+ sys2=GAL_WCS_COORDSYS_ECJ2000;
+ break;
+ case GAL_ARITHMETIC_OP_ECB1950_TO_GALACTIC:
+ sys1=GAL_WCS_COORDSYS_ECB1950;
+ sys2=GAL_WCS_COORDSYS_GALACTIC;
+ break;
+ case GAL_ARITHMETIC_OP_ECB1950_TO_SUPERGALACTIC:
+ sys1=GAL_WCS_COORDSYS_ECB1950;
+ sys2=GAL_WCS_COORDSYS_SUPERGALACTIC;
+ break;
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_EQB1950:
+ sys1=GAL_WCS_COORDSYS_ECJ2000;
+ sys2=GAL_WCS_COORDSYS_EQB1950;
+ break;
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_EQJ2000:
+ sys1=GAL_WCS_COORDSYS_ECJ2000;
+ sys2=GAL_WCS_COORDSYS_EQJ2000;
+ break;
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_ECB1950:
+ sys1=GAL_WCS_COORDSYS_ECJ2000;
+ sys2=GAL_WCS_COORDSYS_ECB1950;
+ break;
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_GALACTIC:
+ sys1=GAL_WCS_COORDSYS_ECJ2000;
+ sys2=GAL_WCS_COORDSYS_GALACTIC;
+ break;
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_SUPERGALACTIC:
+ sys1=GAL_WCS_COORDSYS_ECJ2000;
+ sys2=GAL_WCS_COORDSYS_SUPERGALACTIC;
+ break;
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_EQB1950:
+ sys1=GAL_WCS_COORDSYS_GALACTIC;
+ sys2=GAL_WCS_COORDSYS_EQB1950;
+ break;
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_EQJ2000:
+ sys1=GAL_WCS_COORDSYS_GALACTIC;
+ sys2=GAL_WCS_COORDSYS_EQJ2000;
+ break;
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_ECB1950:
+ sys1=GAL_WCS_COORDSYS_GALACTIC;
+ sys2=GAL_WCS_COORDSYS_ECB1950;
+ break;
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_ECJ2000:
+ sys1=GAL_WCS_COORDSYS_GALACTIC;
+ sys2=GAL_WCS_COORDSYS_ECJ2000;
+ break;
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_SUPERGALACTIC:
+ sys1=GAL_WCS_COORDSYS_GALACTIC;
+ sys2=GAL_WCS_COORDSYS_SUPERGALACTIC;
+ break;
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_EQB1950:
+ sys1=GAL_WCS_COORDSYS_GALACTIC;
+ sys2=GAL_WCS_COORDSYS_EQB1950;
+ break;
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_EQJ2000:
+ sys1=GAL_WCS_COORDSYS_SUPERGALACTIC;
+ sys2=GAL_WCS_COORDSYS_EQJ2000;
+ break;
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_ECB1950:
+ sys1=GAL_WCS_COORDSYS_SUPERGALACTIC;
+ sys2=GAL_WCS_COORDSYS_ECB1950;
+ break;
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_ECJ2000:
+ sys1=GAL_WCS_COORDSYS_SUPERGALACTIC;
+ sys2=GAL_WCS_COORDSYS_ECJ2000;
+ break;
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_GALACTIC:
+ sys1=GAL_WCS_COORDSYS_SUPERGALACTIC;
+ sys2=GAL_WCS_COORDSYS_GALACTIC;
+ 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 "
+ "in this function", __func__, PACKAGE_BUGREPORT, operator);
+ }
+
+ /* Call the main conversion function (and write the values in place). */
+ gal_wcs_coordsys_convert_points(sys1, a->array, b->array,
+ sys2, a->array, b->array, a->size);
+
+ /* Clean up and return. */
+ if( (flags & GAL_ARITHMETIC_FLAG_FREE) && a!=a_in)
+ { gal_data_free(a_in); gal_data_free(b_in); }
+ b->next=a;
+ return b;
+}
+
+
+
+
+
/* Call functions in the 'gnuastro/statistics' library. */
static gal_data_t *
arithmetic_from_statistics(int operator, int flags, gal_data_t *input)
@@ -3316,6 +3493,68 @@ gal_arithmetic_set_operator(char *string, size_t
*num_operands)
else if( !strcmp(string, "au-to-ly"))
{ op=GAL_ARITHMETIC_OP_AU_TO_LY; *num_operands=1; }
+ /* Celestial coordinate conversions. */
+ else if (!strcmp(string, "eq-b1950-to-eq-j2000"))
+ { op=GAL_ARITHMETIC_OP_EQB1950_TO_EQJ2000; *num_operands=2; }
+ else if (!strcmp(string, "eq-b1950-to-ec-b1950"))
+ { op=GAL_ARITHMETIC_OP_EQB1950_TO_ECB1950; *num_operands=2; }
+ else if (!strcmp(string, "eq-b1950-to-ec-j2000"))
+ { op=GAL_ARITHMETIC_OP_EQB1950_TO_ECJ2000; *num_operands=2; }
+ else if (!strcmp(string, "eq-b1950-to-galactic"))
+ { op=GAL_ARITHMETIC_OP_EQB1950_TO_GALACTIC; *num_operands=2; }
+ else if (!strcmp(string, "eq-b1950-to-supergalactic"))
+ { op=GAL_ARITHMETIC_OP_EQB1950_TO_SUPERGALACTIC; *num_operands=2; }
+ else if (!strcmp(string, "eq-j2000-to-eq-b1950"))
+ { op=GAL_ARITHMETIC_OP_EQJ2000_TO_EQB1950; *num_operands=2; }
+ else if (!strcmp(string, "eq-j2000-to-ec-b1950"))
+ { op=GAL_ARITHMETIC_OP_EQJ2000_TO_ECB1950; *num_operands=2; }
+ else if (!strcmp(string, "eq-j2000-to-ec-j2000"))
+ { op=GAL_ARITHMETIC_OP_EQJ2000_TO_ECJ2000; *num_operands=2; }
+ else if (!strcmp(string, "eq-j2000-to-galactic"))
+ { op=GAL_ARITHMETIC_OP_EQJ2000_TO_GALACTIC; *num_operands=2; }
+ else if (!strcmp(string, "eq-j2000-to-supergalactic"))
+ { op=GAL_ARITHMETIC_OP_EQJ2000_TO_SUPERGALACTIC; *num_operands=2; }
+ else if (!strcmp(string, "ec-b1950-to-eq-b1950"))
+ { op=GAL_ARITHMETIC_OP_ECB1950_TO_EQB1950; *num_operands=2; }
+ else if (!strcmp(string, "ec-b1950-to-eq-j2000"))
+ { op=GAL_ARITHMETIC_OP_ECB1950_TO_EQJ2000; *num_operands=2; }
+ else if (!strcmp(string, "ec-b1950-to-ec-j2000"))
+ { op=GAL_ARITHMETIC_OP_ECB1950_TO_ECJ2000; *num_operands=2; }
+ else if (!strcmp(string, "ec-b1950-to-galactic"))
+ { op=GAL_ARITHMETIC_OP_ECB1950_TO_GALACTIC; *num_operands=2; }
+ else if (!strcmp(string, "ec-b1950-to-supergalactic"))
+ { op=GAL_ARITHMETIC_OP_ECB1950_TO_SUPERGALACTIC; *num_operands=2; }
+ else if (!strcmp(string, "ec-j2000-to-eq-b1950"))
+ { op=GAL_ARITHMETIC_OP_ECJ2000_TO_EQB1950; *num_operands=2; }
+ else if (!strcmp(string, "ec-j2000-to-eq-j2000"))
+ { op=GAL_ARITHMETIC_OP_ECJ2000_TO_EQJ2000; *num_operands=2; }
+ else if (!strcmp(string, "ec-j2000-to-ec-b1950"))
+ { op=GAL_ARITHMETIC_OP_ECJ2000_TO_ECB1950; *num_operands=2; }
+ else if (!strcmp(string, "ec-j2000-to-galactic"))
+ { op=GAL_ARITHMETIC_OP_ECJ2000_TO_GALACTIC; *num_operands=2; }
+ else if (!strcmp(string, "ec-j2000-to-supergalactic"))
+ { op=GAL_ARITHMETIC_OP_ECJ2000_TO_SUPERGALACTIC; *num_operands=2; }
+ else if (!strcmp(string, "galactic-to-eq-b1950"))
+ { op=GAL_ARITHMETIC_OP_GALACTIC_TO_EQB1950; *num_operands=2; }
+ else if (!strcmp(string, "galactic-to-eq-j2000"))
+ { op=GAL_ARITHMETIC_OP_GALACTIC_TO_EQJ2000; *num_operands=2; }
+ else if (!strcmp(string, "galactic-to-ec-b1950"))
+ { op=GAL_ARITHMETIC_OP_GALACTIC_TO_ECB1950; *num_operands=2; }
+ else if (!strcmp(string, "galactic-to-ec-j2000"))
+ { op=GAL_ARITHMETIC_OP_GALACTIC_TO_ECJ2000; *num_operands=2; }
+ else if (!strcmp(string, "galactic-to-supergalactic"))
+ { op=GAL_ARITHMETIC_OP_GALACTIC_TO_SUPERGALACTIC; *num_operands=2; }
+ else if (!strcmp(string, "supergalactic-to-eq-b1950"))
+ { op=GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_EQB1950; *num_operands=2; }
+ else if (!strcmp(string, "supergalactic-to-eq-j2000"))
+ { op=GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_EQJ2000; *num_operands=2; }
+ else if (!strcmp(string, "supergalactic-to-ec-b1950"))
+ { op=GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_ECB1950; *num_operands=2; }
+ else if (!strcmp(string, "supergalactic-to-ec-j2000"))
+ { op=GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_ECJ2000; *num_operands=2; }
+ else if (!strcmp(string, "supergalactic-to-galactic"))
+ { op=GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_GALACTIC; *num_operands=2; }
+
/* Statistical/higher-level operators. */
else if (!strcmp(string, "minvalue"))
{ op=GAL_ARITHMETIC_OP_MINVAL; *num_operands=1; }
@@ -3658,6 +3897,68 @@ gal_arithmetic_operator_string(int operator)
case GAL_ARITHMETIC_OP_POOLMEAN: return "pool-mean";
case GAL_ARITHMETIC_OP_POOLMEDIAN: return "pool-median";
+ case GAL_ARITHMETIC_OP_EQB1950_TO_EQJ2000:
+ return "eq-b1950-to-eq-j2000";
+ case GAL_ARITHMETIC_OP_EQB1950_TO_ECB1950:
+ return "eq-b1950-to-eq-b1950";
+ case GAL_ARITHMETIC_OP_EQB1950_TO_ECJ2000:
+ return "eq-b1950-to-ec-j2000";
+ case GAL_ARITHMETIC_OP_EQB1950_TO_GALACTIC:
+ return "eq-b1950-to-galactic";
+ case GAL_ARITHMETIC_OP_EQB1950_TO_SUPERGALACTIC:
+ return "eq-b1950-to-supergalactic";
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_EQB1950:
+ return "eq-j2000-to-eq-b1950";
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_ECB1950:
+ return "eq-j2000-to-ec-b1950";
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_ECJ2000:
+ return "eq-j2000-to-eq-j2000";
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_GALACTIC:
+ return "eq-j2000-to-galactic";
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_SUPERGALACTIC:
+ return "eq-j2000-to-supergalactic";
+ case GAL_ARITHMETIC_OP_ECB1950_TO_EQB1950:
+ return "ec-b1950-to-eq-b1950";
+ case GAL_ARITHMETIC_OP_ECB1950_TO_EQJ2000:
+ return "ec-b1950-to-eq-j2000";
+ case GAL_ARITHMETIC_OP_ECB1950_TO_ECJ2000:
+ return "ec-b1950-to-ec-j2000";
+ case GAL_ARITHMETIC_OP_ECB1950_TO_GALACTIC:
+ return "ec-b1950-to-galactic";
+ case GAL_ARITHMETIC_OP_ECB1950_TO_SUPERGALACTIC:
+ return "ec-b1950-to-supergalacitc";
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_EQB1950:
+ return "ec-j2000-to-eq-b1950";
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_EQJ2000:
+ return "ec-j2000-to-eq-j2000";
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_ECB1950:
+ return "ec-j2000-to-ec-b1950";
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_GALACTIC:
+ return "ec-j2000-to-galactic";
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_SUPERGALACTIC:
+ return "ec-j2000-to-supergalactic";
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_EQB1950:
+ return "galactic-to-eq-b1950";
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_EQJ2000:
+ return "galactic-to-eq-j2000";
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_ECB1950:
+ return "galactic-to-ec-b1950";
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_ECJ2000:
+ return "galactic-to-ec-j2000";
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_SUPERGALACTIC:
+ return "galactic-to-supergalactic";
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_EQB1950:
+ return "supergalactic-to-eq-b1950";
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_EQJ2000:
+ return "supergalactic-to-eq-j2000";
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_ECB1950:
+ return "supergalactic-to-ec-b1950";
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_ECJ2000:
+ return "supergalactic-to-ec-j2000";
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_GALACTIC:
+ return "supergalactic-to-galactic";
+
+ /* Un-recognized! */
default: return NULL;
}
return NULL;
@@ -3751,6 +4052,42 @@ gal_arithmetic(int operator, size_t numthreads, int
flags, ...)
out=arithmetic_function_unary(operator, flags, d1);
break;
+ /* 2-component unit conversion. */
+ case GAL_ARITHMETIC_OP_EQB1950_TO_EQJ2000:
+ case GAL_ARITHMETIC_OP_EQB1950_TO_ECB1950:
+ case GAL_ARITHMETIC_OP_EQB1950_TO_ECJ2000:
+ case GAL_ARITHMETIC_OP_EQB1950_TO_GALACTIC:
+ case GAL_ARITHMETIC_OP_EQB1950_TO_SUPERGALACTIC:
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_EQB1950:
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_ECB1950:
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_ECJ2000:
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_GALACTIC:
+ case GAL_ARITHMETIC_OP_EQJ2000_TO_SUPERGALACTIC:
+ case GAL_ARITHMETIC_OP_ECB1950_TO_EQB1950:
+ case GAL_ARITHMETIC_OP_ECB1950_TO_EQJ2000:
+ case GAL_ARITHMETIC_OP_ECB1950_TO_ECJ2000:
+ case GAL_ARITHMETIC_OP_ECB1950_TO_GALACTIC:
+ case GAL_ARITHMETIC_OP_ECB1950_TO_SUPERGALACTIC:
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_EQB1950:
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_EQJ2000:
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_ECB1950:
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_GALACTIC:
+ case GAL_ARITHMETIC_OP_ECJ2000_TO_SUPERGALACTIC:
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_EQB1950:
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_EQJ2000:
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_ECB1950:
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_ECJ2000:
+ case GAL_ARITHMETIC_OP_GALACTIC_TO_SUPERGALACTIC:
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_EQB1950:
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_EQJ2000:
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_ECB1950:
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_ECJ2000:
+ case GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_GALACTIC:
+ d1 = va_arg(va, gal_data_t *);
+ d2 = va_arg(va, gal_data_t *);
+ out=arithmetic_unit_binary(operator, flags, d1, d2);
+ break;
+
/* Binary function operators. */
case GAL_ARITHMETIC_OP_POW:
case GAL_ARITHMETIC_OP_ATAN2:
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index 1c76b4cb..31881e8b 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -153,8 +153,8 @@ enum gal_arithmetic_operators
GAL_ARITHMETIC_OP_SB_TO_MAG, /* Surface Brightness to Magnitude. */
GAL_ARITHMETIC_OP_COUNTS_TO_SB, /* Counts to Surface Brightness. */
GAL_ARITHMETIC_OP_SB_TO_COUNTS, /* Surface Brightness to counts. */
- GAL_ARITHMETIC_OP_COUNTS_TO_JY, /* Counts to Janskys with AB-mag zeropoint.
*/
- GAL_ARITHMETIC_OP_JY_TO_COUNTS, /* Janskys to Counts with AB-mag zeropoint.
*/
+ GAL_ARITHMETIC_OP_COUNTS_TO_JY, /* Counts to Janskys (AB zeropoint). */
+ GAL_ARITHMETIC_OP_JY_TO_COUNTS, /* Janskys to Counts (AB zeropoint). */
GAL_ARITHMETIC_OP_MAG_TO_JY, /* Magnitude to Janskys. */
GAL_ARITHMETIC_OP_JY_TO_MAG, /* Janskys to Magnitude. */
GAL_ARITHMETIC_OP_COUNTS_TO_NANOMAGGY,/* Counts to SDSS nanomaggies. */
@@ -231,6 +231,38 @@ enum gal_arithmetic_operators
GAL_ARITHMETIC_OP_POOLMEAN, /* The pool-mean of desired pixels. */
GAL_ARITHMETIC_OP_POOLMEDIAN, /* The pool-median of desired pixels. */
+ /* Celestial coordinate conversions (names are clear, no comment). */
+ GAL_ARITHMETIC_OP_EQB1950_TO_EQJ2000,
+ GAL_ARITHMETIC_OP_EQB1950_TO_ECB1950,
+ GAL_ARITHMETIC_OP_EQB1950_TO_ECJ2000,
+ GAL_ARITHMETIC_OP_EQB1950_TO_GALACTIC,
+ GAL_ARITHMETIC_OP_EQB1950_TO_SUPERGALACTIC,
+ GAL_ARITHMETIC_OP_EQJ2000_TO_EQB1950,
+ GAL_ARITHMETIC_OP_EQJ2000_TO_ECB1950,
+ GAL_ARITHMETIC_OP_EQJ2000_TO_ECJ2000,
+ GAL_ARITHMETIC_OP_EQJ2000_TO_GALACTIC,
+ GAL_ARITHMETIC_OP_EQJ2000_TO_SUPERGALACTIC,
+ GAL_ARITHMETIC_OP_ECB1950_TO_EQB1950,
+ GAL_ARITHMETIC_OP_ECB1950_TO_EQJ2000,
+ GAL_ARITHMETIC_OP_ECB1950_TO_ECJ2000,
+ GAL_ARITHMETIC_OP_ECB1950_TO_GALACTIC,
+ GAL_ARITHMETIC_OP_ECB1950_TO_SUPERGALACTIC,
+ GAL_ARITHMETIC_OP_ECJ2000_TO_EQB1950,
+ GAL_ARITHMETIC_OP_ECJ2000_TO_EQJ2000,
+ GAL_ARITHMETIC_OP_ECJ2000_TO_ECB1950,
+ GAL_ARITHMETIC_OP_ECJ2000_TO_GALACTIC,
+ GAL_ARITHMETIC_OP_ECJ2000_TO_SUPERGALACTIC,
+ GAL_ARITHMETIC_OP_GALACTIC_TO_EQB1950,
+ GAL_ARITHMETIC_OP_GALACTIC_TO_EQJ2000,
+ GAL_ARITHMETIC_OP_GALACTIC_TO_ECB1950,
+ GAL_ARITHMETIC_OP_GALACTIC_TO_ECJ2000,
+ GAL_ARITHMETIC_OP_GALACTIC_TO_SUPERGALACTIC,
+ GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_EQB1950,
+ GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_EQJ2000,
+ GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_ECB1950,
+ GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_ECJ2000,
+ GAL_ARITHMETIC_OP_SUPERGALACTIC_TO_GALACTIC,
+
/* Counter for number of operators. */
GAL_ARITHMETIC_OP_LAST_CODE, /* Last code of the library operands. */
};
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 58c03350..35be2a5b 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -203,7 +203,14 @@ gal_wcs_coordsys_identify(struct wcsprm *inwcs);
struct wcsprm *
gal_wcs_coordsys_convert(struct wcsprm *inwcs, int coordsysid);
+void
+gal_wcs_coordsys_convert_points(int sys1, double *lng1_d, double *lat1_d,
+ int sys2, double *lng2_d, double *lat2_d,
+ size_t number);
+void
+gal_wcs_coordsys_sys1_ref_in_sys2(int sys1, int sys2, double *lng2,
+ double *lat2);
/*************************************************************
*********** Distortions ***********
diff --git a/lib/wcs.c b/lib/wcs.c
index 15deaa76..4e8b30a2 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -431,21 +431,22 @@ gal_wcs_read_fitsptr(fitsfile *fptr, int linearmatrix,
size_t hstartwcs,
if( strncmp(fullheader+i, "CRVAL1 = '", 11) == 0 )
fprintf(stderr, "WARNING: WCS Keyword values are not "
"numbers.\n\n"
- "WARNING: The values to the WCS-related keywords are "
- "enclosed in single-quotes. In the FITS standard "
- "this is how string values are stored, therefore "
- "WCSLIB is unable to read them AND WILL PUT ZERO IN "
- "THEIR PLACE (creating a wrong WCS in the output). "
- "Please update the respective keywords of the input "
- "to be numbers (see next line).\n\n"
+ "WARNING: The values to the WCS-related keywords "
+ "are enclosed in single-quotes. In the FITS "
+ "standard this is how string values are stored, "
+ "therefore WCSLIB is unable to read them AND WILL "
+ "PUT ZERO IN THEIR PLACE (creating a wrong WCS in "
+ "the output). Please update the respective keywords "
+ "of the input to be numbers (see next line).\n\n"
"WARNING: You can do this with Gnuastro's 'astfits' "
"program and the '--update' option. The minimal WCS "
- "keywords that need a numerical value are: 'CRVAL1', "
- "'CRVAL2', 'CRPIX1', 'CRPIX2', 'EQUINOX' and "
- "'CD%%_%%' (or 'PC%%_%%', where the %% are integers), "
- "please see the FITS standard, and inspect your FITS "
- "file to identify the full set of keywords that you "
- "need correct (for example PV%%_%% keywords).\n\n");
+ "keywords that need a numerical value are: "
+ "'CRVAL1', 'CRVAL2', 'CRPIX1', 'CRPIX2', 'EQUINOX' "
+ "and 'CD%%_%%' (or 'PC%%_%%', where the %% are "
+ "integers), please see the FITS standard, and "
+ "inspect your FITS file to identify the full set "
+ "of keywords that you need correct (for example "
+ "PV%%_%% keywords).\n\n");
}
/* CTYPE is a mandatory WCS keyword, so if it hasn't been given (its
@@ -917,13 +918,16 @@ gal_wcs_coordsys_identify(struct wcsprm *wcs)
/* Set the pole coordinates (current values taken from the WCSLIB
- manual.
- lng2p1: pole of input (1) system in output (2) system's logitude.
- lat2p1: pole of input (1) system in output (2) system's latitude.
- lng1p2: pole of output (2) system in input (1) system's longitude.
+ manual):
+ lng2p1: longitude (in system 2) of the pole of system 1.
+ lat2p1: latitude (in system 2) of the pole of system 1.
+ lng1p2: longitude (in system 1) of the pole of system 2.
Values from NED (inspired by WCSLIB manual's example).
https://ned.ipac.caltech.edu/coordinate_calculator
+ With the following constraints:
+ - The "Observation epoch" is the same as "Equinox".
+ - The "Position angle (East of North) is set to 0.0.
longi (deg) latit (deg) OUTPUT INPUT
----- ----- ------ -----
@@ -970,13 +974,13 @@ gal_wcs_coordsys_identify(struct wcsprm *wcs)
(------------, -----------) Supgalactic coords. of SupGalactic pole.
*/
static void
-wcs_coordsys_insys_pole_in_outsys(int insys, int outsys, double *lng2p1,
- double *lat2p1, double *lng1p2)
+wcs_coordsys_sys1_pole_in_sys2(int sys1, int sys2, double *lng2p1,
+ double *lat2p1, double *lng1p2)
{
- switch( insys )
+ switch( sys1 )
{
case GAL_WCS_COORDSYS_EQB1950:
- switch( outsys)
+ switch( sys2)
{
case GAL_WCS_COORDSYS_EQB1950:
*lng2p1=NAN; *lat2p1=NAN; *lng1p2=NAN;
return;
@@ -993,12 +997,12 @@ wcs_coordsys_insys_pole_in_outsys(int insys, int outsys,
double *lng2p1,
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
"fix the problem. The code '%d' isn't a recognized WCS "
- "coordinate system ID for 'outsys' (input EQB1950)", __func__,
- PACKAGE_BUGREPORT, outsys);
+ "coordinate system ID for 'sys2' (input EQB1950)",
+ __func__, PACKAGE_BUGREPORT, sys2);
}
break;
case GAL_WCS_COORDSYS_EQJ2000:
- switch( outsys)
+ switch( sys2)
{
case GAL_WCS_COORDSYS_EQB1950:
*lng2p1=359.68621044; *lat2p1=89.72178502; *lng1p2=180.31684301;
return;
@@ -1015,12 +1019,12 @@ wcs_coordsys_insys_pole_in_outsys(int insys, int
outsys, double *lng2p1,
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
"fix the problem. The code '%d' isn't a recognized WCS "
- "coordinate system ID for 'outsys' (input EQJ2000)", __func__,
- PACKAGE_BUGREPORT, outsys);
+ "coordinate system ID for 'sys2' (input EQJ2000)",
+ __func__, PACKAGE_BUGREPORT, sys2);
}
break;
case GAL_WCS_COORDSYS_ECB1950:
- switch( outsys)
+ switch( sys2)
{
case GAL_WCS_COORDSYS_EQB1950:
*lng2p1=270.00000000; *lat2p1=66.55421111; *lng1p2=90.000000000;
return;
@@ -1037,12 +1041,12 @@ wcs_coordsys_insys_pole_in_outsys(int insys, int
outsys, double *lng2p1,
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
"fix the problem. The code '%d' isn't a recognized WCS "
- "coordinate system ID for 'outsys' (input ECB1950)", __func__,
- PACKAGE_BUGREPORT, outsys);
+ "coordinate system ID for 'sys2' (input ECB1950)",
+ __func__, PACKAGE_BUGREPORT, sys2);
}
break;
case GAL_WCS_COORDSYS_ECJ2000:
- switch( outsys)
+ switch( sys2)
{
case GAL_WCS_COORDSYS_EQB1950:
*lng2p1=270.00099211; *lat2p1=66.56069675; *lng1p2=90.699521110;
return;
@@ -1059,12 +1063,12 @@ wcs_coordsys_insys_pole_in_outsys(int insys, int
outsys, double *lng2p1,
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
"fix the problem. The code '%d' isn't a recognized WCS "
- "coordinate system ID for 'outsys' (input ECJ2000)", __func__,
- PACKAGE_BUGREPORT, outsys);
+ "coordinate system ID for 'sys2' (input ECJ2000)",
+ __func__, PACKAGE_BUGREPORT, sys2);
}
break;
case GAL_WCS_COORDSYS_GALACTIC:
- switch( outsys)
+ switch( sys2)
{
case GAL_WCS_COORDSYS_EQB1950:
*lng2p1=192.25000000; *lat2p1=27.40000000; *lng1p2=123.00000000;
return;
@@ -1081,12 +1085,12 @@ wcs_coordsys_insys_pole_in_outsys(int insys, int
outsys, double *lng2p1,
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
"fix the problem. The code '%d' isn't a recognized WCS "
- "coordinate system ID for 'outsys' (input GALACTIC)", __func__,
- PACKAGE_BUGREPORT, outsys);
+ "coordinate system ID for 'sys2' (input GALACTIC)",
+ __func__, PACKAGE_BUGREPORT, sys2);
}
break;
case GAL_WCS_COORDSYS_SUPERGALACTIC:
- switch( outsys)
+ switch( sys2)
{
case GAL_WCS_COORDSYS_EQB1950:
*lng2p1=283.18940711; *lat2p1=15.64407736; *lng1p2=26.731537070;
return;
@@ -1103,17 +1107,186 @@ wcs_coordsys_insys_pole_in_outsys(int insys, int
outsys, double *lng2p1,
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
"fix the problem. The code '%d' isn't a recognized WCS "
- "coordinate system ID for 'outsys' (input SUPERGALACTIC)",
__func__,
- PACKAGE_BUGREPORT, outsys);
+ "coordinate system ID for 'sys2' (input SUPERGALACTIC)",
+ __func__, PACKAGE_BUGREPORT, sys2);
}
break;
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
"fix the problem. The code '%d' isn't a recognized WCS "
- "coordinate system ID for 'insys'", __func__,
- PACKAGE_BUGREPORT, insys);
+ "coordinate system ID for 'sys1'", __func__,
+ PACKAGE_BUGREPORT, sys1);
}
+}
+
+
+
+
+
+/* Return the longitude (RA in equatorial) and latitude (Dec in equatorial)
+ coordinates of the reference point (on the equator) of system 1 in the
+ 2nd system. For example, if you want the galactic center's RA and Dec:
+ 'sys1=GAL_WCS_COORDSYS_GALACTIC' and 'sys2=GAL_WCS_COORDSYS_EQJ2000'.
+
+ Similar to 'wcs_coordsys_sys1_pole_in_sys2', the numbers reported here
+ come from NED. Here, the following extra settings were set (they were
+ not recorded in the comments when 'wcs_coordsys_sys1_pole_in_sys2' was
+ defined).
+
+ - The "Observation epoch" has been set to the same year as the
+ equinox.
+ - The "Position angle (East of North)" is set to 0.
+*/
+void
+gal_wcs_coordsys_sys1_ref_in_sys2(int sys1, int sys2, double *lng2,
+ double *lat2)
+{
+ switch( sys1 )
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ switch( sys2 )
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ *lng2=NAN; *lat2=NAN; return;
+ case GAL_WCS_COORDSYS_EQJ2000:
+ *lng2=0.64069119; *lat2=0.27839567; return;
+ case GAL_WCS_COORDSYS_ECB1950:
+ *lng2=0.00000000; *lat2=0.00000000; return;
+ case GAL_WCS_COORDSYS_ECJ2000:
+ *lng2=0.69855979; *lat2=0.00057803; return;
+ case GAL_WCS_COORDSYS_GALACTIC:
+ *lng2=97.74216087; *lat2=-60.18102400; return;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ *lng2=293.11549599; *lat2=12.69249218; return;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'sys2' (input EQB1950)",
+ __func__, PACKAGE_BUGREPORT, sys2);
+ }
+ break;
+
+ case GAL_WCS_COORDSYS_EQJ2000:
+ switch( sys2 )
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ *lng2=359.35927364; *lat2=-0.27832018; return;
+ case GAL_WCS_COORDSYS_EQJ2000:
+ *lng2=NAN; *lat2=NAN; return;
+ case GAL_WCS_COORDSYS_ECB1950:
+ *lng2=359.30143791; *lat2=-0.00041556; return;
+ case GAL_WCS_COORDSYS_ECJ2000:
+ *lng2=0.0000000000; *lat2=0.000000000; return;
+ case GAL_WCS_COORDSYS_GALACTIC:
+ *lng2=96.33723581; *lat2=-60.18845577; return;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ *lng2=292.65879220; *lat2=13.23092157; return;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'sys2' (input EQJ2000)",
+ __func__, PACKAGE_BUGREPORT, sys2);
+ }
+ break;
+
+ case GAL_WCS_COORDSYS_ECB1950:
+ switch( sys2 )
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ *lng2=0.0000000000; *lat2=0.000000000; return;
+ case GAL_WCS_COORDSYS_EQJ2000:
+ *lng2=0.64069119; *lat2=0.27839567; return;
+ case GAL_WCS_COORDSYS_ECB1950:
+ *lng2=NAN; *lat2=NAN; return;
+ case GAL_WCS_COORDSYS_ECJ2000:
+ *lng2=0.69855979; *lat2=0.00057803; return;
+ case GAL_WCS_COORDSYS_GALACTIC:
+ *lng2=97.74216087; *lat2=-60.18102400; return;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ *lng2=293.11549599; *lat2=12.69249218; return;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'sys2' (input ECB1950)",
+ __func__, PACKAGE_BUGREPORT, sys2);
+ }
+ break;
+
+ case GAL_WCS_COORDSYS_ECJ2000:
+ switch( sys2 )
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ *lng2=359.35927364; *lat2=-0.27832018; return;
+ case GAL_WCS_COORDSYS_EQJ2000:
+ *lng2=0.0000000000; *lat2=0.000000000; return;
+ case GAL_WCS_COORDSYS_ECB1950:
+ *lng2=359.30143791; *lat2=-0.00041556; return;
+ case GAL_WCS_COORDSYS_ECJ2000:
+ *lng2=NAN; *lat2=NAN; return;
+ case GAL_WCS_COORDSYS_GALACTIC:
+ *lng2=96.33723581; *lat2=-60.18845577; return;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ *lng2=292.65879220; *lat2=13.23092157 ; return;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'sys2' (input ECJ2000)",
+ __func__, PACKAGE_BUGREPORT, sys2);
+ }
+ break;
+
+ case GAL_WCS_COORDSYS_GALACTIC:
+ switch( sys2 )
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ *lng2=265.61084403; *lat2=-28.91679035; return;
+ case GAL_WCS_COORDSYS_EQJ2000:
+ *lng2=266.40506655; *lat2=-28.93616241; return;
+ case GAL_WCS_COORDSYS_ECB1950:
+ *lng2=266.14096542; *lat2=-5.52979411; return;
+ case GAL_WCS_COORDSYS_ECJ2000:
+ *lng2=266.83958692; *lat2=-5.53630157; return;
+ case GAL_WCS_COORDSYS_GALACTIC:
+ *lng2=NAN; *lat2=NAN; return;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ *lng2=185.78610785; *lat2=42.31028736; return;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'sys2' (input GALACTIC)",
+ __func__, PACKAGE_BUGREPORT, sys2);
+ }
+ break;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ switch( sys2 )
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ *lng2=41.35517115; *lat2=59.32090820; return;
+ case GAL_WCS_COORDSYS_EQJ2000:
+ *lng2=42.30997710; *lat2=59.52821263; return;
+ case GAL_WCS_COORDSYS_ECB1950:
+ *lng2=59.54957645; *lat2=40.91184040; return;
+ case GAL_WCS_COORDSYS_ECJ2000:
+ *lng2=60.24554909; *lat2=40.91764242; return;
+ case GAL_WCS_COORDSYS_GALACTIC:
+ *lng2=137.37000000; *lat2=0.00000000; return;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ *lng2=NAN; *lat2=NAN; return;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'sys2' (input SUPERGALACTIC)",
+ __func__, PACKAGE_BUGREPORT, sys2);
+ }
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'sys1'", __func__,
+ PACKAGE_BUGREPORT, sys1);
+ }
}
@@ -1176,8 +1349,8 @@ gal_wcs_coordsys_convert(struct wcsprm *wcs, int
outcoordsys)
/* Find the necessary pole coordinates. Note that we have already
accounted for the fact that the input and output coordinate systems
may be the same above, so the NaN outputs will never occur here. */
- wcs_coordsys_insys_pole_in_outsys(incoordsys, outcoordsys,
- &lng2p1, &lat2p1, &lng1p2);
+ wcs_coordsys_sys1_pole_in_sys2(incoordsys, outcoordsys,
+ &lng2p1, &lat2p1, &lng1p2);
/* Find the necessary CTYPE names of the output. */
wcs_coordsys_ctypes(outcoordsys, &clng, &clat, &radesys, &equinox);
@@ -1219,6 +1392,111 @@ gal_wcs_coordsys_convert(struct wcsprm *wcs, int
outcoordsys)
+/* Convert the coordinates of a series of points from one celestial
+ coordinate system to another.
+
+ This is based on the equations in this Wikipedia article:
+
https://en.wikipedia.org/wiki/Galactic_coordinate_system#Conversion_between_equatorial_and_galactic_coordinates
+ A mathematical proof of this equation is present in this StackExchange
+ answer:
+
https://physics.stackexchange.com/questions/88663/converting-between-galactic-and-ecliptic-coordinates
+
+
+ The angle names are taken from this solution (for example "BK" is shown
+ with the variable 'bk_r', since it is in radians). Intermediate angles
+ (that the user never seen by the caller) are in radians and have a '_r'
+ suffix. The angles that the caller gives and recieves are in degrees
+ with a '_d' suffix.
+ */
+void
+gal_wcs_coordsys_convert_points(int sys1, double *lng1_d, double *lat1_d,
+ int sys2, double *lng2_d, double *lat2_d,
+ size_t number)
+{
+ size_t i;
+ double coslat1, coslat2;
+ double lngdiff, coslngdiff, coslat1p2;
+ double lng1_r, lat1_r, lng2_r, lat2_r;
+ double lng2p1diff_r, lng2p1diff_sin_r, lng2p1diff_cos_r;
+
+ /* Format: 'aaaNcM':
+ - aaa either 'lng' (for longitude) or 'lat' (for latitude).
+ - N either 1 or 2 (for the coordinate system of 'aaa').
+ - c either 'p' (for pole) or 'r' (for reference point).
+ - M either 1 or 2 (for the coordinate system of 'c'). */
+ double lng1p2_d, lat1p2_d, lng2p1_d;
+ double lng1p2_r, lat1p2_r, lng2p1_r;
+
+ /* In case the input and output coordinate systems are the same, just
+ copy the first into the second coordinates and return (no need to do
+ any conversion). */
+ if(sys1==sys2)
+ {
+ for(i=0;i<number;++i) {lng2_d[i]=lng1_d[1]; lat2_d[i]=lat1_d[1];}
+ return;
+ }
+
+ /* Get the second system's pole and reference point (on its equator) in
+ the first system. */
+ wcs_coordsys_sys1_pole_in_sys2(sys2, sys1, &lng1p2_d, &lat1p2_d,
+ &lng2p1_d);
+ lng1p2_r=lng1p2_d*M_PI/180.0f;
+ lat1p2_r=lat1p2_d*M_PI/180.0f;
+ lng2p1_r=lng2p1_d*M_PI/180.0f;
+
+ /* Loop over all the coordinates. */
+ coslat1p2=cos(lat1p2_r);
+ for(i=0;i<number;++i)
+ {
+ /* Convert the input coordinates into radians. */
+ lng1_r = lng1_d[i] * M_PI / 180.0f;
+ lat1_r = lat1_d[i] * M_PI / 180.0f;
+
+ /* The latitude in the output coordinate is easy to calculate: */
+ lngdiff=lng1_r-lng1p2_r;
+ coslngdiff=cos(lngdiff);
+ lat2_r = asin( sin(lat1p2_r)*sin(lat1_r)
+ + coslat1p2*cos(lat1_r)*coslngdiff );
+
+ /* We can now calculate the angle between 'PG' and 'GR' ('122.9 - l'
+ in the webpage above), we'll call it 'pggr_r' here. But there is a
+ problem: we need both the sine and cosine of 'pggr_r' to get its
+ position across the whole 0 to 360 degrees ('asin()' only returns
+ a value from -pi/2 to pi/2 and 'acos()' only returns a value
+ between 0 to pi. */
+ coslat1 = cos(lat1_r);
+ coslat2 = cos(lat2_r);
+ lng2p1diff_sin_r = asin( coslat1 * sin(lngdiff) / coslat2 );
+ lng2p1diff_cos_r = acos( ( coslat1p2*sin(lat1_r)
+ - sin(lat1p2_r)*coslat1*coslngdiff )
+ / coslat2 );
+
+ /* We have assumed that 'lng2p1diff_r = lng2p1_r - lng2_r'; so
+ 'lng2_r = lng2p1_r - lng2p1diff_r' (when
+ 'lng2p1diff_sin_r>0'). When 'lng2p1diff_sin_r<0', the cosine needs
+ to be negative (so in the end it becomes 'lng2_r = lng2p1_r +
+ lng2p1diff_r'. */
+ lng2p1diff_r = ( lng2p1diff_sin_r<0
+ ? -1*lng2p1diff_cos_r
+ : lng2p1diff_cos_r );
+ lng2_r = lng2p1_r - lng2p1diff_r;
+
+ /* In case the longitude is larger than 2pi, then we should subtract
+ 2*pi from it and if it is negative, we should add 2pi to it. */
+ if(lng2_r>2*M_PI) lng2_r=lng2_r-2*M_PI;
+ else if (lng2_r<0) lng2_r+=2*M_PI;
+
+ /* Write the values in the output as degrees. */
+ lng2_d[i]=lng2_r*180.0f/M_PI;
+ lat2_d[i]=lat2_r*180.0f/M_PI;
+ }
+}
+
+
+
+
+
+
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master de51acf8: Library (wcs.h): convert points from one coordsys to another,
Mohammad Akhlaghi <=