upgrade ast 8.7.1

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;
+}
+