summaryrefslogtreecommitdiffstats
path: root/ast/tpn.c
diff options
context:
space:
mode:
Diffstat (limited to 'ast/tpn.c')
-rw-r--r--ast/tpn.c393
1 files changed, 393 insertions, 0 deletions
diff --git a/ast/tpn.c b/ast/tpn.c
new file mode 100644
index 0000000..0d30503
--- /dev/null
+++ b/ast/tpn.c
@@ -0,0 +1,393 @@
+#include <math.h>
+#include <stdlib.h>
+#include "wcsmath.h"
+#include "wcstrig.h"
+#include "proj.h"
+
+#define icopysign(X, Y) ((Y) < 0.0 ? -abs(X) : abs(X))
+#define TPN 999
+
+/*============================================================================
+* TAN: gnomonic projection, with correction terms.
+*
+* This projection is no longer part of the FITSWCS standard, but is
+* retained here for use b the AST library as a means of implementing
+* the IRAF TNX projection, and the DSS encoding.
+*
+* Given and/or returned:
+* prj->p Array of latitude coefficients
+* prj->p2 Array of longitude coefficients
+* prj->flag TPN, or -TPN if prj->flag is given < 0.
+* prj->r0 r0; reset to 180/pi if 0.
+* prj->n If zero, only do poly part of transformation (i.e. omit
+* the TAN projection).
+*
+* Returned:
+* prj->code "TPN"
+* prj->phi0 0.0
+* prj->theta0 90.0
+* prj->astPRJfwd Pointer to astTPNfwd().
+* prj->astPRJrev Pointer to astTPNrev().
+* prj->w[ 0 ] Set to 0.0 if a simple tan projection is required
+* (with no polynomial correction). Otherwise, set to 1.0.
+*===========================================================================*/
+
+int astTPNset(prj)
+
+struct AstPrjPrm *prj;
+
+{
+ int m;
+
+ prj->flag = icopysign(TPN, prj->flag);
+ prj->phi0 = 0.0;
+ prj->theta0 = 90.0;
+
+ if (prj->r0 == 0.0) prj->r0 = R2D;
+
+ prj->astPRJfwd = astTPNfwd;
+ prj->astPRJrev = astTPNrev;
+
+/* If all co-efficients have their "unit" values, we do not need to
+ use the polynomial correction. */
+ prj->w[ 0 ] = 0.0;
+
+ if( prj->p[ 0 ] != 0.0 || prj->p2[ 0 ] != 0.0 ) {
+ prj->w[ 0 ] = 1.0;
+
+ } else if( prj->p[ 1 ] != 1.0 || prj->p2[ 1 ] != 1.0 ) {
+ prj->w[ 0 ] = 1.0;
+
+ } else {
+ for( m = 2; m < WCSLIB_MXPAR; m++ ){
+ if( prj->p[ m ] != 0.0 || prj->p2[ m ] != 0.0 ){
+ prj->w[ 0 ] = 1.0;
+ break;
+ }
+ }
+ }
+
+ return 0;
+}
+
+/*--------------------------------------------------------------------------*/
+
+int astTPNfwd(phi, theta, prj, xx, yy)
+
+const double phi, theta;
+double *xx, *yy;
+struct AstPrjPrm *prj;
+
+{
+ double r, xi, eta, x2, xy, y2, r2, x3, x2y, xy2, y3, r3, x4, x3y,
+ x2y2, xy3, y4, x5, x4y, x3y2, x2y3, xy4, y5, r5, x6, x5y, x4y2,
+ x3y3, x2y4, xy5, y6, x7, x6y, x5y2, x4y3, x3y4, x2y5, xy6, y7,
+ r7, tol, f, g, fx, fy, gx, gy, dx, dy, x, y, denom;
+ double *a, *b;
+ int i, ok;
+
+ if (abs(prj->flag) != TPN ) {
+ if (astTPNset(prj)) return 1;
+ }
+
+ if( prj->n ) {
+ double s = astSind(theta);
+ if (prj->flag > 0 && s < 0.0) {
+ return 2;
+ }
+ r = prj->r0*astCosd(theta)/s;
+ xi = r*astSind(phi);
+ eta = -r*astCosd(phi);
+ } else {
+ xi = phi;
+ eta = theta;
+ }
+
+ /* Simple tan */
+ if( prj->w[ 0 ] == 0.0 ){
+ *xx = xi;
+ *yy = eta;
+
+ /* Tan with polynomial corrections: Iterate using Newton's method to
+ get the (x,y) corresponding to the above (xi,eta). */
+ } else {
+ a = prj->p2;
+ b = prj->p;
+
+ /* Initial guess: linear solution assuming a3,... and b3,... are zero. */
+ denom = a[1]*b[1] - a[2]*b[2];
+ if( denom != 0.0 ) {
+ x = ( xi*b[1] - eta*a[2] - a[0]*b[1] + b[0]*a[2] )/denom;
+ y = -( xi*b[2] - eta*a[1] - a[0]*b[2] + b[0]*a[1] )/denom;
+
+ } else {
+ if( a[1] != 0.0 ){
+ x = ( xi - a[0] )/a[1];
+ } else {
+ x = a[0];
+ }
+
+ if( b[1] != 0.0 ){
+ y = ( eta - b[0] )/b[1];
+ } else {
+ y = b[0];
+ }
+ }
+
+/* Iterate up to 50 times, until the required relative accuracy is
+ achieved. */
+ tol = 1.0E-5;
+ ok = 0;
+ for (i = 0; i < 50; i++) {
+
+/* Get required products of the current x and y values */
+ x2 = x*x;
+ xy = x*y;
+ y2 = y*y;
+
+ r2 = x2 + y2;
+ r = sqrt( r2 );
+
+ x3 = x2*x;
+ x2y = x2*y;
+ xy2 = x*y2;
+ y3 = y*y2;
+
+ r3 = r*r2;
+
+ x4 = x3*x;
+ x3y = x3*y;
+ x2y2 = x2*y2;
+ xy3 = x*y3;
+ y4 = y*y3;
+
+ x5 = x4*x;
+ x4y = x4*y;
+ x3y2 = x3*y2;
+ x2y3 = x2*y3;
+ xy4 = x*y4;
+ y5 = y*y4;
+
+ r5 = r3*r2;
+
+ x6 = x5*x;
+ x5y = x5*y;
+ x4y2 = x4*y2;
+ x3y3 = x3*y3;
+ x2y4 = x2*y4;
+ xy5 = x*y5;
+ y6 = y*y5;
+
+ x7 = x6*x;
+ x6y = x6*y;
+ x5y2 = x5*y2;
+ x4y3 = x4*y3;
+ x3y4 = x3*y4;
+ x2y5 = x2*y5;
+ xy6 = x*y6;
+ y7 = y*y6;
+
+ r7 = r5*r2;
+
+/* Get the xi and eta models corresponding to the current x and y values */
+ f = a[0] + a[1]*x + a[2]*y + a[3]*r + a[4]*x2
+ + a[5]*xy + a[6]*y2 + a[7]*x3 + a[8]*x2y + a[9]*xy2
+ + a[10]*y3 + a[11]*r3 + a[12]*x4 + a[13]*x3y + a[14]*x2y2
+ + a[15]*xy3 + a[16]*y4 + a[17]*x5 + a[18]*x4y + a[19]*x3y2
+ + a[20]*x2y3 + a[21]*xy4 + a[22]*y5 + a[23]*r5 + a[24]*x6
+ + a[25]*x5y + a[26]*x4y2 + a[27]*x3y3 + a[28]*x2y4 + a[29]*xy5
+ + a[30]*y6 + a[31]*x7 + a[32]*x6y + a[33]*x5y2 + a[34]*x4y3
+ + a[35]*x3y4 + a[36]*x2y5 + a[37]*xy6 + a[38]*y7 + a[39]*r7;
+
+ g = b[0] + b[1]*y + b[2]*x + b[3]*r + b[4]*y2
+ + b[5]*xy + b[6]*x2 + b[7]*y3 + b[8]*xy2 + b[9]*x2y
+ + b[10]*x3 + b[11]*r3 + b[12]*y4 + b[13]*xy3 + b[14]*x2y2
+ + b[15]*x3y + b[16]*x4 + b[17]*y5 + b[18]*xy4 + b[19]*x2y3
+ + b[20]*x3y2 + b[21]*x4y + b[22]*x5 + b[23]*r5 + b[24]*y6
+ + b[25]*xy5 + b[26]*x2y4 + b[27]*x3y3 + b[28]*x4y2 + b[29]*x5y
+ + b[30]*x6 + b[31]*y7 + b[32]*xy6 + b[33]*x2y5 + b[34]*x3y4
+ + b[35]*x4y3 + b[36]*x5y2 + b[37]*x6y + b[38]*x7 + b[39]*r7;
+
+/* Partial derivative of xi wrt x... */
+ fx = a[1] + a[3]*( (r!=0.0)?(x/r):0.0 ) + 2*a[4]*x +
+ a[5]*y + 3*a[7]*x2 + 2*a[8]*xy + a[9]*y2 +
+ 3*a[11]*r*x + 4*a[12]*x3 + 3*a[13]*x2y + 2*a[14]*xy2 +
+ a[15]*y3 + 5*a[17]*x4 + 4*a[18]*x3y + 3*a[19]*x2y2 +
+ 2*a[20]*xy3 + a[21]*y4 + 5*a[23]*r3*x + 6*a[24]*x5 +
+ 5*a[25]*x4y + 4*a[26]*x3y2 + 3*a[27]*x2y3 + 2*a[28]*xy4 +
+ a[29]*y5 + 7*a[31]*x6 + 6*a[32]*x5y + 5*a[33]*x4y2 +
+ 4*a[34]*x3y3 + 3*a[35]*x2y4 + 2*a[36]*xy5 + a[37]*y6 +
+ 7*a[39]*r5*x;
+
+/* Partial derivative of xi wrt y... */
+ fy = a[2] + a[3]*( (r!=0.0)?(y/r):0.0 ) + a[5]*x +
+ 2*a[6]*y + a[8]*x2 + 2*a[9]*xy + 3*a[10]*y2 +
+ 3*a[11]*r*y + a[13]*x3 + 2*a[14]*x2y + 3*a[15]*xy2 +
+ 4*a[16]*y3 + a[18]*x4 + 2*a[19]*x3y + 3*a[20]*x2y2 +
+ 4*a[21]*xy3 + 5*a[22]*y4 + 5*a[23]*r3*y + a[25]*x5 +
+ 2*a[26]*x4y + 3*a[27]*x3y2 + 4*a[28]*x2y3 + 5*a[29]*xy4 +
+ 6*a[30]*y5 + a[32]*x6 + 2*a[33]*x5y + 3*a[34]*x4y2 +
+ 4*a[35]*x3y3 + 5*a[36]*x2y4 + 6*a[37]*xy5 + 7*a[38]*y6 +
+ 7*a[39]*r5*y;
+
+/* Partial derivative of eta wrt x... */
+ gx = b[2] + b[3]*( (r!=0.0)?(x/r):0.0 ) + b[5]*y +
+ 2*b[6]*x + b[8]*y2 + 2*b[9]*xy + 3*b[10]*x2 +
+ 3*b[11]*r*x + b[13]*y3 + 2*b[14]*xy2 + 3*b[15]*x2y +
+ 4*b[16]*x3 + b[18]*y4 + 2*b[19]*xy3 + 3*b[20]*x2y2 +
+ 4*b[21]*x3y + 5*b[22]*x4 + 5*b[23]*r3*x + b[25]*y5 +
+ 2*b[26]*xy4 + 3*b[27]*x2y3 + 4*b[28]*x3y2 + 5*b[29]*x4y +
+ 6*b[30]*x5 + b[32]*y6 + 2*b[33]*xy5 + 3*b[34]*x2y4 +
+ 4*b[35]*x3y3 + 5*b[36]*x4y2 + 6*b[37]*x5y + 7*b[38]*x6 +
+ 7*b[39]*r5*x;
+
+/* Partial derivative of eta wrt y... */
+ gy = b[1] + b[3]*( (r!=0.0)?(y/r):0.0 ) + 2*b[4]*y +
+ b[5]*x + 3*b[7]*y2 + 2*b[8]*xy + b[9]*x2 +
+ 3*b[11]*r*y + 4*b[12]*y3 + 3*b[13]*xy2 + 2*b[14]*x2y +
+ b[15]*x3 + 5*b[17]*y4 + 4*b[18]*xy3 + 3*b[19]*x2y2 +
+ 2*b[20]*x3y + b[21]*x4 + 5*b[23]*r3*y + 6*b[24]*y5 +
+ 5*b[25]*xy4 + 4*b[26]*x2y3 + 3*b[27]*x3y2 + 2*b[28]*x4y +
+ b[29]*x5 + 7*b[31]*y6 + 6*b[32]*xy5 + 5*b[33]*x2y4 +
+ 4*b[34]*x3y3 + 3*b[35]*x4y2 + 2*b[36]*x5y + b[37]*x6 +
+ 7*b[39]*r5*y;
+
+/* Calculate new x and y values. */
+ f = f - xi;
+ g = g - eta;
+ dx = ( (-f*gy) + (g*fy) ) / ( (fx*gy) - (fy*gx) );
+ dy = ( (-g*fx) + (f*gx) ) / ( (fx*gy) - (fy*gx) );
+ x += dx;
+ y += dy;
+
+/* Check if convergence has been achieved. */
+ if( fabs(dx) <= tol*fabs(x) && fabs(dy) <= tol*fabs(y) ) {
+ ok = 1;
+ break;
+ }
+ }
+
+ *xx = x;
+ *yy = y;
+
+ if( !ok ) return 2;
+
+ }
+
+ return 0;
+
+}
+
+/*--------------------------------------------------------------------------*/
+
+int astTPNrev(x, y, prj, phi, theta)
+
+const double x, y;
+double *phi, *theta;
+struct AstPrjPrm *prj;
+
+{
+ double r, xi, eta, x2, xy, y2, r2, x3, x2y, xy2, y3, r3, x4, x3y,
+ x2y2, xy3, y4, x5, x4y, x3y2, x2y3, xy4, y5, r5, x6, x5y, x4y2,
+ x3y3, x2y4, xy5, y6, x7, x6y, x5y2, x4y3, x3y4, x2y5, xy6, y7,
+ r7;
+ double *a, *b;
+
+ if (abs(prj->flag) != TPN ) {
+ if (astTPNset(prj)) return 1;
+ }
+
+ /* Simple tan */
+ if( prj->w[ 0 ] == 0.0 ){
+ xi = x;
+ eta = y;
+
+ /* Tan with polynomial corrections. */
+ } else {
+ x2 = x*x;
+ xy = x*y;
+ y2 = y*y;
+
+ r2 = x2 + y2;
+ r = sqrt( r2 );
+
+ x3 = x2*x;
+ x2y = x2*y;
+ xy2 = x*y2;
+ y3 = y*y2;
+
+ r3 = r*r2;
+
+ x4 = x3*x;
+ x3y = x3*y;
+ x2y2 = x2*y2;
+ xy3 = x*y3;
+ y4 = y*y3;
+
+ x5 = x4*x;
+ x4y = x4*y;
+ x3y2 = x3*y2;
+ x2y3 = x2*y3;
+ xy4 = x*y4;
+ y5 = y*y4;
+
+ r5 = r3*r2;
+
+ x6 = x5*x;
+ x5y = x5*y;
+ x4y2 = x4*y2;
+ x3y3 = x3*y3;
+ x2y4 = x2*y4;
+ xy5 = x*y5;
+ y6 = y*y5;
+
+ x7 = x6*x;
+ x6y = x6*y;
+ x5y2 = x5*y2;
+ x4y3 = x4*y3;
+ x3y4 = x3*y4;
+ x2y5 = x2*y5;
+ xy6 = x*y6;
+ y7 = y*y6;
+
+ r7 = r5*r2;
+
+ a = prj->p2;
+ xi = a[0] + a[1]*x + a[2]*y + a[3]*r + a[4]*x2
+ + a[5]*xy + a[6]*y2 + a[7]*x3 + a[8]*x2y + a[9]*xy2
+ + a[10]*y3 + a[11]*r3 + a[12]*x4 + a[13]*x3y + a[14]*x2y2
+ + a[15]*xy3 + a[16]*y4 + a[17]*x5 + a[18]*x4y + a[19]*x3y2
+ + a[20]*x2y3 + a[21]*xy4 + a[22]*y5 + a[23]*r5 + a[24]*x6
+ + a[25]*x5y + a[26]*x4y2 + a[27]*x3y3 + a[28]*x2y4 + a[29]*xy5
+ + a[30]*y6 + a[31]*x7 + a[32]*x6y + a[33]*x5y2 + a[34]*x4y3
+ + a[35]*x3y4 + a[36]*x2y5 + a[37]*xy6 + a[38]*y7 + a[39]*r7;
+
+ b = prj->p;
+ eta = b[0] + b[1]*y + b[2]*x + b[3]*r + b[4]*y2
+ + b[5]*xy + b[6]*x2 + b[7]*y3 + b[8]*xy2 + b[9]*x2y
+ + b[10]*x3 + b[11]*r3 + b[12]*y4 + b[13]*xy3 + b[14]*x2y2
+ + b[15]*x3y + b[16]*x4 + b[17]*y5 + b[18]*xy4 + b[19]*x2y3
+ + b[20]*x3y2 + b[21]*x4y + b[22]*x5 + b[23]*r5 + b[24]*y6
+ + b[25]*xy5 + b[26]*x2y4 + b[27]*x3y3 + b[28]*x4y2 + b[29]*x5y
+ + b[30]*x6 + b[31]*y7 + b[32]*xy6 + b[33]*x2y5 + b[34]*x3y4
+ + b[35]*x4y3 + b[36]*x5y2 + b[37]*x6y + b[38]*x7 + b[39]*r7;
+
+ }
+
+ /* Now do the tan projection */
+ if( prj->n ) {
+ r = sqrt(xi*xi + eta*eta);
+ if (r == 0.0) {
+ *phi = 0.0;
+ } else {
+ *phi = astATan2d(xi, -eta);
+ }
+ *theta = astATan2d(prj->r0, r);
+ } else {
+ *phi = xi;
+ *theta = eta;
+ }
+
+ return 0;
+}
+