new AST WCS

3 files changed, 103 insertions, 82 deletions

diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C
index eed9691..12ae087 100644
--- a/tksao/frame/fitsimage.C
+++ b/tksao/frame/fitsimage.C
@@ -2898,7 +2898,6 @@ double FitsImage::getWCSPixelSize(Coord::CoordSystem sys)
     return 0;
 
   astClearStatus; // just to make sure
-  astBegin; // start memory management
 
   Vector cc = center();
   double xx[2], wxx[2];
@@ -2918,35 +2917,46 @@ double FitsImage::getWCSPixelSize(Coord::CoordSystem sys)
   double out = astDistance(ast_[ss],pt0,pt1);
 
   if (astIsASkyFrame(astGetFrame(ast_[ss], AST__CURRENT)))
-    out = radToDeg(out);
-
-  astEnd; // now, clean up memory
-
-  return out;
+    return radToDeg(out);
+  else
+    return out;
 }
 
-Vector FitsImage::getWCScdelt(Coord::CoordSystem sys)
+double FitsImage::getWCSPixelArea(Coord::CoordSystem sys)
 {
   int ss = sys-Coord::WCS;
   if (!(ss>=0 && ast_ && ast_[ss]))
-    return Vector();
+    return 0;
+
+  astClearStatus; // just to make sure
 
   Vector cc = center();
-  Vector wcc;
-  astTran2(ast_[ss], 1, cc.v, cc.v+1, 1, wcc.v, wcc.v+1);
-  Vector oo = cc+Vector(1,1);
-  Vector woo;
-  astTran2(ast_[ss], 1, oo.v, oo.v+1, 1, woo.v, woo.v+1);
-
-  double dd = astDistance(ast_[ss], wcc.v, woo.v);
-  if (dd != AST__BAD) {
-    if (astIsASkyFrame(astGetFrame(ast_[ss], AST__CURRENT)))
-      dd *= 180./M_PI;
-    dd /=sqrt(2);
-    return Vector(dd,dd);
-  }
+  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;
+  astTran2(ast_[ss],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 ll = astDistance(ast_[ss],pt0,pt1);
+  double mm = astDistance(ast_[ss],pt0,pt2);
+
+  if (astIsASkyFrame(astGetFrame(ast_[ss], AST__CURRENT)))
+    return radToDeg(ll)*radToDeg(mm);
   else
-    return Vector();
+    return ll*mm;
 }
 #endif
 
@@ -2990,7 +3000,6 @@ Coord::Orientation FitsImage::getWCSOrientation(Coord::CoordSystem sys,
     return Coord::NORMAL;
 
   astClearStatus; // just to make sure
-  astBegin; // start memory management
 
   if (astIsASkyFrame(astGetFrame(ast_[ss], AST__CURRENT)))
     setAstSkyFrame(ast_[ss],sky);
@@ -3021,7 +3030,6 @@ Coord::Orientation FitsImage::getWCSOrientation(Coord::CoordSystem sys,
     else
       rr = ang<=0 ? Coord::NORMAL : Coord::XX;
   }
-  astEnd; // now, clean up memory
 
   return rr;
 }
@@ -3058,7 +3066,6 @@ double FitsImage::getWCSRotation(Coord::CoordSystem sys, Coord::SkyFrame sky)
     return 0;
 
   astClearStatus; // just to make sure
-  astBegin; // start memory management
 
   if (astIsASkyFrame(astGetFrame(ast_[ss], AST__CURRENT)))
     setAstSkyFrame(ast_[ss],sky);
@@ -3077,7 +3084,6 @@ double FitsImage::getWCSRotation(Coord::CoordSystem sys, Coord::SkyFrame sky)
   bb[0]= wx[1];
   bb[1]= wy[1];
   double ang = astAxAngle(ast_[ss],aa,bb,2);
-  astEnd; // now, clean up memory
 
   if (!(isnan(ang)||isinf(ang)||(ang == -DBL_MAX)||(ang == DBL_MAX)))
     return getWCSOrientation(sys,sky) == Coord::NORMAL ? ang : -ang;
@@ -3303,24 +3309,20 @@ Vector* FitsImage::wcs2pix(Vector* in, int num, Coord::CoordSystem sys,
 
 double FitsImage::wcsdist(Vector a, Vector b, Coord::CoordSystem sys)
 {
-  astClearStatus;
+  int ss = sys-Coord::WCS;
+  if (!(ss>=0 && ast_ && ast_[ss]))
+    return 0;
 
-  int ii = sys-Coord::WCS;
-  double rr=0;
-  if (ii>=0 && ast_ && ast_[ii]) {
-    if (astIsASkyFrame(astGetFrame(ast_[ii], AST__CURRENT))) {
-      Vector aa = a*M_PI/180.;
-      Vector bb = b*M_PI/180.;
-      rr = astDistance(ast_[ii], aa.v, bb.v) *180./M_PI;
-    }
-    else
-      rr = astDistance(ast_[ii], a.v, b.v);
+  astClearStatus; // just to make sure
 
-    if (!astOK) {
-      maperr =1;
-      return 0;
-    }
+  double rr=0;
+  if (astIsASkyFrame(astGetFrame(ast_[ss], AST__CURRENT))) {
+    Vector aa = a*M_PI/180.;
+    Vector bb = b*M_PI/180.;
+    rr = astDistance(ast_[ss], aa.v, bb.v) *180./M_PI;
   }
+  else
+    rr = astDistance(ast_[ss], a.v, b.v);
 
   return rr;
 }
diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h
index bc1c52a..a96c68c 100644
--- a/tksao/frame/fitsimage.h
+++ b/tksao/frame/fitsimage.h
@@ -388,9 +388,11 @@ class FitsImage {
   {return (wcs_ && wcs_[sys-Coord::WCS]) ? wcs_[sys-Coord::WCS]->wcsname : NULL;}
   Coord::Orientation getWCSOrientation(Coord::CoordSystem, Coord::SkyFrame);
   double getWCSRotation(Coord::CoordSystem, Coord::SkyFrame);
+#ifndef NEWWCS
   Vector getWCScdelt(Coord::CoordSystem);
-#ifdef NEWWCS
+#else
   double getWCSPixelSize(Coord::CoordSystem);
+  double getWCSPixelArea(Coord::CoordSystem);
 #endif
   
 #ifdef NEWWCS
diff --git a/tksao/frame/frblt.C b/tksao/frame/frblt.C
index 62755c4..95d812b 100644
--- a/tksao/frame/frblt.C
+++ b/tksao/frame/frblt.C
@@ -374,27 +374,28 @@ int Base::markerAnalysisRadial(Marker* pp, double** x, double** y, double** e,
   if (ptr->hasWCS(sys)) {
 #ifndef NEWWCS
     Vector cdelt= ptr->getWCScdelt(sys);
-    if (ptr->hasWCSCel(sys)) {
-      unit =1;
-      xaxis = fabs(cdelt[0]*60*60);
-    }
-    else {
-      unit =2;
-      xaxis = fabs(cdelt[0]);
-    }
+    double ll = fabs(cdelt[0]);
 #else
-    xaxis = ptr->getWCSPixelSize(sys);
+    double ll = ptr->getWCSPixelSize(sys);
+#endif
+
     if (ptr->hasWCSCel(sys)) {
       unit =1;
-      xaxis *= 60*60;
+      xaxis = ll*60*60;
     }
     else {
       unit =2;
+      xaxis = ll;
     }
-#endif
   }
-
+    
+#ifndef NEWWCS
   Vector cdelt= ptr->getWCScdelt(sys);
+  double aa = fabs(cdelt[0]*cdelt[1]);
+#else
+  double aa = ptr->getWCSPixelArea(sys);
+#endif
+
   for (int kk=0; kk<num; kk++) {
     double err = sqrt(fabs(sum[kk]));
     double area =0;
@@ -407,11 +408,11 @@ int Base::markerAnalysisRadial(Marker* pp, double** x, double** y, double** e,
       break;
     case 1:
       // Cel WCS
-      area = fabs(cdelt[0]*cdelt[1]*60*60*60*60*cnt[kk]);
+      area = aa*60*60*60*60*cnt[kk];
       break;
     case 2:
       // Linear WCS
-      area = fabs(cdelt[0]*cdelt[1]*cnt[kk]);
+      area = aa*cnt[kk];
       break;
     }
 
@@ -436,7 +437,7 @@ int Base::markerAnalysisRadial(Marker* pp, double** x, double** y, double** e,
 // for panda regions
 int Base::markerAnalysisPanda(Marker* pp, double** x, double** y, double** e, 
 			       int num, Vector* annuli, 
-			       int aa, double* angles,
+			       int na, double* angles,
 			       BBox* bb, Coord::CoordSystem sys)
 
 {
@@ -449,10 +450,10 @@ int Base::markerAnalysisPanda(Marker* pp, double** x, double** y, double** e,
   int srcw = ptr->width();
   FitsBound* params = ptr->getDataParams(currentContext->secMode());
 
-  double sum[num][aa];
-  memset(sum,0,num*aa*sizeof(double));
-  int cnt[num][aa];
-  memset(cnt,0,num*aa*sizeof(int));
+  double sum[num][na];
+  memset(sum,0,num*na*sizeof(double));
+  int cnt[num][na];
+  memset(cnt,0,num*na*sizeof(int));
 
   for (int kk=0; kk<num; kk++) {
     // take the bbox and extend to lower/upper pixel boundaries
@@ -461,7 +462,7 @@ int Base::markerAnalysisPanda(Marker* pp, double** x, double** y, double** e,
 
     // main loop
     SETSIGBUS
-      for (int qq=0; qq<aa; qq++) {
+      for (int qq=0; qq<na; qq++) {
 
 	for (int jj=ll[1]; jj<ur[1]; jj++) {
 	  for (int ii=ll[0]; ii<ur[0]; ii++) {
@@ -487,36 +488,37 @@ int Base::markerAnalysisPanda(Marker* pp, double** x, double** y, double** e,
     CLEARSIGBUS
   }
 
-  *x = (double*)malloc(num*aa*sizeof(double));
-  *y = (double*)malloc(num*aa*sizeof(double));
-  *e = (double*)malloc(num*aa*sizeof(double));
+  *x = (double*)malloc(num*na*sizeof(double));
+  *y = (double*)malloc(num*na*sizeof(double));
+  *e = (double*)malloc(num*na*sizeof(double));
 
   int unit =0;
   double xaxis =1;
   if (ptr->hasWCS(sys)) {
 #ifndef NEWWCS
     Vector cdelt= ptr->getWCScdelt(sys);
-    if (ptr->hasWCSCel(sys)) {
-      unit =1;
-      xaxis = fabs(cdelt[0]*60*60);
-    }
-    else {
-      unit =2;
-      xaxis = fabs(cdelt[0]);
-    }
+    double ll = fabs(cdelt[0]);
 #else
-    xaxis = ptr->getWCSPixelSize(sys);
+    double ll = ptr->getWCSPixelSize(sys);
+#endif
+
     if (ptr->hasWCSCel(sys)) {
       unit =1;
-      xaxis *= 60*60;
+      xaxis = ll*60*60;
     }
     else {
       unit =2;
+      xaxis = ll;
     }
-#endif
   }
 
+#ifndef NEWWCS
   Vector cdelt= ptr->getWCScdelt(sys);
+  double aa = fabs(cdelt[0]*cdelt[1]);
+#else
+  double aa= ptr->getWCSPixelArea(sys);
+#endif
+
   for (int qq=0; qq<aa; qq++) {
     for (int kk=0; kk<num; kk++) {
       double err = sqrt(fabs(sum[kk][qq]));
@@ -530,11 +532,11 @@ int Base::markerAnalysisPanda(Marker* pp, double** x, double** y, double** e,
 	break;
       case 1:
 	// Cel WCS
-	area = fabs(cdelt[0]*cdelt[1]*60*60*60*60*cnt[kk][qq]);
+	area = aa*60*60*60*60*cnt[kk][qq];
 	break;
       case 2:
 	// Linear WCS
-	area = fabs(cdelt[0]*cdelt[1]*cnt[kk][qq]);
+	area = aa*cnt[kk][qq];
 	break;
       }
 
@@ -830,9 +832,14 @@ int Base::markerAnalysisStats1(Marker* pp,FitsImage* ptr, ostream& str,
     return 0;
   default: 
     {
+#ifndef NEWWCS
       Vector cdelt= ptr->getWCScdelt(sys);
+      double ll = fabs(cdelt[0]);
+#else
+      double ll = ptr->getWCSPixelSize(sys);
+#endif
       if (ptr->hasWCSCel(sys)) {
-	str << "1 pixel = "<< fabs(cdelt[0]*60*60) << " arcsec";
+	str << "1 pixel = "<< ll*60*60 << " arcsec";
 	str << endl << endl;
 	str << "reg\t" << "sum\t" << "error\t\t" 
 	    << "area\t\t" << "surf_bri\t\t" << "surf_err" << endl
@@ -843,7 +850,7 @@ int Base::markerAnalysisStats1(Marker* pp,FitsImage* ptr, ostream& str,
 	return 1;
       }
       else {
-	str << "1 pixel = "<< fabs(cdelt[0]);
+	str << "1 pixel = "<< ll;
 	str << endl << endl;
 	str << "reg\t" << "sum\t" << "error\t\t" 
 	    << "area\t\t" << "surf_bri\t\t" << "surf_err" << endl
@@ -871,15 +878,25 @@ void Base::markerAnalysisStats2(FitsImage* ptr, ostream& str,
   case 1:
     {
       // Cel WCS
+#ifndef NEWWCS
       Vector cdelt= ptr->getWCScdelt(sys);
-      area = fabs(cdelt[0]*cdelt[1]*60*60*60*60*cnt);
+      double aa = fabs(cdelt[0]*cdelt[1]);
+#else
+      double aa = ptr->getWCSPixelArea(sys);
+#endif
+      area = aa*60*60*60*60*cnt;
     }
     break;
   case 2:
     {
       // Linear WCS
+#ifndef NEWWCS
       Vector cdelt= ptr->getWCScdelt(sys);
-      area = fabs(cdelt[0]*cdelt[1]*cnt);
+      double aa = fabs(cdelt[0]*cdelt[1]);
+#else
+      double aa = ptr->getWCSPixelArea(sys);
+#endif
+      area = aa*cnt;
     }
     break;
   }