From 75488234583582af8356db0c54efe806851aa791 Mon Sep 17 00:00:00 2001
From: William Joye <wjoye@cfa.harvard.edu>
Date: Wed, 13 Dec 2017 14:06:43 -0500
Subject: upgrade AST

---
 tksao/frame/fitsimage.C | 78 +++++++++++++++++++++++++++++++------------------
 tksao/frame/fitsimage.h | 23 +++++++--------
 2 files changed, 59 insertions(+), 42 deletions(-)

diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C
index 6ba97ad..28de1dc 100644
--- a/tksao/frame/fitsimage.C
+++ b/tksao/frame/fitsimage.C
@@ -2971,10 +2971,7 @@ double FitsImage::getWCSPixelSize(Coord::CoordSystem sys)
   wcsTran(3, in, 1, out);
   double dd = (wcsDistance(out[0],out[1]) + wcsDistance(out[0],out[2]))/2.;
 
-  if (wcsIsASkyFrame())
-    return radToDeg(dd);
-  else
-    return dd;
+  return wcsIsASkyFrame() ? radToDeg(dd) : dd;
 }
 
 double FitsImage::getWCSPixelArea(Coord::CoordSystem sys)
@@ -2994,10 +2991,7 @@ double FitsImage::getWCSPixelArea(Coord::CoordSystem sys)
   double ll = wcsDistance(out[0], out[1]);
   double mm = wcsDistance(out[0], out[2]);
 
-  if (wcsIsASkyFrame())
-    return radToDeg(ll)*radToDeg(mm);
-  else
-    return ll*mm;
+  return wcsIsASkyFrame() ? radToDeg(ll)*radToDeg(mm) : ll*mm;
 }
 #endif
 
@@ -3102,17 +3096,16 @@ double FitsImage::getWCSRotation(Coord::CoordSystem sys, Coord::SkyFrame sky)
   in[1] = center()+Vector(0,1);
   wcsTran(2, in, 1, out);
   double ang = wcsAxAngle(out[0], out[1]);
-
-  //  {
-  //    Vector npix = wcsTran(out[0]+Vector(0,.01),0);
-  //    Vector epix = wcsTran(out[0]+Vector(.01,0),0);
-  //    Vector north = (npix-in[0]).normalize();
-  //    Vector east = (epix-in[0]).normalize();
-  //    Vector diff = (north-east).abs();
-  //    cerr << diff << endl;
-  //    double bb = -(north.angle()-M_PI_2);
-  //    double aa = (getWCSOrientation(sys,sky) == Coord::NORMAL) ? ang : -ang;
-  //  }  
+  {
+    //    Vector npix = wcsTran(out[0]+Vector(0,.1),0);
+    //    Vector epix = wcsTran(out[0]+Vector(.1,0),0);
+    //    Vector north = (npix-in[0]).normalize();
+    //    Vector east = (epix-in[0]).normalize();
+    //    Vector diff = (north-east).abs();
+    //    double bb = -(north.angle()-M_PI_2);
+    //    double aa = getWCSOrientation(sys,sky) == Coord::NORMAL ? ang : -ang;
+    //    cerr << radToDeg(aa) << '=' << radToDeg(bb)  << endl;
+  }  
 
   if (!(isnan(ang)||isinf(ang)||(ang == -DBL_MAX)||(ang == DBL_MAX)))
     return getWCSOrientation(sys,sky) == Coord::NORMAL ? ang : -ang;
@@ -3265,15 +3258,15 @@ char* FitsImage::pix2wcs(const Vector& in, Coord::CoordSystem sys,
 	case Coord::FK4_NO_E:
 	case Coord::FK5:
 	case Coord::ICRS:
-	  setWCSFormat(ast_,1,"hms.3");
-	  setWCSFormat(ast_,2,"+dms.3");
+	  setWCSFormat(1,"hms.3");
+	  setWCSFormat(2,"+dms.3");
 	  break;
 	case Coord::GALACTIC:
 	case Coord::SUPERGALACTIC:
 	case Coord::ECLIPTIC:
 	case Coord::HELIOECLIPTIC:
-	  setWCSFormat(ast_,1,"+dms.3");
-	  setWCSFormat(ast_,2,"+dms.3");
+	  setWCSFormat(1,"+dms.3");
+	  setWCSFormat(2,"+dms.3");
 	  break;
 	}
 	str << astFormat(ast_,1,out[0]) << ' '
@@ -3342,15 +3335,15 @@ char* FitsImage::pix2wcs(const Vector3d& in, Coord::CoordSystem sys,
 	case Coord::FK4_NO_E:
 	case Coord::FK5:
 	case Coord::ICRS:
-	  setWCSFormat(ast_,1,"hms.3");
-	  setWCSFormat(ast_,2,"+dms.3");
+	  setWCSFormat(1,"hms.3");
+	  setWCSFormat(2,"+dms.3");
 	  break;
 	case Coord::GALACTIC:
 	case Coord::SUPERGALACTIC:
 	case Coord::ECLIPTIC:
 	case Coord::HELIOECLIPTIC:
-	  setWCSFormat(ast_,1,"+dms.3");
-	  setWCSFormat(ast_,2,"+dms.3");
+	  setWCSFormat(1,"+dms.3");
+	  setWCSFormat(2,"+dms.3");
 	  break;
 	}
 	str << astFormat(ast_,1,out[0]) << ' '
@@ -3682,6 +3675,11 @@ void FitsImage::astInit(FitsHead* hd, FitsHead* prim)
   case 4:
     break;
   }
