summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rwxr-xr-xtksao/configure1
-rw-r--r--tksao/configure.ac1
-rw-r--r--tksao/frame/base.C9
-rw-r--r--tksao/frame/base.h6
-rw-r--r--tksao/frame/contourparser.C4
-rw-r--r--tksao/frame/contourparser.Y4
-rw-r--r--tksao/frame/ds9parser.C4
-rw-r--r--tksao/frame/ds9parser.Y4
-rw-r--r--tksao/frame/fitsimage.C383
-rw-r--r--tksao/frame/fitsimage.h11
-rw-r--r--tksao/frame/grid25d.C9
-rw-r--r--tksao/frame/grid2d.C9
-rw-r--r--tksao/frame/grid3d.C9
-rw-r--r--tksao/frame/wcsast.C384
-rw-r--r--tksao/frame/wcsast.h22
15 files changed, 434 insertions, 426 deletions
diff --git a/tksao/configure b/tksao/configure
index ff3056c..629c597 100755
--- a/tksao/configure
+++ b/tksao/configure
@@ -5808,6 +5808,7 @@ frame/text.C
frame/tnglex.C
frame/tngparser.C
frame/vect.C
+frame/wcsast.C
frame/xylex.C
frame/xyparser.C
list/list.C
diff --git a/tksao/configure.ac b/tksao/configure.ac
index 069a4f9..c5bfeba 100644
--- a/tksao/configure.ac
+++ b/tksao/configure.ac
@@ -237,6 +237,7 @@ frame/text.C
frame/tnglex.C
frame/tngparser.C
frame/vect.C
+frame/wcsast.C
frame/xylex.C
frame/xyparser.C
list/list.C
diff --git a/tksao/frame/base.C b/tksao/frame/base.C
index afacd6b..de5afa6 100644
--- a/tksao/frame/base.C
+++ b/tksao/frame/base.C
@@ -8,6 +8,7 @@
#include "context.h"
#include "fitsimage.h"
#include "marker.h"
+#include "wcsast.h"
#include "ps.h"
#include "circle.h"
@@ -388,13 +389,13 @@ Matrix Base::calcAlignWCS(FitsImage* fits1, FitsImage* fits2,
astBegin; // start memory management
AstFrameSet* wcs1 = (AstFrameSet*)astCopy(fits1->ast_);
- fits1->wcsSystem(wcs1,sys1);
- fits1->wcsSkyFrame(wcs1,sky);
+ wcsSystem(wcs1,sys1);
+ wcsSkyFrame(wcs1,sky);
astInvert(wcs1);
AstFrameSet* wcs2 = (AstFrameSet*)astCopy(fits2->ast_);
- fits2->wcsSystem(wcs2,sys1);
- fits2->wcsSkyFrame(wcs2,sky);
+ wcsSystem(wcs2,sys1);
+ wcsSkyFrame(wcs2,sky);
astInvert(wcs2);
AstFrameSet* cvt = (AstFrameSet*)astConvert(wcs2, wcs1, "");
diff --git a/tksao/frame/base.h b/tksao/frame/base.h
index 931681b..4de5d18 100644
--- a/tksao/frame/base.h
+++ b/tksao/frame/base.h
@@ -550,9 +550,9 @@ public:
void resetCompositeMarker() {compositeMarker = NULL;}
- Coord::CoordSystem wcsSystem() {return wcsSystem_;}
- Coord::SkyFrame wcsSkyFrame() {return wcsSkyFrame_;}
- Coord::SkyFormat wcsSkyFormat() {return wcsSkyFormat_;}
+ Coord::CoordSystem getWCSSystem() {return wcsSystem_;}
+ Coord::SkyFrame getWCSSkyFrame() {return wcsSkyFrame_;}
+ Coord::SkyFormat getWCSSkyFormat() {return wcsSkyFormat_;}
Coord::CoordSystem xySystem() {return xySystem_;}
Coord::SkyFrame xySky() {return xySky_;}
diff --git a/tksao/frame/contourparser.C b/tksao/frame/contourparser.C
index f98f807..9f5adb8 100644
--- a/tksao/frame/contourparser.C
+++ b/tksao/frame/contourparser.C
@@ -1879,8 +1879,8 @@ yyreduce:
cl = NULL;
cc = NULL;
globalSystem = Coord::WCS;
- globalWCS = fr->wcsSystem();
- globalSky = fr->wcsSkyFrame();
+ globalWCS = fr->getWCSSystem();
+ globalSky = fr->getWCSSkyFrame();
strcpy(globalColor,"green");
globalDash = 0;
globalDashList[0] = 8;
diff --git a/tksao/frame/contourparser.Y b/tksao/frame/contourparser.Y
index 07d4d09..9f7742c 100644
--- a/tksao/frame/contourparser.Y
+++ b/tksao/frame/contourparser.Y
@@ -259,8 +259,8 @@ initGlobal: {
cl = NULL;
cc = NULL;
globalSystem = Coord::WCS;
- globalWCS = fr->wcsSystem();
- globalSky = fr->wcsSkyFrame();
+ globalWCS = fr->getWCSSystem();
+ globalSky = fr->getWCSSkyFrame();
strcpy(globalColor,"green");
globalDash = 0;
globalDashList[0] = 8;
diff --git a/tksao/frame/ds9parser.C b/tksao/frame/ds9parser.C
index 9433062..d4c26a0 100644
--- a/tksao/frame/ds9parser.C
+++ b/tksao/frame/ds9parser.C
@@ -3788,8 +3788,8 @@ yyreduce:
{
// global properties
globalSystem = Coord::PHYSICAL;
- globalWCS = fr->wcsSystem();
- globalSky = fr->wcsSkyFrame();
+ globalWCS = fr->getWCSSystem();
+ globalSky = fr->getWCSSkyFrame();
globalTile = 1;
globalProps =
Marker::SELECT | Marker::EDIT | Marker::MOVE |
diff --git a/tksao/frame/ds9parser.Y b/tksao/frame/ds9parser.Y
index 07055a0..2319f3b 100644
--- a/tksao/frame/ds9parser.Y
+++ b/tksao/frame/ds9parser.Y
@@ -775,8 +775,8 @@ globalCompass : coordSystem skyFrame
initGlobal:{
// global properties
globalSystem = Coord::PHYSICAL;
- globalWCS = fr->wcsSystem();
- globalSky = fr->wcsSkyFrame();
+ globalWCS = fr->getWCSSystem();
+ globalSky = fr->getWCSSkyFrame();
globalTile = 1;
globalProps =
Marker::SELECT | Marker::EDIT | Marker::MOVE |
diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C
index e3d6a4a..4ce6608 100644
--- a/tksao/frame/fitsimage.C
+++ b/tksao/frame/fitsimage.C
@@ -26,6 +26,7 @@
#include "analysis.h"
#include "photo.h"
#include "colorscale.h"
+#include "wcsast.h"
FitsImage::FitsImage(Context* cx, Tcl_Interp* pp)
{
@@ -3256,388 +3257,6 @@ void FitsImage::setWCSSkyFrame(Coord::CoordSystem sys, Coord::SkyFrame sky)
}
}
-int FitsImage::wcsSystem(AstFrameSet* ast, Coord::CoordSystem sys)
-{
- int nn = astGetI(ast,"nframe");
- char cc = ' ';
- int ww = sys-Coord::WCS;
- if (ww < 0)
- return 0;
- if (ww)
- cc = ww+'@';
-
- for (int ss=0; ss<nn; ss++) {
- const char* id = astGetC(astGetFrame(ast,ss+1),"Ident");
- if (cc == id[0]) {
- astSetI(ast,"Current",ss+1);
- return 1;
- }
- }
-
- return 0;
-}
-
-int FitsImage::wcsSkyFrame(AstFrameSet* ast, Coord::SkyFrame sky)
-{
- AstFrameSet* fs =
- (AstFrameSet*)astFindFrame(ast, astSkyFrame(" MaxAxes=4")," ");
-
- if (!fs)
- return 0;
-
- switch (sky) {
- case Coord::FK4:
- astSet(fs, "System=FK4, Equinox=B1950");
- break;
- case Coord::FK5:
- astSet(fs, "System=FK5, Equinox=J2000");
- break;
- case Coord::ICRS:
- astSet(fs, "System=ICRS");
- break;
- case Coord::GALACTIC:
- astSet(fs, "System=GALACTIC");
- break;
- case Coord::ECLIPTIC:
- astSet(fs, "System=ECLIPTIC");
- // get AST to agree with WCSSUBS
- astSetD(fs, "EQUINOX", astGetD(fs, "EPOCH"));
- break;
- }
-
- return 1;
-}
-
-Vector FitsImage::wcsTran(AstFrameSet* ast, const Vector& in, int forward)
-{
- int naxes = astGetI(ast,"Naxes");
- switch (naxes) {
- case 1:
- // error
- break;
- case 2:
- double xout, yout;
- astTran2(ast, 1, in.v, in.v+1, forward, &xout, &yout);
- return Vector(xout, yout);
- case 3:
- {
- double pin[3];
- double pout[3];
- pin[0] = in[0];
- pin[1] = in[1];
- pin[2] = forward ? context_->slice(2) : 0;
- astTranN(ast, 1, 3, 1, pin, forward, 3, 1, pout);
- return Vector(pout[0],pout[1]);
- }
- break;
- case 4:
- {
- double pin[4];
- double pout[4];
- pin[0] = in[0];
- pin[1] = in[1];
- pin[2] = forward ? context_->slice(2) : 0;
- pin[3] = forward ? context_->slice(3) : 0;
- astTranN(ast, 1, 4, 1, pin, forward, 4, 1, pout);
- return Vector(pout[0],pout[1]);
- }
- break;
- }
- return Vector();
-}
-
-void FitsImage::wcsTran(AstFrameSet* ast, int npoint,
- Vector* in, int forward, Vector* out)
-{
- int naxes = astGetI(ast,"Naxes");
- switch (naxes) {
- case 1:
- // error
- break;
- case 2:
- {
- double* xin = new double[npoint];
- double* yin = new double[npoint];
- double* xout = new double[npoint];
- double* yout = new double[npoint];
-
- for (int ii=0; ii<npoint; ii++) {
- xin[ii] = in[ii][0];
- yin[ii] = in[ii][1];
- }
- astTran2(ast, npoint, xin, yin, forward, xout, yout);
- for (int ii=0; ii<npoint; ii++)
- out[ii] = Vector(xout[ii],yout[ii]);
-
- if (xin)
- delete [] xin;
- if (yin)
- delete [] yin;
- if (xout)
- delete [] xout;
- if (yout)
- delete [] yout;
- }
- break;
- case 3:
- {
- double* ptr_in[3];
- ptr_in[0] = new double[npoint];
- ptr_in[1] = new double[npoint];
- ptr_in[2] = new double[npoint];
- double* ptr_out[3];
- ptr_out[0] = new double[npoint];
- ptr_out[1] = new double[npoint];
- ptr_out[2] = new double[npoint];
-
- for (int kk=0; kk<npoint; kk++) {
- ptr_in[0][kk] = in[kk][0];
- ptr_in[1][kk] = in[kk][1];
- ptr_in[2][kk] = forward ? context_->slice(2) : 0;
- }
- astTranP(ast, npoint, 3, (const double**)ptr_in, forward, 3, ptr_out);
- for (int kk=0; kk<npoint; kk++)
- out[kk] = Vector(ptr_out[0][kk], ptr_out[1][kk]);
-
- if (ptr_in[0])
- delete [] ptr_in[0];
- if (ptr_in[1])
- delete [] ptr_in[1];
- if (ptr_in[2])
- delete [] ptr_in[2];
-
- if (ptr_out[0])
- delete [] ptr_out[0];
- if (ptr_out[1])
- delete [] ptr_out[1];
- if (ptr_out[2])
- delete [] ptr_out[2];
- }
- break;
- case 4:
- {
- double* ptr_in[4];
- ptr_in[0] = new double[npoint];
- ptr_in[1] = new double[npoint];
- ptr_in[2] = new double[npoint];
- ptr_in[3] = new double[npoint];
- double* ptr_out[4];
- ptr_out[0] = new double[npoint];
- ptr_out[1] = new double[npoint];
- ptr_out[2] = new double[npoint];
- ptr_out[3] = new double[npoint];
-
- for (int kk=0; kk<npoint; kk++) {
- ptr_in[0][kk] = in[kk][0];
- ptr_in[1][kk] = in[kk][1];
- ptr_in[2][kk] = forward ? context_->slice(2) : 0;
- ptr_in[3][kk] = forward ? context_->slice(3) : 0;
- }
- astTranP(ast, npoint, 4, (const double**)ptr_in, forward, 4, ptr_out);
- for (int kk=0; kk<npoint; kk++)
- out[kk] = Vector(ptr_out[0][kk], ptr_out[1][kk]);
-
- if (ptr_in[0])
- delete [] ptr_in[0];
- if (ptr_in[1])
- delete [] ptr_in[1];
- if (ptr_in[2])
- delete [] ptr_in[2];
- if (ptr_in[3])
- delete [] ptr_in[3];
-
- if (ptr_out[0])
- delete [] ptr_out[0];
- if (ptr_out[1])
- delete [] ptr_out[1];
- if (ptr_out[2])
- delete [] ptr_out[2];
- if (ptr_out[3])
- delete [] ptr_out[3];
- }
- break;
- }
-}
-
-Vector3d FitsImage::wcsTran(AstFrameSet* ast, const Vector3d& in, int forward)
-{
- int naxes = astGetI(ast,"Naxes");
- switch (naxes) {
- case 1:
- case 2:
- double pin[2];
- double pout[2];
- pin[0] = in[0];
- pin[1] = in[1];
- astTranN(ast, 1, 2, 1, pin, forward, 2, 1, pout);
- return Vector3d(pout[0],pout[1],forward ? 1 : 0);
- break;
- case 3:
- {
- double pin[3];
- double pout[3];
- pin[0] = in[0];
- pin[1] = in[1];
- pin[2] = in[2];
- astTranN(ast, 1, 3, 1, pin, forward, 3, 1, pout);
- return Vector3d(pout[0],pout[1],pout[2]);
- }
- break;
- case 4:
- {
- double pin[4];
- double pout[4];
- pin[0] = in[0];
- pin[1] = in[1];
- pin[2] = in[2];
- pin[3] = forward ? context_->slice(3) : 0;
- astTranN(ast, 1, 4, 1, pin, forward, 4, 1, pout);
- return Vector3d(pout[0],pout[1],pout[2]);
- }
- break;
- }
- return Vector3d();
-}
-
-double FitsImage::wcsDistance(AstFrameSet* ast,
- const Vector& vv1, const Vector& vv2)
-{
- int naxes = astGetI(ast,"Naxes");
- switch (naxes) {
- case 1:
- // error
- break;
- case 2:
- return astDistance(ast, vv1.v, vv2.v);
- 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);
- }
- 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);
- }
- }
-
- return 0;
-}
-
-double FitsImage::wcsAngle(AstFrameSet* ast,
- const Vector& vv1, const Vector& vv2,
- const Vector& vv3)
-{
- int naxes = astGetI(ast,"Naxes");
- switch (naxes) {
- case 1:
- // error
- break;
- case 2:
- return astAngle(ast,vv1.v,vv2.v,vv3.v);
- 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;
- double ptr3[3];
- ptr3[0] = vv3[0];
- ptr3[1] = vv3[1];
- ptr3[2] = 0;
-
- return astAngle(ast, ptr1, ptr2, ptr3);
- }
- 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;
- double ptr3[4];
- ptr3[0] = vv3[0];
- ptr3[1] = vv3[1];
- ptr3[2] = 0;
- ptr3[3] = 0;
-
- return astAngle(ast, ptr1, ptr2, ptr3);
- }
- }
-
- return 0;
-}
-
-double FitsImage::wcsAxAngle(AstFrameSet* ast,
- const Vector& vv1, const Vector& vv2)
-{
- int naxes = astGetI(ast,"Naxes");
- switch (naxes) {
- case 1:
- // error
- break;
- case 2:
- return astAxAngle(ast, vv1.v, vv2.v, 2);
- 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 astAxAngle(ast, ptr1, ptr2, 2);
- }
- 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 astAxAngle(ast, ptr1, ptr2, 2);
- }
- }
-
- return 0;
-}
-
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 e44ee31..9a521c0 100644
--- a/tksao/frame/fitsimage.h
+++ b/tksao/frame/fitsimage.h
@@ -10,7 +10,6 @@
#include "fitsdata.h"
#include "coord.h"
#include "file.h"
-#include "wcs.h"
#include "head.h"
// WCS ' ','A' to 'Z', WCS0
@@ -389,16 +388,6 @@ class FitsImage {
double getWCSDist(const Vector&, const Vector&, Coord::CoordSystem);
const char* getWCSName(Coord::CoordSystem);
- Vector wcsTran(AstFrameSet*, const Vector&, int);
- void wcsTran(AstFrameSet*, int, Vector*, int, Vector*);
- Vector3d wcsTran(AstFrameSet*, const Vector3d&, int);
- double wcsDistance(AstFrameSet*, const Vector&, const Vector&);
- double wcsAngle(AstFrameSet*, const Vector&, const Vector&, const Vector&);
- double wcsAxAngle(AstFrameSet*, const Vector&, const Vector&);
-
- int wcsSystem(AstFrameSet* ast, Coord::CoordSystem sys);
- int wcsSkyFrame(AstFrameSet* ast, Coord::SkyFrame sky);
-
void setWCSSkyFrame(Coord::CoordSystem, Coord::SkyFrame);
void setWCSFormat(int, const char*);
diff --git a/tksao/frame/grid25d.C b/tksao/frame/grid25d.C
index 17956e7..5d45d9a 100644
--- a/tksao/frame/grid25d.C
+++ b/tksao/frame/grid25d.C
@@ -6,10 +6,7 @@
#include "context.h"
#include "frame3dbase.h"
#include "fitsimage.h"
-
-extern "C" {
- #include "ast.h"
-}
+#include "wcsast.h"
extern Grid25dBase* astGrid25dPtr;
@@ -64,8 +61,8 @@ int Grid25d::doit(RenderMode rm)
}
AstFrameSet* ast = (AstFrameSet*)astCopy(fits->ast_);
- fits->wcsSystem(ast,system_);
- fits->wcsSkyFrame(ast,sky_);
+ wcsSystem(ast,system_);
+ wcsSkyFrame(ast,sky_);
int naxes = astGetI(ast,"Naxes");
switch (naxes) {
diff --git a/tksao/frame/grid2d.C b/tksao/frame/grid2d.C
index 02a80c7..38c09db 100644
--- a/tksao/frame/grid2d.C
+++ b/tksao/frame/grid2d.C
@@ -6,10 +6,7 @@
#include "context.h"
#include "framebase.h"
#include "fitsimage.h"
-
-extern "C" {
- #include "ast.h"
-}
+#include "wcsast.h"
extern Grid2dBase* astGrid2dPtr;
@@ -65,8 +62,8 @@ int Grid2d::doit(RenderMode rm)
}
AstFrameSet* ast = (AstFrameSet*)astCopy(fits->ast_);
- fits->wcsSystem(ast,system_);
- fits->wcsSkyFrame(ast,sky_);
+ wcsSystem(ast,system_);
+ wcsSkyFrame(ast,sky_);
int naxes = astGetI(ast,"Naxes");
switch (naxes) {
diff --git a/tksao/frame/grid3d.C b/tksao/frame/grid3d.C
index 8dad14b..6fa9b4e 100644
--- a/tksao/frame/grid3d.C
+++ b/tksao/frame/grid3d.C
@@ -6,10 +6,7 @@
#include "context.h"
#include "frame3dbase.h"
#include "fitsimage.h"
-
-extern "C" {
- #include "ast.h"
-}
+#include "wcsast.h"
extern Grid3dBase* astGrid3dPtr;
@@ -65,8 +62,8 @@ int Grid3d::doit(RenderMode rm)
}
AstFrameSet* ast = (AstFrameSet*)astCopy(fits->ast_);
- fits->wcsSystem(ast,system_);
- fits->wcsSkyFrame(ast,sky_);
+ wcsSystem(ast,system_);
+ wcsSkyFrame(ast,sky_);
int naxes = astGetI(ast,"Naxes");
switch (naxes) {
diff --git a/tksao/frame/wcsast.C b/tksao/frame/wcsast.C
new file mode 100644
index 0000000..0ee2540
--- /dev/null
+++ b/tksao/frame/wcsast.C
@@ -0,0 +1,384 @@
+// Copyright (C) 1999-2018
+// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+// For conditions of distribution and use, see copyright notice in "copyright"
+
+#include "wcsast.h"
+
+int wcsSystem(AstFrameSet* ast, Coord::CoordSystem sys)
+{
+ int nn = astGetI(ast,"nframe");
+ char cc = ' ';
+ int ww = sys-Coord::WCS;
+ if (ww < 0)
+ return 0;
+ if (ww)
+ cc = ww+'@';
+
+ for (int ss=0; ss<nn; ss++) {
+ const char* id = astGetC(astGetFrame(ast,ss+1),"Ident");
+ if (cc == id[0]) {
+ astSetI(ast,"Current",ss+1);
+ return 1;
+ }
+ }
+
+ return 0;
+}
+
+int wcsSkyFrame(AstFrameSet* ast, Coord::SkyFrame sky)
+{
+ AstFrameSet* fs =
+ (AstFrameSet*)astFindFrame(ast, astSkyFrame(" MaxAxes=4")," ");
+
+ if (!fs)
+ return 0;
+
+ switch (sky) {
+ case Coord::FK4:
+ astSet(fs, "System=FK4, Equinox=B1950");
+ break;
+ case Coord::FK5:
+ astSet(fs, "System=FK5, Equinox=J2000");
+ break;
+ case Coord::ICRS:
+ astSet(fs, "System=ICRS");
+ break;
+ case Coord::GALACTIC:
+ astSet(fs, "System=GALACTIC");
+ break;
+ case Coord::ECLIPTIC:
+ astSet(fs, "System=ECLIPTIC");
+ // get AST to agree with WCSSUBS
+ astSetD(fs, "EQUINOX", astGetD(fs, "EPOCH"));
+ break;
+ }
+
+ return 1;
+}
+
+Vector wcsTran(AstFrameSet* ast, const Vector& in, int forward)
+{
+ int naxes = astGetI(ast,"Naxes");
+ switch (naxes) {
+ case 1:
+ // error
+ break;
+ case 2:
+ double xout, yout;
+ astTran2(ast, 1, in.v, in.v+1, forward, &xout, &yout);
+ return Vector(xout, yout);
+ case 3:
+ {
+ double pin[3];
+ double pout[3];
+ pin[0] = in[0];
+ pin[1] = in[1];
+ pin[2] = forward ? 1 : 0;
+ astTranN(ast, 1, 3, 1, pin, forward, 3, 1, pout);
+ return Vector(pout[0],pout[1]);
+ }
+ break;
+ case 4:
+ {
+ double pin[4];
+ double pout[4];
+ pin[0] = in[0];
+ pin[1] = in[1];
+ pin[2] = forward ? 1 : 0;
+ pin[3] = forward ? 1 : 0;
+ astTranN(ast, 1, 4, 1, pin, forward, 4, 1, pout);
+ return Vector(pout[0],pout[1]);
+ }
+ break;
+ }
+ return Vector();
+}
+
+void wcsTran(AstFrameSet* ast, int npoint, Vector* in, int forward, Vector* out)
+{
+ int naxes = astGetI(ast,"Naxes");
+ switch (naxes) {
+ case 1:
+ // error
+ break;
+ case 2:
+ {
+ double* xin = new double[npoint];
+ double* yin = new double[npoint];
+ double* xout = new double[npoint];
+ double* yout = new double[npoint];
+
+ for (int ii=0; ii<npoint; ii++) {
+ xin[ii] = in[ii][0];
+ yin[ii] = in[ii][1];
+ }
+ astTran2(ast, npoint, xin, yin, forward, xout, yout);
+ for (int ii=0; ii<npoint; ii++)
+ out[ii] = Vector(xout[ii],yout[ii]);
+
+ if (xin)
+ delete [] xin;
+ if (yin)
+ delete [] yin;
+ if (xout)
+ delete [] xout;
+ if (yout)
+ delete [] yout;
+ }
+ break;
+ case 3:
+ {
+ double* ptr_in[3];
+ ptr_in[0] = new double[npoint];
+ ptr_in[1] = new double[npoint];
+ ptr_in[2] = new double[npoint];
+ double* ptr_out[3];
+ ptr_out[0] = new double[npoint];
+ ptr_out[1] = new double[npoint];
+ ptr_out[2] = new double[npoint];
+
+ for (int kk=0; kk<npoint; kk++) {
+ ptr_in[0][kk] = in[kk][0];
+ ptr_in[1][kk] = in[kk][1];
+ ptr_in[2][kk] = forward ? 1 : 0;
+ }
+ astTranP(ast, npoint, 3, (const double**)ptr_in, forward, 3, ptr_out);
+ for (int kk=0; kk<npoint; kk++)
+ out[kk] = Vector(ptr_out[0][kk], ptr_out[1][kk]);
+
+ if (ptr_in[0])
+ delete [] ptr_in[0];
+ if (ptr_in[1])
+ delete [] ptr_in[1];
+ if (ptr_in[2])
+ delete [] ptr_in[2];
+
+ if (ptr_out[0])
+ delete [] ptr_out[0];
+ if (ptr_out[1])
+ delete [] ptr_out[1];
+ if (ptr_out[2])
+ delete [] ptr_out[2];
+ }
+ break;
+ case 4:
+ {
+ double* ptr_in[4];
+ ptr_in[0] = new double[npoint];
+ ptr_in[1] = new double[npoint];
+ ptr_in[2] = new double[npoint];
+ ptr_in[3] = new double[npoint];
+ double* ptr_out[4];
+ ptr_out[0] = new double[npoint];
+ ptr_out[1] = new double[npoint];
+ ptr_out[2] = new double[npoint];
+ ptr_out[3] = new double[npoint];
+
+ for (int kk=0; kk<npoint; kk++) {
+ ptr_in[0][kk] = in[kk][0];
+ ptr_in[1][kk] = in[kk][1];
+ ptr_in[2][kk] = forward ? 1 : 0;
+ ptr_in[3][kk] = forward ? 1 : 0;
+ }
+ astTranP(ast, npoint, 4, (const double**)ptr_in, forward, 4, ptr_out);
+ for (int kk=0; kk<npoint; kk++)
+ out[kk] = Vector(ptr_out[0][kk], ptr_out[1][kk]);
+
+ if (ptr_in[0])
+ delete [] ptr_in[0];
+ if (ptr_in[1])
+ delete [] ptr_in[1];
+ if (ptr_in[2])
+ delete [] ptr_in[2];
+ if (ptr_in[3])
+ delete [] ptr_in[3];
+
+ if (ptr_out[0])
+ delete [] ptr_out[0];
+ if (ptr_out[1])
+ delete [] ptr_out[1];
+ if (ptr_out[2])
+ delete [] ptr_out[2];
+ if (ptr_out[3])
+ delete [] ptr_out[3];
+ }
+ break;
+ }
+}
+
+Vector3d wcsTran(AstFrameSet* ast, const Vector3d& in, int forward)
+{
+ int naxes = astGetI(ast,"Naxes");
+ switch (naxes) {
+ case 1:
+ case 2:
+ double pin[2];
+ double pout[2];
+ pin[0] = in[0];
+ pin[1] = in[1];
+ astTranN(ast, 1, 2, 1, pin, forward, 2, 1, pout);
+ return Vector3d(pout[0],pout[1],forward ? 1 : 0);
+ break;
+ case 3:
+ {
+ double pin[3];
+ double pout[3];
+ pin[0] = in[0];
+ pin[1] = in[1];
+ pin[2] = in[2];
+ astTranN(ast, 1, 3, 1, pin, forward, 3, 1, pout);
+ return Vector3d(pout[0],pout[1],pout[2]);
+ }
+ break;
+ case 4:
+ {
+ double pin[4];
+ double pout[4];
+ pin[0] = in[0];
+ pin[1] = in[1];
+ pin[2] = in[2];
+ pin[3] = forward ? 1 : 0;
+ astTranN(ast, 1, 4, 1, pin, forward, 4, 1, pout);
+ return Vector3d(pout[0],pout[1],pout[2]);
+ }
+ break;
+ }
+ return Vector3d();
+}
+
+double wcsDistance(AstFrameSet* ast, const Vector& vv1, const Vector& vv2)
+{
+ int naxes = astGetI(ast,"Naxes");
+ switch (naxes) {
+ case 1:
+ // error
+ break;
+ case 2:
+ return astDistance(ast, vv1.v, vv2.v);
+ 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);
+ }
+ 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);
+ }
+ }
+
+ return 0;
+}
+
+double wcsAngle(AstFrameSet* ast, const Vector& vv1, const Vector& vv2,
+ const Vector& vv3)
+{
+ int naxes = astGetI(ast,"Naxes");
+ switch (naxes) {
+ case 1:
+ // error
+ break;
+ case 2:
+ return astAngle(ast,vv1.v,vv2.v,vv3.v);
+ 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;
+ double ptr3[3];
+ ptr3[0] = vv3[0];
+ ptr3[1] = vv3[1];
+ ptr3[2] = 0;
+
+ return astAngle(ast, ptr1, ptr2, ptr3);
+ }
+ 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;
+ double ptr3[4];
+ ptr3[0] = vv3[0];
+ ptr3[1] = vv3[1];
+ ptr3[2] = 0;
+ ptr3[3] = 0;
+
+ return astAngle(ast, ptr1, ptr2, ptr3);
+ }
+ }
+
+ return 0;
+}
+
+double wcsAxAngle(AstFrameSet* ast, const Vector& vv1, const Vector& vv2)
+{
+ int naxes = astGetI(ast,"Naxes");
+ switch (naxes) {
+ case 1:
+ // error
+ break;
+ case 2:
+ return astAxAngle(ast, vv1.v, vv2.v, 2);
+ 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 astAxAngle(ast, ptr1, ptr2, 2);
+ }
+ 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 astAxAngle(ast, ptr1, ptr2, 2);
+ }
+ }
+
+ return 0;
+}
+
diff --git a/tksao/frame/wcsast.h b/tksao/frame/wcsast.h
new file mode 100644
index 0000000..258c198
--- /dev/null
+++ b/tksao/frame/wcsast.h
@@ -0,0 +1,22 @@
+// Copyright (C) 1999-2018
+// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+// For conditions of distribution and use, see copyright notice in "copyright"
+
+#include "vector.h"
+#include "vector3d.h"
+#include "coord.h"
+
+extern "C" {
+#include "ast.h"
+}
+
+int wcsSystem(AstFrameSet* ast, Coord::CoordSystem sys);
+int wcsSkyFrame(AstFrameSet* ast, Coord::SkyFrame sky);
+
+Vector wcsTran(AstFrameSet*, const Vector&, int);
+void wcsTran(AstFrameSet*, int, Vector*, int, Vector*);
+Vector3d wcsTran(AstFrameSet*, const Vector3d&, int);
+
+double wcsDistance(AstFrameSet*, const Vector&, const Vector&);
+double wcsAngle(AstFrameSet*, const Vector&, const Vector&, const Vector&);
+double wcsAxAngle(AstFrameSet*, const Vector&, const Vector&);