fix block WCS

6 files changed, 36 insertions, 82 deletions

diff --git a/tksao/fitsy++/block.C b/tksao/fitsy++/block.C
index cea9293..2d0c477 100644
--- a/tksao/fitsy++/block.C
+++ b/tksao/fitsy++/block.C
@@ -81,8 +81,7 @@ void FitsBlock::initHeader(FitsFile* fits, Vector& block)
   initKeySEC("DATASEC",block);
   initKeySEC("TRIMSEC",block);
 
-  // WCS
-  initWCS(block);
+  // WCS blocking done later in FitsImage::resetWCS()
 
   // clear cards
   if (head_->find("BZERO"))
@@ -193,77 +192,3 @@ void FitsBlock::initLTMV(Vector& block)
   else
     head_->appendReal("LTM2_2", ltm22/block[1], 8, "");
 }
-
-void FitsBlock::initWCS(Vector& block)
-{
-  // check for WCS
-  if (!head_->find("CRPIX1") && !head_->find("CRPIX2"))
-    return;
-
-  double crpix1 = head_->getReal("CRPIX1", 0);
-  double crpix2 = head_->getReal("CRPIX2", 0);
-  Matrix pp = Matrix(1,0,0,1,crpix1,crpix2) *
-    Translate(-.5,-.5) *
-    Scale(1/block[0],1/block[1]) *
-    Translate(.5,.5);
-
-  if (head_->find("CRPIX1"))
-    head_->setReal("CRPIX1", pp.matrix(2,0), 8, "");
-  else
-    head_->appendReal("CRPIX1", pp.matrix(2,0), 8, "");
-
-  if (head_->find("CRPIX2"))
-    head_->setReal("CRPIX2", pp.matrix(2,1), 8, "");
-  else
-    head_->appendReal("CRPIX2", pp.matrix(2,1), 8, "");
-
-  // CD
-  if (head_->find("CD1_1") || head_->find("CD1_2") ||
-      head_->find("CD2_1") || head_->find("CD2_2")) {
-    double cd11 = head_->getReal("CD1_1", 1);
-    double cd12 = head_->getReal("CD1_2", 0);
-    double cd21 = head_->getReal("CD2_1", 0);
-    double cd22 = head_->getReal("CD2_2", 1);
-    Matrix mm = Matrix(cd11,cd12,cd21,cd22,0,0) *
-      Scale(block[0],block[1]);
-
-    if (head_->find("CD1_1"))
-      head_->setReal("CD1_1", mm.matrix(0,0), 8, "");
-    else
-      head_->appendReal("CD1_1", mm.matrix(0,0), 8, "");
-
-    if (head_->find("CD1_2"))
-      head_->setReal("CD1_2", mm.matrix(0,1), 8, "");
-    else
-      head_->appendReal("CD1_2", mm.matrix(0,1), 8, "");
-
-    if (head_->find("CD2_1"))
-      head_->setReal("CD2_1", mm.matrix(1,0), 8, "");
-    else
-      head_->appendReal("CD2_1", mm.matrix(1,0), 8, "");
-
-    if (head_->find("CD2_2"))
-      head_->setReal("CD2_2", mm.matrix(1,1), 8, "");
-    else
-      head_->appendReal("CD2_2", mm.matrix(1,1), 8, "");
-
-  }
-  // CDELT
-  else if (head_->find("CDELT1") || head_->find("CDELT2")) {
-    double cdelt1 = head_->getReal("CDELT1", 1);
-    double cdelt2 = head_->getReal("CDELT2", 0);
-    Matrix mm = Matrix(cdelt1,0,0,cdelt2,0,0) *
-      Scale(block[0],block[1]);
-
-    if (head_->find("CDELT1"))
-      head_->setReal("CDELT1", mm.matrix(0,0), 8, "");
-    else
-      head_->appendReal("CDELT1", mm.matrix(0,0), 8, "");
-
-    if (head_->find("CDELT2"))
-      head_->setReal("CDELT2", mm.matrix(1,1), 8, "");
-    else
-      head_->appendReal("CDELT2", mm.matrix(1,1), 8, "");
-  }
-}
-
diff --git a/tksao/fitsy++/block.h b/tksao/fitsy++/block.h
index 6da9e3f..5a5d0ef 100644
--- a/tksao/fitsy++/block.h
+++ b/tksao/fitsy++/block.h
@@ -20,7 +20,6 @@ class FitsBlock : public FitsFile {
   void initCCDSUM(Vector&);
   void initKeySEC(const char*, Vector&);
   void initLTMV(Vector&);
-  void initWCS(Vector&);
 
  public:
   FitsBlock(FitsFile*, Vector&);
diff --git a/tksao/frame/context.C b/tksao/frame/context.C
index e900df4..05dff3d 100644
--- a/tksao/frame/context.C
+++ b/tksao/frame/context.C
@@ -1249,7 +1249,7 @@ int Context::loadMosaicWFPC2(Base::MemType which, const char* fn,
 
   // Its legal, save it
   bfits_ = img;
-  loadInit(1, Base::WFPC2,Coord::WCS);
+  loadInit(1, Base::WFPC2, Coord::WCS);
 
   // remember in case of compress
   Base::MemType sav = which;
@@ -1418,7 +1418,7 @@ int Context::loadMosaicWFPC2(Base::MemType which, const char* fn,
 	   << ends;
 
       // fix fitsimage params
-      ptr->wfpc2WCS(istr);
+      ptr->wfpc2WCS(bfits_->head(), istr);
 
       Matrix mm = parent_->calcAlignWCS(bfits_, ptr, Coord::WCS,
 					Coord::WCS, Coord::FK5);
diff --git a/tksao/frame/fitsblock.C b/tksao/frame/fitsblock.C
index 49629a8..bb0a844 100644
--- a/tksao/frame/fitsblock.C
+++ b/tksao/frame/fitsblock.C
@@ -5,6 +5,7 @@
 #include "fitsimage.h"
 #include "block.h"
 #include "context.h"
+#include "wcsast.h"
 
 void* blockproc(void* tt);
 
@@ -187,4 +188,3 @@ void* blockproc(void* tt)
 
   return NULL;
 }
-
diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C
index aaf972e..889befd 100644
--- a/tksao/frame/fitsimage.C
+++ b/tksao/frame/fitsimage.C
@@ -889,10 +889,22 @@ void FitsImageIIS::iisWCS(const Matrix& mm, const Vector& z, int zt)
 
 // FitsImage
 
-void FitsImage::wfpc2WCS(istream& str)
+void FitsImage::wfpc2WCS(FitsHead* pp, istream& str)
 {
   FitsHead* hh = parseWCS(str);
 
+  // EQUINOX
+  if (pp->find("EQUINOX")) {
+    char* equ = pp->getString("EQUINOX");
+    hh->appendString("EQUINOX", equ, NULL);
+  }
+
+  // DATE-OBS
+  if (pp->find("DATE-OBS")) {
+    char* equ = pp->getString("DATE-OBS");
+    hh->appendString("DATE-OBS", equ, NULL);
+  }
+
   // Process OBJECT keyword
   if (objectKeyword_)
     delete [] objectKeyword_;
@@ -1178,6 +1190,24 @@ void FitsImage::resetWCS()
     initWCS(wfpc2Header_);
   else
     initWCS(image_->head());
+
+
+  // apply block factor
+  if (ast_) {
+    Vector block = context_->blockFactor();
+    if (block[0] != 1 && block[1] != 1) {
+      astClearStatus; // just to make sure
+      astBegin; // start memory management
+
+      Vector ll(.5,.5);
+      Vector ur(1.5,1.5);
+      Vector rr = ur*Translate(-.5,-.5)*Scale(block)*Translate(.5,.5);
+      AstWinMap* winmap = astWinMap(2, ll.vv(), rr.vv(), ll.vv(), ur.vv(), "");
+      astRemapFrame(ast_, AST__BASE, winmap);
+
+      astEnd;
+    }
+  }
 }
 
 void FitsImage::initWCS0(const Vector& pix)
diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h
index aad99b9..1acc5bf 100644
--- a/tksao/frame/fitsimage.h
+++ b/tksao/frame/fitsimage.h
@@ -408,7 +408,7 @@ class FitsImage {
   Vector3d wcs2pix(const Vector3d&, Coord::CoordSystem, Coord::SkyFrame);
   int wcsInv() {return wcsInv_;}
 
-  void wfpc2WCS(istream&);
+  void wfpc2WCS(FitsHead*, istream&);
   void appendWCS(istream&);
   void listWCS(ostream&, Coord::CoordSystem);
   void resetWCS();