[gnuastro-commits] master edd1686d: Library (wcs.h): gal_wcs_to_cd accounts for altlin==3

[Mohammad Akhlaghi]

branch: master
commit edd1686d551a4fae64175624f9448b5a98325452
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Library (wcs.h): gal_wcs_to_cd accounts for altlin==3
    
    Until now, this function would only consider scenarios where the input file
    had a single rotation matrix (in particular, 'CDi_j' and 'PCi_j'). However,
    it can happen that an input FITS file has both! In this case, it would
    abort with a bug notice.
    
    With this commit, when both representations are present, it will check if
    the values are the same or not. If they are, then it will just re-write
    'altlin' so WCSLIB uses the 'CDi_j' representation (as requested from this
    function). But if they are not the same, then it will abort with an
    informative error to let the user know about the problem and how they can
    fix it.
---
 NEWS                         |  2 ++
 THANKS                       |  1 +
 doc/announce-acknowledge.txt |  1 +
 lib/wcs.c                    | 33 ++++++++++++++++++++++++++++-----
 4 files changed, 32 insertions(+), 5 deletions(-)

diff --git a/NEWS b/NEWS
index 3780a020..b4b2698b 100644
--- a/NEWS
+++ b/NEWS
@@ -127,6 +127,8 @@ See the end of the file for license conditions.
               blanks present.
   bug #64420: Error message of wrongly specified HDU doesn't give the
               problematic option name; reported by Aaron Watkins.
+  bug #64431: gal_wcs_to_cd does not account for altlin==3. Reported by
+              Colin Orion Chandler.
 
 
 
diff --git a/THANKS b/THANKS
index 91312fc8..a9713424 100644
--- a/THANKS
+++ b/THANKS
@@ -43,6 +43,7 @@ support in Gnuastro. The list is ordered alphabetically (by 
family name).
     Carlos Morales Socorro               cmorsoc@gmail.com
     Christopher Willmer                  cnaw@as.arizona.edu
     Clotilde Laigle                      laigle@iap.fr
+    Colin Orion Chandler                 coc123@uw.edu
     Craig Gordon                         craig.a.gordon@nasa.gov
     David Shupe                          shupe@ipac.caltech.edu
     David Valls-Gabaud                   david.valls-gabaud@obspm.fr
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 9d7a6bdd..3765fb41 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -2,6 +2,7 @@ Alphabetically ordered list to acknowledge in the next release.
 
 Aaron Watkins
 Agata Rożek
+Colin Orion Chandler
 Irene Pintos Castro
 Rashid Yaaqib
 Raul Infante-Sainz
diff --git a/lib/wcs.c b/lib/wcs.c
index 542ec8b8..1c8fc697 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -1771,7 +1771,8 @@ gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs)
 void
 gal_wcs_to_cd(struct wcsprm *wcs)
 {
-  size_t i, j, naxis;
+  size_t i, j, n;
+  double er=1e-10;
 
   /* If there is on WCS, then don't do anything. */
   if(wcs==NULL) return;
@@ -1779,7 +1780,7 @@ gal_wcs_to_cd(struct wcsprm *wcs)
   /* 'wcs->altlin' identifies which rotation element is being used (PCi_j,
      CDi_J or CROTAi). For PCi_j, the first bit will be 1 (==1), for CDi_j,
      the second bit is 1 (==2) and for CROTAi, the third bit is 1 (==4). */
-  naxis=wcs->naxis;
+  n=wcs->naxis;
   switch(wcs->altlin)
     {
    /* PCi_j: Convert it to CDi_j. */
@@ -1791,10 +1792,10 @@ gal_wcs_to_cd(struct wcsprm *wcs)
          multiply the PC and CDELT matrices, then set cdelt=1 in all
          dimensions. This is actually how WCSLIB reads a FITS header with
          only a CD matrix. */
-      for(i=0;i<naxis;++i)
+      for(i=0;i<n;++i)
         {
-          for(j=0;j<naxis;++j)
-            wcs->cd[i*naxis+j] = wcs->pc[i*naxis+j] *= wcs->cdelt[i];
+          for(j=0;j<n;++j)
+            wcs->cd[i*n+j] = wcs->pc[i*n+j] *= wcs->cdelt[i];
           wcs->cdelt[i]=1;
         }
 
@@ -1805,6 +1806,28 @@ gal_wcs_to_cd(struct wcsprm *wcs)
     /* CDi_j: No need to do any conversion. */
     case 2: return; break;
 
+    /* Both PCi_j and CDi_j are present! If they are the same (within
+       floating point errors), we'll just set WCS to use the 'CD' matrix
+       (as demanded from this function). If they are not the same, then
+       print an error. */
+    case 3:
+      for(i=0;i<n;++i)
+        for(j=0;j<n;++j)
+          if(wcs->cd[i*n+j] - wcs->pc[i*n+j] * wcs->cdelt[i] > er)
+            error(EXIT_FAILURE, 0, "%s: the given WCS has both the "
+                  "CDi_j and PCi_j+CDELTi conventions for defining "
+                  "the rotation and scale. However, they do not match! "
+                  "Please inspect the file and remove the wrong set of "
+                  "keywords (you can use 'astfits file.fits "
+                  "--delete=KEYNAME' to delete the keyword 'KEYNAME'; "
+                  "and you can call '--delete' multiple times). For "
+                  "more on the definition of the different "
+                  "representations, see the FITS standard: "
+                  "https://fits.gsfc.nasa.gov/fits_standard.html";,
+                  __func__);
+      wcs->altlin=2;
+      break;
+
     /* CROTAi: not yet supported. */
     case 4:
       error(0, 0, "%s: WARNING: Conversion of 'CROTAi' keywords to the CD "