simplify mapLen code

4 files changed, 33 insertions, 56 deletions

diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C
index 2655be9..71ec4b1 100644
--- a/tksao/frame/fitsimage.C
+++ b/tksao/frame/fitsimage.C
@@ -1434,7 +1434,7 @@ void FitsImage::initWCS0(const Vector& pix)
   hd->appendReal("CRVAL2", 0, 8, NULL);
 
   // CD
-  Vector ss = getWCSSize(Coord::WCS);
+  double ss = getWCSSize(Coord::WCS);
   double ang = getWCSRotation(Coord::WCS,Coord::FK5);
   Matrix flip;
   switch (getWCSOrientation(Coord::WCS,Coord::FK5)) {
@@ -3056,7 +3056,7 @@ Vector FitsImage::getWCScdelt(Coord::CoordSystem sys)
 #endif
 
 #ifndef OLDWCS
-Vector FitsImage::getWCSSize(Coord::CoordSystem sys)
+double FitsImage::getWCSSize(Coord::CoordSystem sys)
 {
   if (!hasWCS(sys))
     return 0;
@@ -3064,14 +3064,14 @@ Vector FitsImage::getWCSSize(Coord::CoordSystem sys)
   astClearStatus; // just to make sure
   setWCSSkyFrame(sys, Coord::FK5);
   
+  // only check lon axis as it will not change near poles
   Vector in[3];
   Vector out[3];
   in[0] = center();
-  in[1] = center()+Vector(1,0);
-  in[2] = center()+Vector(0,1);
-  wcsTran(3, in, 1, out);
+  in[1] = center()+Vector(0,1);
+  wcsTran(2, in, 1, out);
 
-  Vector dd = Vector(wcsDistance(out[0],out[1]), wcsDistance(out[0],out[2]));
+  double dd = wcsDistance(out[0],out[1]);
   return hasWCSCel(sys) ? radToDeg(dd) : dd;
 }
 #endif
diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h
index 071d77c..2ad85b8 100644
--- a/tksao/frame/fitsimage.h
+++ b/tksao/frame/fitsimage.h
@@ -451,7 +451,7 @@ class FitsImage {
   void setWCSSkyFrame(Coord::CoordSystem, Coord::SkyFrame);
   void setWCSFormat(int, const char*);
 
-  Vector getWCSSize(Coord::CoordSystem);
+  double getWCSSize(Coord::CoordSystem);
 #endif
 
   int hasWCS(Coord::CoordSystem);
diff --git a/tksao/frame/fitsmap.C b/tksao/frame/fitsmap.C
index e1cc229..17352f1 100644
--- a/tksao/frame/fitsmap.C
+++ b/tksao/frame/fitsmap.C
@@ -247,14 +247,14 @@ Vector FitsImage::mapLenFromRef(const Vector& vv, Coord::CoordSystem sys,
 }
 #endif
 
-#ifdef OLDWCS
 double FitsImage::mapLenToRef(double dd, Coord::CoordSystem sys, 
 			      Coord::DistFormat dist)
 {
-  Vector rr = mapLenToRef(Vector(dd,0), sys, dist);
-  return rr[0];
+  Vector rr = mapLenToRef(Vector(0,dd), sys, dist);
+  return rr[1];
 }
 
+#ifdef OLDWCS
 Vector FitsImage::mapLenToRef(const Vector& vv, Coord::CoordSystem sys,
 			      Coord::DistFormat dist)
 {
@@ -292,63 +292,43 @@ Vector FitsImage::mapLenToRef(const Vector& vv, Coord::CoordSystem sys,
   return Vector();
 }
 #else
-double FitsImage::mapLenToRef(double dd, Coord::CoordSystem sys, 
+Vector FitsImage::mapLenToRef(const Vector& vv, Coord::CoordSystem sys,
 			      Coord::DistFormat dist)
 {
+  // this is correct
+  //   we are looking for a length, so in image coords
   switch (sys) {
   case Coord::IMAGE:
-    return dd*imageToRef[1].length();
+    return vv;
   case Coord::PHYSICAL:
-    return dd*physicalToRef[1].length();
+    return vv*physicalToImage;
   case Coord::AMPLIFIER:
-    return dd*(amplifierToPhysical * physicalToRef)[1].length();
+    return vv*amplifierToPhysical*physicalToImage;
   case Coord::DETECTOR:
-    return dd*(detectorToPhysical * physicalToRef)[1].length();
+    return vv*detectorToPhysical*physicalToImage;
   default:
     {
       if (!hasWCS(sys) || !astInv_)
 	return 0;
-      
-      astClearStatus; // just to make sure
-      setWCSSkyFrame(sys, Coord::FK5);
 
-      double rdd = dd;
+      Vector svv = vv;
       if (hasWCSCel(sys)) {
-	rdd = dd;
 	switch (dist) {
 	case Coord::DEGREE:
 	  break;
 	case Coord::ARCMIN:
-	  rdd /= 60.;
+	  svv /= 60.;
 	  break;
 	case Coord::ARCSEC:
-	  rdd /= 60.*60.;
+	  svv /= 60.*60.;
 	  break;
 	}
-	rdd = degToRad(rdd); // no zeroTWOPI since this is a length
       }
-
-      Vector cc = center();
-      Vector wcc = wcsTran(cc,1);
-      Vector wpp = wcc+Vector(0,rdd);
-      Vector pp = wcsTran(wpp,0);
-      astInvert(ast_);
-      double rr = wcsDistance(cc, pp);
-      astInvert(ast_);
-
-      return rr;
+      return svv/getWCSSize(sys);
     }
   }
 
-  return 0;
-}
-
-Vector FitsImage::mapLenToRef(const Vector& vv, Coord::CoordSystem sys,
-			      Coord::DistFormat dist)
-{
-  double rx = mapLenToRef(((Vector)vv)[0],sys,dist);
-  double ry = mapLenToRef(((Vector)vv)[1],sys,dist);
-  return Vector(rx,ry);
+  return Vector();
 }
 #endif
 
diff --git a/tksao/frame/frblt.C b/tksao/frame/frblt.C
index f9aa4fc..f79df87 100644
--- a/tksao/frame/frblt.C
+++ b/tksao/frame/frblt.C
@@ -376,8 +376,7 @@ int Base::markerAnalysisRadial(Marker* pp, double** x, double** y, double** e,
     Vector cdelt= ptr->getWCScdelt(sys);
     double ll = fabs(cdelt[0]);
 #else
-    Vector rr = ptr->getWCSSize(sys);
-    double ll = (rr[0]+rr[1])/2.;
+    double ll = ptr->getWCSSize(sys);
 #endif
 
     if (ptr->hasWCSCel(sys)) {
@@ -394,8 +393,8 @@ int Base::markerAnalysisRadial(Marker* pp, double** x, double** y, double** e,
   Vector cdelt= ptr->getWCScdelt(sys);
   double aa = fabs(cdelt[0]*cdelt[1]);
 #else
-  Vector rr = ptr->getWCSSize(sys);
-  double aa = rr.area();
+  double rr = ptr->getWCSSize(sys);
+  double aa = rr*rr;
 #endif
 
   for (int kk=0; kk<num; kk++) {
@@ -501,8 +500,7 @@ int Base::markerAnalysisPanda(Marker* pp, double** x, double** y, double** e,
     Vector cdelt= ptr->getWCScdelt(sys);
     double ll = fabs(cdelt[0]);
 #else
-    Vector rr = ptr->getWCSSize(sys);
-    double ll = (rr[0]+rr[1])/2.;
+    double ll = ptr->getWCSSize(sys);
 #endif
 
     if (ptr->hasWCSCel(sys)) {
@@ -519,8 +517,8 @@ int Base::markerAnalysisPanda(Marker* pp, double** x, double** y, double** e,
   Vector cdelt= ptr->getWCScdelt(sys);
   double aa = fabs(cdelt[0]*cdelt[1]);
 #else
-  Vector rr = ptr->getWCSSize(sys);
-  double aa = rr.area();
+  double rr = ptr->getWCSSize(sys);
+  double aa = rr*rr;
 #endif
 
   for (int qq=0; qq<aa; qq++) {
@@ -840,8 +838,7 @@ int Base::markerAnalysisStats1(Marker* pp,FitsImage* ptr, ostream& str,
       Vector cdelt= ptr->getWCScdelt(sys);
       double ll = fabs(cdelt[0]);
 #else
-      Vector rr = ptr->getWCSSize(sys);
-      double ll = (rr[0]+rr[1])/2.;
+      double ll = ptr->getWCSSize(sys);
 #endif
       if (ptr->hasWCSCel(sys)) {
 	str << "1 pixel = "<< ll*60*60 << " arcsec";
@@ -887,8 +884,8 @@ void Base::markerAnalysisStats2(FitsImage* ptr, ostream& str,
       Vector cdelt= ptr->getWCScdelt(sys);
       double aa = fabs(cdelt[0]*cdelt[1]);
 #else
-      Vector rr = ptr->getWCSSize(sys);
-      double aa = rr.area();
+      double rr = ptr->getWCSSize(sys);
+      double aa = rr*rr;
 #endif
       area = aa*60*60*60*60*cnt;
     }
@@ -900,8 +897,8 @@ void Base::markerAnalysisStats2(FitsImage* ptr, ostream& str,
       Vector cdelt= ptr->getWCScdelt(sys);
       double aa = fabs(cdelt[0]*cdelt[1]);
 #else
-      Vector rr = ptr->getWCSSize(sys);
-      double aa = rr.area();
+      double rr = ptr->getWCSSize(sys);
+      double aa = rr*rr;
 #endif
       area = aa*cnt;
     }