From d81a364b0f8e7e847dae9274ee89fe6cdcffaa0e Mon Sep 17 00:00:00 2001
From: William Joye <wjoye@cfa.harvard.edu>
Date: Wed, 7 Nov 2018 14:37:21 -0500
Subject: fix save fits in case of wcs axis>2

---
 tksao/frame/frsave.C | 14 ++++----------
 tksao/frame/wcsast.C | 50 ++++++++++++++++++++++++++++++++++++++++++++++++++
 tksao/frame/wcsast.h |  1 +
 3 files changed, 55 insertions(+), 10 deletions(-)

diff --git a/tksao/frame/frsave.C b/tksao/frame/frsave.C
index 4a3aee1..467239b 100644
--- a/tksao/frame/frsave.C
+++ b/tksao/frame/frsave.C
@@ -10,6 +10,7 @@
 #include "outchannel.h"
 #include "outsocket.h"
 #include "sigbus.h"
+#include "wcsast.h"
 
 void Base::saveFits(OutFitsStream& str)
 {
@@ -435,17 +436,10 @@ void FrameBase::saveFitsResampleKeyword(OutFitsStream& str, FitsHead& dst)
 
   // ast
   AstFrameSet* ast = (AstFrameSet*)astCopy(currentContext->fits->ast_);
-
   Matrix mx = currentContext->fits->imageToRef * refToWidget * cc;
-  double ss[] = {mx.matrix(0,0),mx.matrix(1,0),
-		 mx.matrix(0,1),mx.matrix(1,1)};
-  double tt[] = {mx.matrix(2,0),mx.matrix(2,1)};
-    
-  AstMatrixMap* mm = astMatrixMap(2, 2, 0, ss, "");
-  AstShiftMap* sm = astShiftMap(2, tt, "");
-  AstCmpMap* cmp = astCmpMap(mm, sm, 1, "");
-
-  astRemapFrame(ast, AST__BASE, cmp);
+  AstCmpMap* cmp = wcsMatrixMap(ast, mx);
+  if (cmp)
+    astRemapFrame(ast, AST__BASE, cmp);
 
   // write to channel
   if (!astWrite(chan, ast)) {
diff --git a/tksao/frame/wcsast.C b/tksao/frame/wcsast.C
index 5ae4ac7..d1b75a6 100644
--- a/tksao/frame/wcsast.C
+++ b/tksao/frame/wcsast.C
@@ -438,4 +438,54 @@ AstWinMap* wcsWinMap(AstFrameSet* ast, Vector& ll, Vector& ur, Vector& rr)
     }
     break;
   }
+
+  return NULL;
+}
+
+AstCmpMap* wcsMatrixMap(AstFrameSet* ast, Matrix& mx)
+{
+  AstMatrixMap* mm =NULL;
+  AstShiftMap* sm =NULL;
+
+  int naxes = astGetI(ast,"Naxes");
+  switch (naxes) {
+  case 1:
+    // error
+    return NULL;
+  case 2:
+    {
+      double ss[] = {mx.matrix(0,0),mx.matrix(1,0),
+		     mx.matrix(0,1),mx.matrix(1,1)};
+      double tt[] = {mx.matrix(2,0),mx.matrix(2,1)};
+      mm = astMatrixMap(2, 2, 0, ss, "");
+      sm = astShiftMap(2, tt, "");
+    }
+    break;
+  case 3:
+    {
+      double ss[] = {mx.matrix(0,0),mx.matrix(1,0),0,
+		     mx.matrix(0,1),mx.matrix(1,1),0,
+		     0,0,1};
+      double tt[] = {mx.matrix(2,0),mx.matrix(2,1),0};
+      mm = astMatrixMap(3, 3, 0, ss, "");
+      sm = astShiftMap(3, tt, "");
+    }
+    break;
+  case 4:
+    {
+      double ss[] = {mx.matrix(0,0),mx.matrix(1,0),0,0,
+		     mx.matrix(0,1),mx.matrix(1,1),0,0,
+		     0,0,1,0,
+		     0,0,0,1};
+      double tt[] = {mx.matrix(2,0),mx.matrix(2,1),0,0};
+      mm = astMatrixMap(4, 4, 0, ss, "");
+      sm = astShiftMap(4, tt, "");
+    }
+    break;
+  }
+
+  if (mm && sm)
+    return astCmpMap(mm, sm, 1, "");
+  else
+    return NULL;
 }
diff --git a/tksao/frame/wcsast.h b/tksao/frame/wcsast.h
index e1c85fe..79dc101 100644
--- a/tksao/frame/wcsast.h
+++ b/tksao/frame/wcsast.h
@@ -23,3 +23,4 @@ double wcsAngle(AstFrameSet*, const Vector&, const Vector&, const Vector&);
 double wcsAxAngle(AstFrameSet*, const Vector&, const Vector&);
 
 AstWinMap* wcsWinMap(AstFrameSet*, Vector&, Vector&, Vector&);
+AstCmpMap* wcsMatrixMap(AstFrameSet* ast, Matrix&);
-- 
cgit v0.12