update AST WCS

3 files changed, 105 insertions, 16 deletions

diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C
index 81ff0a7..0c7eba5 100644
--- a/tksao/frame/fitsimage.C
+++ b/tksao/frame/fitsimage.C
@@ -3346,7 +3346,9 @@ char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys,
 
   return lbuf;
 }
+
 #else
+
 char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, 
 			 Coord::SkyFrame sky, Coord::SkyFormat format,
 			 char* lbuf)
@@ -3925,8 +3927,8 @@ void FitsImage::astinit(int ss, FitsHead* hd, FitsHead* prim)
 
   setAstWCSSystem(ast_[ss], (Coord::CoordSystem)(ss+Coord::WCS));
 
-  astClearStatus; // just to make sure
-  astBegin; // start memory management
+  //  astClearStatus; // just to make sure
+  //  astBegin; // start memory management
 
   int naxes = astGetI(ast_[ss],"Naxes");
   switch (naxes) {
@@ -3942,6 +3944,7 @@ void FitsImage::astinit(int ss, FitsHead* hd, FitsHead* prim)
   case 3:
   case 4:
     {
+      if (0) {
       AstFrameSet* ast = ast_[ss];
 
       int pickc[2] = {1,2};
@@ -3959,11 +3962,12 @@ void FitsImage::astinit(int ss, FitsHead* hd, FitsHead* prim)
       int idata =  astGetI(ast, "Current");
       astSetI(ast, "Current", isky);
       astSetI(ast, "Base", idata);
+      }
     }
     break;
   }
 
-  astEnd; // now, clean up memory
+  //  astEnd; // now, clean up memory
 #endif
 
   // set default skyframe
@@ -3982,9 +3986,6 @@ void FitsImage::astinit(FitsHead* hd, FitsHead* prim)
   if (!newast_)
     return;
 
-  astClearStatus; // just to make sure
-  astBegin; // start memory management
-
   int naxes = astGetI(newast_,"Naxes");
   switch (naxes) {
   case 1:
@@ -3999,6 +4000,10 @@ void FitsImage::astinit(FitsHead* hd, FitsHead* prim)
   case 3:
   case 4:
     {
+      if (0) {
+      astClearStatus; // just to make sure
+      astBegin; // start memory management
+
       AstFrameSet* ast = newast_;
 
       int pickc[2] = {1,2};
@@ -4016,13 +4021,15 @@ void FitsImage::astinit(FitsHead* hd, FitsHead* prim)
       int idata =  astGetI(ast, "Current");
       astSetI(ast, "Current", isky);
       astSetI(ast, "Base", idata);
+
+      astEnd; // now, clean up memory
+      }
     }
     break;
   }
   
-  astEnd; // now, clean up memory
-
-  setAstWCSSkyFrame(newast_,Coord::FK5);
+  if (astWCSIsASkyFrame(newast_))
+    setAstWCSSkyFrame(newast_,Coord::FK5);
 }
 #endif
 
@@ -4205,7 +4212,40 @@ void FitsImage::astWCSTran(AstFrameSet* ast, int npoint,
 			   int forward,
 			   double* xout, double* yout)
 {
-  astTran2(ast, npoint, xin, yin, forward, xout, yout);
+  //  astTran2(ast, npoint, xin, yin, forward, xout, yout);
+
+  int naxes = astGetI(ast,"Naxes");
+  switch (naxes) {
+  case 1:
+    // error
+    break;
+  case 2:
+    astTran2(ast, npoint, xin, yin, forward, xout, yout);
+    break;
+  case 3:
+    {
+      double* ptr_in[3];
+      ptr_in[0] = (double*)xin;
+      ptr_in[1] = (double*)yin;
+      ptr_in[2] = new double[npoint];
+      for (int kk=0; kk<npoint; kk++)
+	ptr_in[2][kk] = 1;
+      
+
+      double* ptr_out[3];
+      ptr_out[0] = (double*)xout;
+      ptr_out[1] = (double*)yout;
+      ptr_out[2] = new double[npoint];
+
+      astTranP(ast, npoint, 3, (const double**)ptr_in, forward, 3, ptr_out);
+
+      if (ptr_in[2])
+	delete [] ptr_in[2];
+      if (ptr_out[2])
+	delete [] ptr_out[2];
+    }
+    break;
+  }
 }
 
 #endif
diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h
index c14c193..3759285 100644
--- a/tksao/frame/fitsimage.h
+++ b/tksao/frame/fitsimage.h
@@ -409,9 +409,13 @@ class FitsImage {
 #endif
   void setAstWCSSkyFrame(AstFrameSet*, Coord::SkyFrame);
   void setAstWCSFormat(AstFrameSet*, int, const char*);
+#ifndef NEWWCS
   AstFrameSet* getAST(Coord::CoordSystem sys) 
     {return (ast_ && ast_[sys-Coord::WCS]) ? ast_[sys-Coord::WCS] : NULL;}
-
+#else
+  AstFrameSet* getAST(Coord::CoordSystem sys) {return newast_;}
+#endif
+  
   int hasWCS(Coord::CoordSystem);
   int hasWCSEqu(Coord::CoordSystem);
   int hasWCSCel(Coord::CoordSystem);
diff --git a/tksao/frame/grid2d.C b/tksao/frame/grid2d.C
index 6d5bf22..aea9320 100644
--- a/tksao/frame/grid2d.C
+++ b/tksao/frame/grid2d.C
@@ -68,9 +68,9 @@ int Grid2d::doit(RenderMode rm)
     break;
   default:
     {
+#ifndef NEWWCS
       AstFrameSet* ast = (AstFrameSet*)astCopy(fits->getAST(system_));
 
-#ifndef NEWWCS
       // set desired skyformat
       if (fits->astWCSIsASkyFrame(ast))
 	fits->setAstWCSSkyFrame(ast, sky_);
@@ -81,14 +81,59 @@ int Grid2d::doit(RenderMode rm)
       astInvert(ast);
       astAddFrame(frameSet,2,astUnitMap(2,""),ast);
       astSetI(frameSet,"current",astGetI(frameSet,"nframe"));
-
 #else
+      astClearStatus; // just to make sure
+      astBegin; // start memory management
+
+      AstFrameSet* ast = (AstFrameSet*)astCopy(fits->getAST(system_));
+
+      fits->setAstWCSSystem(ast, system_);
+      fits->setAstWCSSkyFrame(ast, sky_);
+      
+      int offset =0;
+      int naxes = astGetI(ast,"Naxes");
+      switch (naxes) {
+      case 1:
+      case 2:
+	break;
+      case 3:
+      case 4:
+	{
+	  int pick[2] = {1, 2};
+	  astInvert(ast);
+	  AstMapping* mapb =NULL;
+	  AstFrame* permb = (AstFrame*)astPickAxes(ast, 2, pick, &mapb);
+	  astAddFrame(ast, AST__CURRENT, mapb, permb);
+	  astInvert(ast);
+	  
+	  AstMapping** mapc =NULL;
+	  AstFrame* permc = (AstFrame*)astPickAxes(ast, 2, pick, &mapc);
+	  astAddFrame(ast, AST__CURRENT, mapc, permc);
+	  
+	  if (0) {
+	  int isky = astGetI(ast, "Current");
+	  int pickb[4] = {1, 2, 0, 0};
+	  AstMapping* mapb;
+	  AstFrame* foo = astFrame(2,"Domain=DATA");
+	  astPickAxes(foo, naxes, pickb, &mapb);
+	  astInvert(mapb);
+	  astAddFrame(ast, AST__BASE, mapb, foo);
+	  int idata =  astGetI(ast, "Current");
+	  astSetI(ast, "Current", isky);
+	  astSetI(ast, "Base", idata);
+	  offset =1;
+	  }
+	}
+	break;
+      }
+
       // add wcs to frameset
       // this will link frameset to wcs with unitMap
       astInvert(ast);
-      astAddFrame(frameSet,AST__CURRENT,astUnitMap(2,""),ast);
-      fits->setAstWCSSystem(frameSet,system_);
-      fits->setAstWCSSkyFrame(frameSet,sky_);
+      astAddFrame(frameSet, AST__CURRENT, astUnitMap(2,""), ast);
+      astSetI(frameSet,"Current",astGetI(frameSet,"nframe")-offset);
+
+      astEnd; // now, clean up memory
 #endif
     }
   }