[gnuastro-commits] master 042b78f 2/4: Pixel scale function added, angular distance renamed

[Mohammad Akhlaghi]

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

    Pixel scale function added, angular distance renamed
    
    `gal_wcs_pixel_scale_deg' is a new function added to Gnuastro's WCS
    library. As the name suggests, it will give the pixel scale of the input
    WCS structure in degrees. It has also been documented in the book.
    
    `gal_wcs_angular_distance_deg' is the new name of the old function called
    `gal_wcs_angular_distance'. That function would take the angles in radians
    and also return its value in radians. Radians are very inconvenient and
    uncommon in celestial coordinates, so its input and output was modified to
    accept angles in degrees and also return the angular distance in
    degrees. To be more explicit about the units, the `_deg' suffix was added
    to the name.
---
 doc/gnuastro.texi  |   20 ++++++++++------
 lib/gnuastro/wcs.h |    5 +++-
 lib/wcs.c          |   65 ++++++++++++++++++++++++++++------------------------
 3 files changed, 52 insertions(+), 38 deletions(-)

diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 5effba9..98fc4a0 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -17072,22 +17072,28 @@ Similar to @code{gal_wcs_xy_array_to_radec} but 
convert the world
 coordinates in @code{radec} to image coordinates in @code{xy}.
 @end deftypefun
 
address@hidden double gal_wcs_angular_distance (double @code{r1}, double 
@code{d1}, double @code{r2}, double @code{d2})
-Return the angular distance between a point located at (@code{r1},
address@hidden) to (@code{r2}, @code{d2}). The distance (along a great circle)
-on a sphere between two points is calculated with the equation below. Since
-the pixel sides are usually very small, we won't be using the direct
-formula:
address@hidden double gal_wcs_angular_distance_deg (double @code{r1}, double 
@code{d1}, double @code{r2}, double @code{d2})
+Return the angular distance (in degrees) between a point located at
+(@code{r1}, @code{d1}) to (@code{r2}, @code{d2}). All input coordinates are
+in degrees. The distance (along a great circle) on a sphere between two
+points is calculated with the equation below. Since the pixel sides are
+usually very small, we won't be using the direct formula:
 
 @dispmath {\cos(d)=\sin(d_1)\sin(d_2)+\cos(d_1)\cos(d_2)\cos(r_1-r_2)}
 
 Here, we will be using the
 @url{https://en.wikipedia.org/wiki/Haversine_formula, Haversine formula}
-which better considering floating point errors:
+which is better considering floating point errors:
 
 @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
 
address@hidden void gal_wcs_pixel_scale_deg (struct wcsprm @code{*wcs}, double 
@code{*dx}, double @code{*dy})
+Return the pixel scale of the WCS structure @code{wcs} for the two image
+dimensions @code{dx} and @code{dy} in degrees. Here, the X-axis is the
+first FITS axis, or the horizontal axis when viewed in SAO ds9.
address@hidden deftypefun
+
 @deftypefun double gal_wcs_pixel_area_arcsec2 (struct wcsprm @code{*wcs})
 Given the WCS structure @code{wcs} calculate the pixel area in
 arcsec-squared.
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 1041a8f..aa81d29 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -57,7 +57,10 @@ gal_wcs_radec_array_to_xy(struct wcsprm *wcs, double *radec, 
double *xy,
                           size_t number, size_t width);
 
 double
-gal_wcs_angular_distance(double r1, double d1, double r2, double d2);
+gal_wcs_angular_distance_deg(double r1, double d1, double r2, double d2);
+
+void
+gal_wcs_pixel_scale_deg(struct wcsprm *wcs, double *dx, double *dy);
 
 double
 gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs);
diff --git a/lib/wcs.c b/lib/wcs.c
index 5d17828..0c716df 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -126,14 +126,40 @@ gal_wcs_radec_array_to_xy(struct wcsprm *wcs, double 
*radec, double *xy,
    floating point errors (from Wikipedia:)
 
    sin^2(distance)/2=sin^2( (d1-d2)/2 )+cos(d1)*cos(d2)*sin^2( (r1-r2)/2 )
+
+   Inputs and outputs are all in degrees.
 */
 double
-gal_wcs_angular_distance(double r1, double d1, double r2, double d2)
+gal_wcs_angular_distance_deg(double r1, double d1, double r2, double d2)
 {
-  double a=sin( (d1-d2)/2 );
-  double b=sin( (r1-r2)/2 );
+  /* Convert degrees to radians. */
+  double r1r=r1*M_PI/180, d1r=d1*M_PI/180;
+  double r2r=r2*M_PI/180, d2r=d2*M_PI/180;
+
+  /* To make things easier to read: */
+  double a=sin( (d1r-d2r)/2 );
+  double b=sin( (r1r-r2r)/2 );
 
-  return 2*asin( sqrt( a*a + cos(d1)*cos(d2)*b*b) );
+  /* Return the result: */
+  return 2*asin( sqrt( a*a + cos(d1r)*cos(d2r)*b*b) ) * 180/M_PI;
+}
+
+
+
+
+/* Return the pixel scale of the image in both dimentions in degrees. */
+void
+gal_wcs_pixel_scale_deg(struct wcsprm *wcs, double *dx, double *dy)
+{
+  double radec[6], xy[]={0,0,1,0,0,1};
+
+  /* Get the RA and Dec of the bottom left, bottom right and top left
+     sides of the first pixel in the image. */
+  gal_wcs_xy_array_to_radec(wcs, xy, radec, 3, 2);
+
+  /* Calculate the distances and convert back to degrees: */
+  *dx = gal_wcs_angular_distance_deg(radec[0], radec[1], radec[2], radec[3]);
+  *dy = gal_wcs_angular_distance_deg(radec[0], radec[1], radec[4], radec[5]);
 }
 
 
@@ -148,32 +174,11 @@ gal_wcs_angular_distance(double r1, double d1, double r2, 
double d2)
 double
 gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs)
 {
-  double xy[]={0,0,1,0,0,1};
-  double st, *d, *df, radec[6];
-
-  /* Get the RA and Dec of the bottom left, bottom right and top left
-     sides of the first pixel in the image. */
-  gal_wcs_xy_array_to_radec(wcs, xy, radec, 3, 2);
-
-  /* Covert the RA and dec values to radians for easy calculation: */
-  df=(d=radec)+6; do *d++ *= M_PI/180.0f; while(d<df);
-
-  /* For a check:
-  printf("\n\nAlong first axis: %g\nAlong second axis: %g\n\n",
-         ( angulardistance(radec[0], radec[1], radec[2], radec[3])
-           *180/M_PI*3600 ),
-         ( angulardistance(radec[0], radec[1], radec[4], radec[5])
-           *180/M_PI*3600 ) );
-  */
-
-  /* Get the area in stradians. */
-  st= ( gal_wcs_angular_distance(radec[0], radec[1], radec[2], radec[3]) *
-        gal_wcs_angular_distance(radec[0], radec[1], radec[4], radec[5]) );
+  double dx, dy;
 
-  /* Convert the stradians to arcsec^2:
+  /* Get the pixel scales along each axis in degrees. */
+  gal_wcs_pixel_scale_deg(wcs, &dx, &dy);
 
-     1deg^2 = (180/PI)^2 * 1stradian.
-     1arcsec^2 = (3600*3600) * 1degree^2
-   */
-  return st*180.0f*180.0f*3600.0f*3600.0f/(M_PI*M_PI);
+  /* Return the result */
+  return dx * dy * 3600.0f * 3600.0f;
 }