+
+  if (wcsIsASkyFrame()) {
+    //    int base = astGetI(ast_,"Base");
+    //    int current = astGetI(ast_,"Current");
+  }
 }
 
 void FitsImage::wcsInit()
@@ -3704,6 +3702,7 @@ void FitsImage::wcsInit()
 
   // since we have ast_
   wcs_[0] =1; 
+
   int nn = astGetI(ast_, "Nframe");
   for (int ii=0; ii<nn; ii++) {
     const char* id = astGetC(astGetFrame(ast_,ii+1), "Ident");
@@ -3733,8 +3732,9 @@ void FitsImage::wcsEquInit()
   astClearStatus;
   astBegin;
 
-  // since we have ast_
+  // for AST WCS
   wcsEqu_[0] = astIsASkyFrame(astGetFrame(ast_, AST__CURRENT));
+
   int nn = astGetI(ast_, "Nframe");
   for (int ii=0; ii<nn; ii++) {
     AstFrame* ff = (AstFrame*)astGetFrame(ast_,ii+1);
@@ -3783,8 +3783,9 @@ void FitsImage::wcsCelInit()
   astClearStatus;
   astBegin;
 
-  // since we have ast_
+  // for AST WCS
   wcsCel_[0] = astIsASkyFrame(astGetFrame(ast_, AST__CURRENT));
+
   int nn = astGetI(ast_, "Nframe");
   for (int ii=0; ii<nn; ii++) {
     AstFrame* ff = (AstFrame*)astGetFrame(ast_,ii+1);
@@ -3875,6 +3876,7 @@ int FitsImage::checkWCS(Vector3d& vv)
 	  fabs(vv[2]) < FLT_MAX ) ? 1 : 0;
 }
 
+#ifndef NEWWCS
 void FitsImage::setWCSFormat(AstFrameSet* aa, int id, const char* format)
 {
   // is it already set?
@@ -3891,6 +3893,24 @@ void FitsImage::setWCSFormat(AstFrameSet* aa, int id, const char* format)
   str << "Format(" << id << ")=" << format << ends;
   astSet(aa, str.str().c_str());
 }
+#else
+void FitsImage::setWCSFormat(int id, const char* format)
+{
+  // is it already set?
+  // ast is very slow when changing params
+  {
+    ostringstream str;
+    str << "Format(" << id << ")" << ends;
+    const char* out = astGetC(ast_, str.str().c_str());
+    if (!strcmp(out,format))
+      return;
+  }
+
+  ostringstream str;
+  str << "Format(" << id << ")=" << format << ends;
+  astSet(ast_, str.str().c_str());
+}
+#endif
 
 #ifndef NEWWCS
 void FitsImage::setWCSSkyFrame(AstFrameSet* ast, Coord::SkyFrame sky)
diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h
index 4a57975..03ed620 100644
--- a/tksao/frame/fitsimage.h
+++ b/tksao/frame/fitsimage.h
@@ -401,30 +401,25 @@ class FitsImage {
   void processKeywordsFitsSection();
   int processKeywordsIRAF(FitsImage*);
 
-#ifndef NEWWCS
-  WorldCoor* getWCS(Coord::CoordSystem sys) 
-  {return (wcs_ && wcs_[sys-Coord::WCS]) ? wcs_[sys-Coord::WCS] : NULL;}
-#endif
-  const char* getWCSName(Coord::CoordSystem);
   Coord::Orientation getWCSOrientation(Coord::CoordSystem, Coord::SkyFrame);
   double getWCSRotation(Coord::CoordSystem, Coord::SkyFrame);
   double getWCSDist(const Vector&, const Vector&, Coord::CoordSystem);
+  const char* getWCSName(Coord::CoordSystem);
+
 #ifndef NEWWCS
+  WorldCoor* getWCS(Coord::CoordSystem sys) 
+  {return (wcs_ && wcs_[sys-Coord::WCS]) ? wcs_[sys-Coord::WCS] : NULL;}
   Vector getWCScdelt(Coord::CoordSystem);
-#else
-  double getWCSPixelSize(Coord::CoordSystem);
-  double getWCSPixelArea(Coord::CoordSystem);
-#endif
-  
-#ifndef NEWWCS
   AstFrameSet* getAST(Coord::CoordSystem sys) 
     {return (ast_ && ast_[sys-Coord::WCS]) ? ast_[sys-Coord::WCS] : NULL;}
+
   Vector wcsTran(AstFrameSet*, const Vector&, int);
   void wcsTran(AstFrameSet*, int, Vector*, int, Vector*);
   double wcsDistance(AstFrameSet*, const Vector&, const Vector&);
 
   int wcsIsASkyFrame(AstFrameSet*);
   void setWCSSkyFrame(AstFrameSet*, Coord::SkyFrame);
+  void setWCSFormat(AstFrameSet*, int, const char*);
 #else
   AstFrameSet* wcsCopy() {return (AstFrameSet*)astCopy(ast_);}
   Vector wcsTran(const Vector&, int);
@@ -439,9 +434,11 @@ class FitsImage {
   int wcsIsASkyFrame();
   void setWCSSystem(Coord::CoordSystem);
   void setWCSSkyFrame(Coord::SkyFrame);
-#endif
+  void setWCSFormat(int, const char*);
 
-  void setWCSFormat(AstFrameSet*, int, const char*);
+  double getWCSPixelSize(Coord::CoordSystem);
+  double getWCSPixelArea(Coord::CoordSystem);
+#endif
 
   int hasWCS(Coord::CoordSystem);
   int hasWCSEqu(Coord::CoordSystem);
-- 
cgit v0.12