fix savefits wcs for dss

3 files changed, 83 insertions, 159 deletions

diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C
index 889befd..474326d 100644
--- a/tksao/frame/fitsimage.C
+++ b/tksao/frame/fitsimage.C
@@ -84,6 +84,7 @@ FitsImage::FitsImage(Context* cx, Tcl_Interp* pp)
 
   manageWCS_ =1;
   ast_ =NULL;
+  encoding_ =NULL;
   wcs_ =NULL;
   wcsNaxes_ =NULL;
 
@@ -159,6 +160,8 @@ FitsImage::~FitsImage()
   if (manageWCS_) {
     if (ast_)
       astAnnul(ast_);
+    if (encoding_)
+      delete [] encoding_;
 
     if (wcs_)
       delete [] wcs_;
@@ -1062,6 +1065,9 @@ void FitsImage::initWCS(FitsHead* hd)
     if (ast_)
       astAnnul(ast_);
     ast_ =NULL;
+    if (encoding_)
+      delete [] encoding_;
+    encoding_ =NULL;
 
     if (wcs_)
       delete [] wcs_;
@@ -1104,6 +1110,7 @@ void FitsImage::initWCS(FitsHead* hd)
       while (sptr) {
 	if (sptr == this) {
 	  ast_ = ptr->ast_;
+	  encoding_ = ptr->encoding_;
 
 	  wcs_ = ptr->wcs_;
 	  wcsNaxes_ = ptr->wcsNaxes_;
@@ -1134,6 +1141,9 @@ void FitsImage::initWCS(FitsHead* hd)
   if (ast_)
     astAnnul(ast_);
   ast_ =NULL;
+  if (encoding_)
+    delete [] encoding_;
+  encoding_ =NULL;
   
   ast_ = fits2ast(hd);
   if (!ast_)
@@ -3464,6 +3474,12 @@ AstFrameSet* FitsImage::fits2ast(FitsHead* hd)
       astClearStatus;
   }
 
+  // encoding
+  // must come after astPutsFits and before astRead
+  const char* encode = astGetC(chan, "Encoding");
+  if (encode)
+    encoding_ = dupstr(encode);
+
   // we may have an error, just reset
   astClearStatus;
   astClear(chan, "Card");
diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h
index 1acc5bf..5c4fe87 100644
--- a/tksao/frame/fitsimage.h
+++ b/tksao/frame/fitsimage.h
@@ -134,6 +134,7 @@ class FitsImage {
 
  public:
   AstFrameSet* ast_;  // ast frameset;
+  const char* encoding_; // ast encoding of original header
 
  private:
   char* root(const char*);
diff --git a/tksao/frame/frsave.C b/tksao/frame/frsave.C
index b0f854e..4a3aee1 100644
--- a/tksao/frame/frsave.C
+++ b/tksao/frame/frsave.C
@@ -376,57 +376,42 @@ void FrameBase::saveFitsResampleKeyword(OutFitsStream& str, FitsHead& dst)
   FitsHead* src = currentContext->fits->head();
   Vector center = Vector(options->width, options->height)/2.;
 
+  // center mx
+  Matrix cc =
+    Translate(-center) *
+    Translate(1,0) *
+    FlipY() *
+    Translate(center);
+
   // OBJECT
   char* object = src->getString("OBJECT");
   if (object)
     dst.appendString("OBJECT", object, NULL);
 
-  // DATE-OBS
-  char* date = src->getString("DATE");
-  if (date)
-    dst.appendString("DATE", date, NULL);
-
-  char* dateobs = src->getString("DATE-OBS");
-  if (dateobs)
-    dst.appendString("DATE-OBS", dateobs, NULL);
-
-  char* timeobs = src->getString("TIME-OBS");
-  if (timeobs)
-    dst.appendString("TIME-OBS", timeobs, NULL);
-
-  char* dateend = src->getString("DATE-END");
-  if (dateend)
-    dst.appendString("DATE-END", dateend, NULL);
-
-  char* timeend = src->getString("TIME-END");
-  if (timeend)
-    dst.appendString("TIME-END", timeend, NULL);
+  // LTMV
+  if (currentContext->fits->hasLTMV()) {
+    Matrix ltmv = currentContext->fits->physicalToRef * refToWidget * cc;
+    dst.appendReal("LTM1_1", ltmv[0][0], 9, NULL);
+    dst.appendReal("LTM1_2", ltmv[0][1], 9, NULL);
+    dst.appendReal("LTM2_1", ltmv[1][0], 9, NULL);
+    dst.appendReal("LTM2_2", ltmv[1][1], 9, NULL);
+    dst.appendReal("LTV1",   ltmv[2][0], 9, NULL);
+    dst.appendReal("LTV2",   ltmv[2][1], 9, NULL);
+  }
 
-  // LTMV,DTMV
   if (!isMosaic()) {
-    if (currentContext->fits->hasLTMV()) {
-      Matrix ltmv = currentContext->fits->physicalToRef * refToWidget *
-	Translate(-center) *
-	Translate(1,0) *
-	FlipY() *
-	Translate(center);
-
-      dst.appendReal("LTM1_1", ltmv[0][0], 9, NULL);
-      dst.appendReal("LTM1_2", ltmv[0][1], 9, NULL);
-      dst.appendReal("LTM2_1", ltmv[1][0], 9, NULL);
-      dst.appendReal("LTM2_2", ltmv[1][1], 9, NULL);
-      dst.appendReal("LTV1",   ltmv[2][0], 9, NULL);
-      dst.appendReal("LTV2",   ltmv[2][1], 9, NULL);
+    if (currentContext->fits->hasATMV()) {
+      Matrix dtmv = currentContext->fits->amplifierToRef * refToWidget * cc;
+      dst.appendReal("ATM1_1", dtmv[0][0], 9, NULL);
+      dst.appendReal("ATM1_2", dtmv[0][1], 9, NULL);
+      dst.appendReal("ATM2_1", dtmv[1][0], 9, NULL);
+      dst.appendReal("ATM2_2", dtmv[1][1], 9, NULL);
+      dst.appendReal("ATV1",   dtmv[2][0], 9, NULL);
+      dst.appendReal("ATV2",   dtmv[2][1], 9, NULL);
     }
-  }
-  else {
-    if (currentContext->fits->hasDTMV()) {
-      Matrix dtmv = currentContext->fits->detectorToRef * refToWidget *
-	Translate(-center) *
-	Translate(1,0) *
-	FlipY() *
-	Translate(center);
 
+    if (currentContext->fits->hasDTMV()) {
+      Matrix dtmv = currentContext->fits->detectorToRef * refToWidget * cc;
       dst.appendReal("DTM1_1", dtmv[0][0], 9, NULL);
       dst.appendReal("DTM1_2", dtmv[0][1], 9, NULL);
       dst.appendReal("DTM2_1", dtmv[1][0], 9, NULL);
@@ -437,124 +422,46 @@ void FrameBase::saveFitsResampleKeyword(OutFitsStream& str, FitsHead& dst)
   }
 
   // WCS
-  if (currentContext->fits->hasWCS(Coord::WCS)) {
-    if (src->find("RADESYS"))
-      dst.appendString("RADESYS", src->getString("RADESYS"), NULL);
-    if (src->find("EQUINOX"))
-      dst.appendReal("EQUINOX", src->getReal("EQUINOX",2000), 9, NULL);
-    if (src->find("EPOCH"))
-      dst.appendReal("EPOCH", src->getReal("EPOCH",2000), 9, NULL);
-    if (src->find("MJD-OBS"))
-      dst.appendReal("MJD-OBS", src->getReal("MJD-OBS",51544), 9, NULL);
-    if (src->find("CTYPE1"))
-      dst.appendString("CTYPE1", src->getString("CTYPE1"), NULL);
-    if (src->find("CTYPE2"))
-      dst.appendString("CTYPE2", src->getString("CTYPE2"), NULL);
-    if (src->find("CRVAL1"))
-      dst.appendReal("CRVAL1", src->getReal("CRVAL1",0), 9, NULL);
-    if (src->find("CRVAL2"))
-      dst.appendReal("CRVAL2", src->getReal("CRVAL2",0), 9, NULL);
-    if (src->find("CUNIT1"))
-      dst.appendString("CUNIT1", src->getString("CUNIT1"), NULL);
-    if (src->find("CUNIT2"))
-      dst.appendString("CUNIT2", src->getString("CUNIT2"), NULL);
-
-    // crpix
-    if (src->find("CRPIX1") || src->find("CRPIX2")) {
-      double crpix1 = src->getReal("CRPIX1",0);
-      double crpix2 = src->getReal("CRPIX2",0);
-
-      Vector crpix = Vector(crpix1,crpix2) * 
-	currentContext->fits->imageToWidget *
-	Translate(-center) *
-	Translate(1,0) *
-	FlipY() *
-	Translate(center);
-
-      dst.appendReal("CRPIX1", crpix[0], 9, NULL);
-      dst.appendReal("CRPIX2", crpix[1], 9, NULL);
-    }
-
-    // cd 
-    if (src->find("CD1_1") || src->find("CD1_2") || 
-	src->find("CD2_1") || src->find("CD2_2")) {
-      // cd keywords
-      double cd11 = src->getReal("CD1_1",0);
-      double cd12 = src->getReal("CD1_2",0);
-      double cd21 = src->getReal("CD2_1",0);
-      double cd22 = src->getReal("CD2_2",0);
-
-      Matrix cd = Matrix(cd11, cd12, cd21, cd22,0,0) *
-	currentContext->fits->imageToRef * refToUser *
-	wcsOrientationMatrix *
-	Rotate(wcsRotation) *
-	orientationMatrix *
-	Scale(zoom_.invert()) *
-	Rotate(rotation) *
-	Translate(center) *
-	Translate(-center) *
-	Translate(1,0) *
-	FlipY() * 
-	Translate(center);
-
-      dst.appendReal("CD1_1", cd.matrix(0,0), 9, NULL);
-      dst.appendReal("CD1_2", cd.matrix(0,1), 9, NULL);
-      dst.appendReal("CD2_1", cd.matrix(1,0), 9, NULL);
-      dst.appendReal("CD2_2", cd.matrix(1,1), 9, NULL);
-    }
-    else if (src->find("PC1_1") || src->find("PC1_2") || 
-	     src->find("PC2_1") || src->find("PC2_2")) {
-      // pc keywords
-      double pc11 = src->getReal("PC1_1",1);
-      double pc12 = src->getReal("PC1_2",0);
-      double pc21 = src->getReal("PC2_1",0);
-      double pc22 = src->getReal("PC2_2",1);
-      double cdelt1 = src->getReal("CDELT1",1);
-      double cdelt2 = src->getReal("CDELT2",1);
-
-      Matrix cd = Scale(cdelt1,cdelt2) * Matrix(pc11,pc12,pc21,pc22,0,0) *
-	currentContext->fits->imageToRef * refToUser *
-	wcsOrientationMatrix *
-	Rotate(wcsRotation) *
-	orientationMatrix *
-	Scale(zoom_.invert()) *
-	Rotate(rotation) *
-	Translate(center) *
-	Translate(-center) *
-	Translate(1,0) *
-	FlipY() * 
-	Translate(center);
-
-      dst.appendReal("CD1_1", cd.matrix(0,0), 9, NULL);
-      dst.appendReal("CD1_2", cd.matrix(0,1), 9, NULL);
-      dst.appendReal("CD2_1", cd.matrix(1,0), 9, NULL);
-      dst.appendReal("CD2_2", cd.matrix(1,1), 9, NULL);
-    }
-    else if (src->find("CDELT1") || src->find("CDELT2")) {
-      // crota2
-      double cdelt1 = src->getReal("CDELT1",1);
-      double cdelt2 = src->getReal("CDELT2",1);
-      double crot2 = src->getReal("CROT2",0);
-
-      Matrix cd = Scale(cdelt1,cdelt2) * Rotate(crot2) *
-	currentContext->fits->imageToRef * refToUser *
-	wcsOrientationMatrix *
-	Rotate(wcsRotation) *
-	orientationMatrix *
-	Scale(zoom_.invert()) *
-	Rotate(rotation) *
-	Translate(center) *
-	Translate(-center) *
-	Translate(1,0) *
-	FlipY() * 
-	Translate(center);
-
-      dst.appendReal("CD1_1", cd.matrix(0,0), 9, NULL);
-      dst.appendReal("CD1_2", cd.matrix(0,1), 9, NULL);
-      dst.appendReal("CD2_1", cd.matrix(1,0), 9, NULL);
-      dst.appendReal("CD2_2", cd.matrix(1,1), 9, NULL);
-    }
+  astClearStatus; // just to make sure
+  astBegin; // start memory management
+
+  // create channel and set encoding
+  AstFitsChan* chan = astFitsChan(NULL, NULL, "");
+  const char* fitswcs = "FITS-WCS";
+  const char* encode = currentContext->fits->encoding_;
+  if (!encode || !*encode)
+    encode = fitswcs;
+  astSet (chan, "Card=1, Encoding=%s", encode);
+
+  // 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);
+
+  // write to channel
+  if (!astWrite(chan, ast)) {
+    // try again
+    encode = fitswcs;
+    astSet (chan, "Card=1, Encoding=%s", encode);
+    astWrite(chan, ast);
   }
+    
+  // dump cards from channel
+  astClear(chan, "Card");
+  char card[81];
+  while (astFindFits(chan, "%f", card, 1))
+    dst.cardins(card,NULL);
+
+  astEnd; // now, clean up memory
 }
 
 void FrameBase::saveFitsResampleFits(OutFitsStream& str)