update AST WCS

2 files changed, 59 insertions, 26 deletions

diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C
index 805f779..6cb6b72 100644
--- a/tksao/frame/fitsimage.C
+++ b/tksao/frame/fitsimage.C
@@ -2929,32 +2929,19 @@ double FitsImage::getWCSPixelSize(Coord::CoordSystem sys)
   astClearStatus; // just to make sure
   setWCSSystem(newast_, sys);
   
-  Vector cc = center();
-  double xx[3], wxx[3];
-  xx[0] = cc[0];
-  xx[1] = cc[0];
-  xx[2] = cc[0];
-  double yy[3], wyy[3];
-  yy[0] = cc[1];
-  yy[1] = cc[1]+1;
-  yy[2] = cc[1]+1;
-  wcsTran(newast_,3,xx,yy,1,wxx,wyy);
-
-  double pt0[2];
-  pt0[0] = wxx[0];
-  pt0[1] = wyy[0];
-  double pt1[2];
-  pt1[0] = wxx[1];
-  pt1[1] = wyy[1];
-  double pt2[2];
-  pt2[0] = wxx[2];
-  pt2[1] = wyy[2];
-  double out = (wcsDistance(newast_,pt0,pt1)+wcsDistance(newast_,pt0,pt2))/2.;
+  Vector in[3];
+  Vector out[3];
+  in[0] = center();
+  in[1] = center()+Vector(1,0);
+  in[2] = center()+Vector(0,1);
+  wcsTran(newast_, 3, in, 1, out);
+  double dd = (wcsDistance(newast_,out[0],out[1]) +
+	       wcsDistance(newast_,out[0],out[2]))/2.;
 
   if (wcsIsASkyFrame(newast_))
-    return radToDeg(out);
+    return radToDeg(dd);
   else
-    return out;
+    return dd;
 }
 
 double FitsImage::getWCSPixelArea(Coord::CoordSystem sys)
@@ -4117,8 +4104,8 @@ void FitsImage::wcsTran(AstFrameSet* ast, int npoint,
 	delete [] yin;
       if (xout)
 	delete [] xout;
-      if (xout)
-	delete [] xout;
+      if (yout)
+	delete [] yout;
     }
     break;
   case 3:
@@ -4323,6 +4310,51 @@ double FitsImage::wcsDistance(AstFrameSet* ast, double* point1, double* point2)
   return 0;
 }
 
+double FitsImage::wcsDistance(AstFrameSet* ast, Vector vv1, Vector vv2)
+{
+  int naxes = astGetI(ast,"Naxes");
+  switch (naxes) {
+  case 1:
+    // error
+    break;
+  case 2:
+    return astDistance(ast, vv1.v, vv2.v);
+    break;
+  case 3:
+    {
+      double ptr1[3];
+      ptr1[0] = vv1[0];
+      ptr1[1] = vv1[1];
+      ptr1[2] = 0;
+      double ptr2[3];
+      ptr2[0] = vv2[0];
+      ptr2[1] = vv2[1];
+      ptr2[2] = 0;
+
+      return astDistance(ast, ptr1, ptr2);
+    }
+    break;
+  case 4:
+    {
+      double ptr1[4];
+      ptr1[0] = vv1[0];
+      ptr1[1] = vv1[1];
+      ptr1[2] = 0;
+      ptr1[3] = 0;
+      double ptr2[4];
+      ptr2[0] = vv2[0];
+      ptr2[1] = vv2[1];
+      ptr2[2] = 0;
+      ptr2[3] = 0;
+
+      return astDistance(ast, ptr1, ptr2);
+    }
+    break;
+  }
+
+  return 0;
+}
+
 #endif
 
 AstFrameSet* FitsImage::fits2ast(FitsHead* hd) 
diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h
index e949184..22624fb 100644
--- a/tksao/frame/fitsimage.h
+++ b/tksao/frame/fitsimage.h
@@ -399,12 +399,13 @@ class FitsImage {
 #endif
   
   int wcsIsASkyFrame(AstFrameSet*);
-  double wcsDistance(AstFrameSet*, double*, double*);
   Vector wcsTran(AstFrameSet*, Vector, int);
   void wcsTran(AstFrameSet*, int, Vector*, int, Vector*);
+  double wcsDistance(AstFrameSet*, double*, double*);
 #ifdef NEWWCS
   void wcsTran(AstFrameSet*, int, const double*, const double*, int,
 	       double*, double*);
+  double wcsDistance(AstFrameSet*, Vector, Vector);
 #endif
   
 #ifdef NEWWCS