update AST WCS

3 files changed, 134 insertions, 6 deletions

diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C
index 49eb5df..ec1387e 100644
--- a/tksao/frame/fitsimage.C
+++ b/tksao/frame/fitsimage.C
@@ -105,6 +105,9 @@ FitsImage::FitsImage(Context* cx, Tcl_Interp* pp)
   imageToData3d = Translate3d(-.5, -.5, -.5);
   dataToImage3d = Translate3d( .5,  .5,  .5);
 
+  imageToRef3d = imageToData3d;
+  refToImage3d = dataToImage3d;
+
   manageWCS_ =1;
   wcsx_ =NULL;
 #ifndef NEWWCS
@@ -3119,7 +3122,7 @@ Vector FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys,
   return Vector();
 }
 #else
-Vector FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, 
+Vector FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys,
 			  Coord::SkyFrame sky)
 {
   astClearStatus; // just to make sure
@@ -3140,10 +3143,32 @@ Vector FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys,
 
   return Vector();
 }
+
+Vector3d FitsImage::pix2wcs(Vector3d in, Coord::CoordSystem sys,
+			    Coord::SkyFrame sky)
+{
+  astClearStatus; // just to make sure
+
+  if (!hasWCS(sys))
+    return Vector();
+
+  setWCSSystem(newast_,sys);
+  setWCSSkyFrame(newast_,sky);
+
+  Vector out = wcsTran(newast_, in, 1);
+  if (astOK && checkWCS(out)) {
+    if (wcsIsASkyFrame(newast_))
+      return out.radToDeg();
+    else
+      return out;
+  }
+
+  return Vector3d();
+}
 #endif
 
 #ifndef NEWWCS
-char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, 
+char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys,
 			 Coord::SkyFrame sky, Coord::SkyFormat format,
 			 char* lbuf)
 {
@@ -3204,7 +3229,7 @@ char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys,
   return lbuf;
 }
 #else
