simplify wcs code

4 files changed, 83 insertions, 59 deletions

diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C
index b6b0a1a..0504708 100644
--- a/tksao/frame/fitsimage.C
+++ b/tksao/frame/fitsimage.C
@@ -1333,11 +1333,11 @@ void FitsImage::match(const char* xxname1, const char* yyname1,
 
   Vector* in1 = new Vector[nxx1];
   for (int ii=0; ii<nxx1; ii++)
-    in1[ii] = degToRad(Vector(ixx1[ii],iyy1[ii]));
+    in1[ii] = vDegToRad(Vector(ixx1[ii],iyy1[ii]),sys1);
 
   Vector* ptr2 = new Vector[nxx2];
   for (int ii=0; ii<nxx2; ii++)
-    ptr2[ii] = degToRad(Vector(ixx2[ii],iyy2[ii]));
+    ptr2[ii] = vDegToRad(Vector(ixx2[ii],iyy2[ii]),sys2);
 
   double rr = rad;
   switch (dist) {
@@ -2708,7 +2708,7 @@ double FitsImage::getWCSDist(const Vector& vv1, const Vector& vv2,
   astEnd;
   
   return hasWCSCel(sys) ?
-    radToDeg(wcsDistance(ast_,degToRad(vv1),degToRad(vv2))) :
+    radToDeg(wcsDistance(ast_,vDegToRad(vv1,sys),vDegToRad(vv2,sys))) :
     wcsDistance(ast_,vv1,vv2);
 }
 
@@ -2779,15 +2779,7 @@ Vector FitsImage::pix2wcs(const Vector& in, Coord::CoordSystem sys,
   astNorm(ast_, out.v);
   astEnd;
 
-  if (hasWCSCel(sys)) {
-    int id = sys-Coord::WCS;
-    if (wcsCelLon_[id] == 1 || wcsCelLat_[id] == 1)
-      out[0] = radToDeg(out[0]);
-    if (wcsCelLon_[id] == 2 || wcsCelLat_[id] == 2)
-      out[1] = radToDeg(out[1]);
-  }
-
-  return out;
+  return vRadToDeg(out,sys);
 }
 
 VectorStr FitsImage::pix2wcs(const Vector& in, Coord::CoordSystem sys,
@@ -2832,15 +2824,7 @@ Vector3d FitsImage::pix2wcs(const Vector3d& in, Coord::CoordSystem sys,
   astNorm(ast_, out.v);
   astEnd;
   
-  if (hasWCSCel(sys)) {
-    int id = sys-Coord::WCS;
-    if (wcsCelLon_[id] == 1 || wcsCelLat_[id] == 1)
-      out[0] = radToDeg(out[0]);
-    if (wcsCelLon_[id] == 2 || wcsCelLat_[id] == 2)
-      out[1] = radToDeg(out[1]);
-  }
-
-  return out;
+  return vRadToDeg(out,sys);
 }
 
 VectorStr3d FitsImage::pix2wcs(const Vector3d& in, Coord::CoordSystem sys,
@@ -2871,34 +2855,34 @@ VectorStr3d FitsImage::pix2wcs(const Vector3d& in, Coord::CoordSystem sys,
 Vector FitsImage::wcs2pix(const Vector& vv, Coord::CoordSystem sys,
 			  Coord::SkyFrame sky)
 {
-  astClearStatus; // just to make sure
-  astBegin; // start memory management
-
-  if (!(hasWCS(sys) && wcsInv_)) {
+  if (!hasWCS(sys) || !wcsInv_) {
     maperr =1;
     return Vector();
   }
 
+  astClearStatus; // just to make sure
+  astBegin; // start memory management
+
   setWCSSystem(sys);
   setWCSSkyFrame(sky);
 
-  Vector in = hasWCSCel(sys) ? degToRad(vv) : vv;
+  Vector in = vDegToRad(vv,sys);
   Vector out = wcsTran(ast_, in, 0);
-
   astEnd;
-    
-  if (astOK && checkWCS(out))
-    return out;
+  if (!astOK || !checkWCS(out)) {
+    maperr =1;
+    return Vector();
+  }
   
-  maperr =1;
-  return Vector();
+  maperr =0;
+  return out;
 }
 
 Vector3d FitsImage::wcs2pix(const Vector3d& vv, Coord::CoordSystem sys,
 			    Coord::SkyFrame sky)
 {
-  if (!(hasWCS(sys) && hasWCS3D(sys) && wcsInv_))
-  return Vector3d();
+  if (!hasWCS(sys) || !hasWCS3D(sys) || !wcsInv_)
+    return Vector3d();
     
   astClearStatus; // just to make sure
   astBegin; // start memory management
@@ -2906,15 +2890,13 @@ Vector3d FitsImage::wcs2pix(const Vector3d& vv, Coord::CoordSystem sys,
   setWCSSystem(sys);
   setWCSSkyFrame(sky);
 
-  Vector3d in = hasWCSCel(sys) ? degToRad(vv) : vv;
+  Vector3d in = vDegToRad(vv,sys);
   Vector3d out = wcsTran(ast_, in, 0);
-
   astEnd;
+  if (!astOK || !checkWCS(out))
+    return Vector3d();
   
-  if (astOK && checkWCS(out))
-    return out;
-  
-  return Vector3d();
+  return out;
 }
 
 int FitsImage::hasWCS(Coord::CoordSystem sys)
@@ -3236,6 +3218,62 @@ void FitsImage::setWCSFormat(Coord::CoordSystem sys, Coord::SkyFrame sky,
       wcsFormat(ast_, ii+1, str.str().c_str());
 }
 
+Vector FitsImage::vRadToDeg(const Vector& vv, Coord::CoordSystem sys)
+{
+  Vector out = vv;
+  if (hasWCSCel(sys)) {
+    int id = sys-Coord::WCS;
+    if (wcsCelLon_[id] == 1 || wcsCelLat_[id] == 1)
+      out[0] = radToDeg(out[0]);
+    if (wcsCelLon_[id] == 2 || wcsCelLat_[id] == 2)
+      out[1] = radToDeg(out[1]);
+  }
+  return out;
+}
+
+Vector3d FitsImage::vRadToDeg(const Vector3d& vv, Coord::CoordSystem sys)
+{
+  Vector3d out = vv;
+  if (hasWCSCel(sys)) {
+    int id = sys-Coord::WCS;
+    if (wcsCelLon_[id] == 1 || wcsCelLat_[id] == 1)
+      out[0] = radToDeg(out[0]);
+    if (wcsCelLon_[id] == 2 || wcsCelLat_[id] == 2)
+      out[1] = radToDeg(out[1]);
+    if (wcsCelLon_[id] == 2 || wcsCelLat_[id] == 2)
+      out[2] = radToDeg(out[2]);
+  }
+  return out;
+}
+
+Vector FitsImage::vDegToRad(const Vector& vv, Coord::CoordSystem sys)
+{
+  Vector out = vv;
+  if (hasWCSCel(sys)) {
+    int id = sys-Coord::WCS;
+    if (wcsCelLon_[id] == 1 || wcsCelLat_[id] == 1)
+      out[0] = degToRad(out[0]);
+    if (wcsCelLon_[id] == 2 || wcsCelLat_[id] == 2)
+      out[1] = degToRad(out[1]);
+  }
+  return out;
+}
+
+Vector3d FitsImage::vDegToRad(const Vector3d& vv, Coord::CoordSystem sys)
+{
+  Vector3d out = vv;
+  if (hasWCSCel(sys)) {
+    int id = sys-Coord::WCS;
+    if (wcsCelLon_[id] == 1 || wcsCelLat_[id] == 1)
+      out[0] = degToRad(out[0]);
+    if (wcsCelLon_[id] == 2 || wcsCelLat_[id] == 2)
+      out[1] = degToRad(out[1]);
+    if (wcsCelLon_[id] == 2 || wcsCelLat_[id] == 2)
+      out[2] = degToRad(out[2]);
+  }
+  return out;
+}
+
 static void fits2TAB(AstFitsChan* chan, const char* extname,
 		     int extver, int extlevel, int* status)
 {
diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h
index 621773a..5ae1afc 100644
--- a/tksao/frame/fitsimage.h
+++ b/tksao/frame/fitsimage.h
@@ -167,6 +167,11 @@ class FitsImage {
   void setWCSSkyFrame(Coord::SkyFrame);
   void setWCSFormat(Coord::CoordSystem, Coord::SkyFrame, Coord::SkyFormat);
 
+  Vector vRadToDeg(const Vector& vv, Coord::CoordSystem sys);
+  Vector vDegToRad(const Vector& vv, Coord::CoordSystem sys);
+  Vector3d vRadToDeg(const Vector3d& vv, Coord::CoordSystem sys);
+  Vector3d vDegToRad(const Vector3d& vv, Coord::CoordSystem sys);
+
  public:
   char* fileName;
   char* rootBaseFileName;
diff --git a/tksao/util/util.C b/tksao/util/util.C
index 11d692f..faade52 100644
--- a/tksao/util/util.C
+++ b/tksao/util/util.C
@@ -150,22 +150,6 @@ double degToRad(double aa)
   return  M_PI*aa/180.;
 }
 
-Vector degToRad(const Vector& vv)
-{
-  Vector out =vv;
-  out[0] *= M_PI/180.;
-  out[1] *= M_PI/180.;
-  return out;
-}
-
-Vector3d degToRad(const Vector3d& vv)
-{
-  Vector3d out =vv;
-  out[0] *= M_PI/180.;
-  out[1] *= M_PI/180.;
-  return out;
-}
-
 double dmsToDegree(int sign, int degree, int min, double sec)
 {
   // sign is needed because of -00 vs +00
diff --git a/tksao/util/util.h b/tksao/util/util.h
index 8b5247c..3609ae6 100644
--- a/tksao/util/util.h
+++ b/tksao/util/util.h
@@ -92,9 +92,6 @@ extern double m180To180(double);
 extern double degToRad(double);
 extern double radToDeg(double);
 
-extern Vector degToRad(const Vector&);
-extern Vector3d degToRad(const Vector3d&);
-
 extern int parseSection(char*, Vector*, Vector*);
 
 extern double dmsToDegree(int, int, int, double);