diff options
-rwxr-xr-x | tksao/configure | 1 | ||||
-rw-r--r-- | tksao/configure.ac | 1 | ||||
-rw-r--r-- | tksao/frame/base.C | 9 | ||||
-rw-r--r-- | tksao/frame/base.h | 6 | ||||
-rw-r--r-- | tksao/frame/contourparser.C | 4 | ||||
-rw-r--r-- | tksao/frame/contourparser.Y | 4 | ||||
-rw-r--r-- | tksao/frame/ds9parser.C | 4 | ||||
-rw-r--r-- | tksao/frame/ds9parser.Y | 4 | ||||
-rw-r--r-- | tksao/frame/fitsimage.C | 383 | ||||
-rw-r--r-- | tksao/frame/fitsimage.h | 11 | ||||
-rw-r--r-- | tksao/frame/grid25d.C | 9 | ||||
-rw-r--r-- | tksao/frame/grid2d.C | 9 | ||||
-rw-r--r-- | tksao/frame/grid3d.C | 9 | ||||
-rw-r--r-- | tksao/frame/wcsast.C | 384 | ||||
-rw-r--r-- | tksao/frame/wcsast.h | 22 |
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&); |