-char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, 
+char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys,
 			 Coord::SkyFrame sky, Coord::SkyFormat format,
 			 char* lbuf)
 {
@@ -3266,7 +3291,7 @@ char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys,
 #endif
 
 #ifndef NEWWCS
-Vector FitsImage::wcs2pix(Vector in, Coord::CoordSystem sys, 
+Vector FitsImage::wcs2pix(Vector in, Coord::CoordSystem sys,
 			  Coord::SkyFrame sky)
 {
   astClearStatus;
@@ -3290,7 +3315,7 @@ Vector FitsImage::wcs2pix(Vector in, Coord::CoordSystem sys,
   return Vector();
 }
 #else
-Vector FitsImage::wcs2pix(Vector in, Coord::CoordSystem sys, 
+Vector FitsImage::wcs2pix(Vector in, Coord::CoordSystem sys,
 			  Coord::SkyFrame sky)
 {
   astClearStatus; // just to make sure
@@ -3310,6 +3335,26 @@ Vector FitsImage::wcs2pix(Vector in, Coord::CoordSystem sys,
   maperr =1;
   return Vector();
 }
+
+Vector3d FitsImage::wcs2pix(Vector3d in, Coord::CoordSystem sys,
+			    Coord::SkyFrame sky)
+{
+  astClearStatus; // just to make sure
+
+  if (hasWCS(sys)) {
+    setWCSSystem(newast_,sys);
+    setWCSSkyFrame(newast_,sky);
+
+    if (wcsIsASkyFrame(newast_))
+      in *= M_PI/180.;
+
+    Vector out = wcsTran(newast_, in, 0);
+    if (astOK && checkWCS(out))
+      return out;
+  }
+
+  return Vector3d();
+}
 #endif
 
 #ifndef NEWWCS
@@ -4097,6 +4142,42 @@ void FitsImage::wcsTran(AstFrameSet* ast, int npoint,
     break;
   }
 }
+
+Vector3d FitsImage::wcsTran(AstFrameSet* ast, Vector3d& in, int forward)
+{
+  int naxes = astGetI(ast,"Naxes");
+  switch (naxes) {
+  case 1:
+  case 2:
+    // error
+    break;
+  case 3:
+    {
+      double pin[3];
+      double pout[3];
+      pin[0] = in[0];
+      pin[1] = in[1];
+      pin[2] = in[3];
+      astTranN(ast, 1, 3, 1, pin, forward, 3, 1, pout);
+      return Vector3d(pout[0],pout[1],pout[2]);
+    }
+    break;
+  case 4:
+    {
+      double pin[4];
+      double pout[4];
+      pin[0] = in[0];
+      pin[1] = in[1];
+      pin[2] = in[3];
+      pin[3] = forward ? context_->slice(3) : 0;
+      astTranN(ast, 1, 4, 1, pin, forward, 4, 1, pout);
+      return Vector3d(pout[0],pout[1],pout[2]);
+    }
+    break;
+  }
+  return Vector3d();
+}
+
 #endif
 
 #ifndef NEWWCS
diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h
index 38872f6..44ea871 100644
--- a/tksao/frame/fitsimage.h
+++ b/tksao/frame/fitsimage.h
@@ -209,6 +209,8 @@ class FitsImage {
 
   Matrix3d dataToImage3d;
   Matrix3d imageToData3d;
+  Matrix3d refToImage3d;
+  Matrix3d imageToRef3d;
 
   Matrix3d dataToRef3d;
   Matrix3d refToData3d;
@@ -366,6 +368,10 @@ class FitsImage {
 
   Vector wcs2pix(Vector, Coord::CoordSystem, Coord::SkyFrame);
 
+#ifdef NEWWCS
+  Vector3d pix2wcs(Vector3d, Coord::CoordSystem, Coord::SkyFrame);
+  Vector3d wcs2pix(Vector3d, Coord::CoordSystem, Coord::SkyFrame);
+#endif
   double pix2wcsx(double, Coord::CoordSystem, int);
   double wcs2pixx(double, Coord::CoordSystem, int);
 
@@ -403,6 +409,7 @@ class FitsImage {
   void wcsTran(AstFrameSet*, int, Vector*, int, Vector*);
   double wcsDistance(AstFrameSet*, Vector, Vector);
 #ifdef NEWWCS
+  Vector3d wcsTran(AstFrameSet*, Vector3d&, int);
   double wcsAngle(AstFrameSet*, Vector&, Vector&, Vector&);
   double wcsAxAngle(AstFrameSet*, Vector&, Vector&);
 #endif
@@ -444,6 +451,10 @@ class FitsImage {
   Vector mapFromRef(const Vector&, Coord::CoordSystem, Coord::SkyFrame =Coord::FK5);
   void mapFromRef(const Vector&, Coord::CoordSystem, Coord::SkyFrame, Coord::SkyFormat, char*);
   Vector mapToRef(const Vector&, Coord::CoordSystem, Coord::SkyFrame =Coord::FK5);
+#ifdef NEWWCS
+  Vector3d mapFromRef(const Vector3d&, Coord::CoordSystem, Coord::SkyFrame =Coord::FK5);
+  Vector3d mapToRef(const Vector3d&, Coord::CoordSystem, Coord::SkyFrame =Coord::FK5);
+#endif
   double mapFromRef3axis(double, Coord::CoordSystem, int);
   double mapToRef3axis(double, Coord::CoordSystem, int);
   double mapLenFromRef(double, Coord::CoordSystem, Coord::DistFormat =Coord::DEGREE);
diff --git a/tksao/frame/fitsmap.C b/tksao/frame/fitsmap.C
index 60a453a..618aad7 100644
--- a/tksao/frame/fitsmap.C
+++ b/tksao/frame/fitsmap.C
@@ -20,7 +20,7 @@ Vector FitsImage::mapFromRef(const Vector& vv, Coord::CoordSystem out,
     return vv * refToDetector;
   default:
     if (hasWCS(out))
-      return pix2wcs(vv * refToImage, out, sky);
+      return pix2wcs(vv * refToImage3d, out, sky);
   }
 
   return Vector();
@@ -60,6 +60,42 @@ Vector FitsImage::mapToRef(const Vector& vv, Coord::CoordSystem in,
   return Vector();
 }
 
+#ifdef NEWWCS
+Vector3d FitsImage::mapFromRef(const Vector3d& vv, Coord::CoordSystem out,
+			       Coord::SkyFrame sky)
+{
+  switch (out) {
+  case Coord::IMAGE:
+  case Coord::PHYSICAL:
+  case Coord::AMPLIFIER:
+  case Coord::DETECTOR:
+    return vv * refToImage3d;
+  default:
+    if (hasWCS(out))
+      return pix2wcs(vv * refToImage3d, out, sky);
+  }
+
+  return Vector3d();
+}
+
+Vector3d FitsImage::mapToRef(const Vector3d& vv, Coord::CoordSystem in,
+			     Coord::SkyFrame sky)
+{
+  switch (in) {
+  case Coord::IMAGE:
+  case Coord::PHYSICAL:
+  case Coord::AMPLIFIER:
+  case Coord::DETECTOR:
+    return vv * imageToRef3d;
+  default:
+    if (hasWCS(in))
+      return wcs2pix(vv, in, sky) * imageToRef3d;
+  }
+
+  return Vector3d();
+}
+#endif
+
 void FitsImage::listFromRef(ostream& str, const Vector& vv,
 			    Coord::CoordSystem sys, Coord::SkyFrame sky, 
 			    Coord::SkyFormat format)