upgrade ast 8.7.1

1 files changed, 0 insertions, 6359 deletions

diff --git a/ast/polymap.c b/ast/polymap.c
deleted file mode 100644
index f03a631..0000000
--- a/ast/polymap.c
+++ /dev/null
@@ -1,6359 +0,0 @@
-/*
-*class++
-*  Name:
-*     PolyMap
-
-*  Purpose:
-*     Map coordinates using polynomial functions.
-
-*  Constructor Function:
-c     astPolyMap
-f     AST_POLYMAP
-
-*  Description:
-*     A PolyMap is a form of Mapping which performs a general polynomial
-*     transformation.  Each output coordinate is a polynomial function of
-*     all the input coordinates. The coefficients are specified separately
-*     for each output coordinate. The forward and inverse transformations
-*     are defined independantly by separate sets of coefficients. If no
-*     inverse transformation is supplied, the default behaviour is to use
-*     an iterative method to evaluate the inverse based only on the forward
-*     transformation (see attribute IterInverse).
-
-*  Inheritance:
-*     The PolyMap class inherits from the Mapping class.
-
-*  Attributes:
-*     In addition to those attributes common to all Mappings, every
-*     PolyMap also has the following attributes:
-*
-*     - IterInverse: Provide an iterative inverse transformation?
-*     - NiterInverse: Maximum number of iterations for iterative inverse
-*     - TolInverse: Target relative error for iterative inverse
-
-*  Functions:
-c     In addition to those functions applicable to all Objects, the
-c     following functions may also be applied to all Mappings:
-f     In addition to those routines applicable to all Objects, the
-f     following routines may also be applied to all Mappings:
-*
-c     - astPolyCoeffs: Retrieve the coefficients of a PolyMap transformation
-c     - astPolyTran: Fit a PolyMap inverse or forward transformation
-f     - AST_POLYCOEFFS: Retrieve the coefficients of a PolyMap transformation
-f     - AST_POLYTRAN: Fit a PolyMap inverse or forward transformation
-
-*  Copyright:
-*     Copyright (C) 2017 East Asian Observatory.
-*     Copyright (C) 1997-2006 Council for the Central Laboratory of the
-*     Research Councils
-*     Copyright (C) 2011 Science & Technology Facilities Council.
-*     All Rights Reserved.
-
-*  Licence:
-*     This program is free software: you can redistribute it and/or
-*     modify it under the terms of the GNU Lesser General Public
-*     License as published by the Free Software Foundation, either
-*     version 3 of the License, or (at your option) any later
-*     version.
-*
-*     This program is distributed in the hope that it will be useful,
-*     but WITHOUT ANY WARRANTY; without even the implied warranty of
-*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-*     GNU Lesser General Public License for more details.
-*
-*     You should have received a copy of the GNU Lesser General
-*     License along with this program.  If not, see
-*     <http://www.gnu.org/licenses/>.
-
-*  Authors:
-*     DSB: D.S. Berry (Starlink)
-
-*  History:
-*     27-SEP-2003 (DSB):
-*        Original version.
-*     13-APR-2005 (DSB):
-*        Changed the keys used by the Dump/astLoadPolyMap functions. They
-*        used to exceed 8 characters and consequently caused problems for
-*        FitsChans.
-*     20-MAY-2005 (DSB):
-*        Correct the indexing of keywords produced in the Dump function.
-*     20-APR-2006 (DSB):
-*        Guard against undefined transformations in Copy.
-*     10-MAY-2006 (DSB):
-*        Override astEqual.
-*     4-JUL-2008 (DSB):
-*        Fixed loop indexing problems in Equal function.
-*     27-MAY-2011 (DSB):
-*        Added public method astPolyTran.
-*     18-JUL-2011 (DSB):
-*        - Added attributes IterInverse, NiterInverse and TolInverse.
-*        - Do not report an error if astPolyTran fails to fit an inverse.
-*     15-OCT-2011 (DSB):
-*        Improve argument checking and error reporting in PolyTran
-*     8-MAY-2014 (DSB):
-*        Move to using CMinPack for minimisations.
-*     11-NOV-2016 (DSB):
-*        - Fix bug in MapMerge that could cause a seg fault. It did not
-*        check that the PolyMap had a defined transformation before accessing
-*        the transformation's coefficient array.
-*        - Fix similar bugs in Equal that could cause seg faults.
-*     15-MAR-2017 (DSB):
-*        - Change the GetTranForward and GetTranInverse functions so that they
-*        take into account the state of the Invert attribute.
-*        - Improve docs for the IterInverse attribute to explain that the
-*        inverse transformation replaced is always the original inverse
-*        transformation, as defined by the arguments supplied to the PolyMap
-*        constructor, regardless of the state of the Invert attribute.
-*     17-MAR-2017 (DSB):
-*        - Add the astPolyCoeffs method.
-*     30-MAR-2017 (DSB):
-*        Modify the astPolyTran method so that it can be used by the
-*        ChebyMap class to determine new transformations implemented as
-*        Chebyshev polynomials.
-*     27-JUN-2017 (DSB):
-*        In SamplePoly1D/2D ensure the final sample on each axis does not
-*        go above the supplied upper bound. This can happen due to rounding
-*        error. This is important for ChebyMaps since points outside the
-*        bounds are set bad when transformed using a ChebyMap, causing NaNs
-*        to be generated in lmder1 (cminpack minimisation function).
-*     3-JUL-2017 (DSB):
-*        Within FitPoly1D and FitPoly2D, use an initial guess that represents
-*        a unit mapping between normalised input and output values, rather
-*        than unnormalised PolyMap values. This is in case the PolyMap being
-*        fitted includes a change of scale (e.g. the PolyMap input is in "mm"
-*        but the output is in "rads" and includes some large scaling factor
-*        to do the conversion).
-*     9-JAN-2018 (DSB):
-*        Correct Transform to take account of AST__BAD coeffs correctly.
-*        Previously bad coeffs could generate NaN output values.
-*     27-APR-2018 (DSB):
-*        When calculating the iterative inverse, use an initial guess based
-*        on the linear truncation of the PolyMap rather than a UnitMap.
-*     3-MAY-2018 (DSB):
-*        Clear the IterInverse attribute in the PolyMap returned by
-*        astPolyTran (why would you want to use the iterative inverse if
-*        you have just created an inverse fit?).
-*     4-MAY-2018 (DSB):
-*        Fix various places where there was confusion over whether the
-*        "forward transformation" refers to the forward transformation
-*        of the original uninverted PolyMap, or the current forward
-*        transformation of the PolyMap (i.e. taking the "Invert" flag into
-*        account).
-*class--
-*/
-
-/* Module Macros. */
-/* ============== */
-/* Set the name of the class we are implementing. This indicates to
-   the header files that define class interfaces that they should make
-   "protected" symbols available. */
-#define astCLASS PolyMap
-
-/* Include files. */
-/* ============== */
-/* Interface definitions. */
-/* ---------------------- */
-
-#include "globals.h"             /* Thread-safe global data access */
-#include "error.h"               /* Error reporting facilities */
-#include "memory.h"              /* Memory allocation facilities */
-#include "object.h"              /* Base Object class */
-#include "pointset.h"            /* Sets of points/coordinates */
-#include "matrixmap.h"           /* Matrix multiplication mappings */
-#include "shiftmap.h"            /* Shift of origin mappings */
-#include "mapping.h"             /* Coordinate mappings (parent class) */
-#include "cmpmap.h"              /* Compound mappings */
-#include "polymap.h"             /* Interface definition for this class */
-#include "unitmap.h"             /* Unit mappings */
-#include "cminpack/cminpack.h"   /* Levenberg - Marquardt minimization */
-#include "pal.h"                 /* SLALIB function definitions */
-
-/* Error code definitions. */
-/* ----------------------- */
-#include "ast_err.h"             /* AST error codes */
-
-/* C header files. */
-/* --------------- */
-#include <ctype.h>
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <limits.h>
-#include <float.h>
-
-/* Module Variables. */
-/* ================= */
-
-/* Address of this static variable is used as a unique identifier for
-   member of this class. */
-static int class_check;
-
-/* Pointers to parent class methods which are extended by this class. */
-static AstPointSet *(* parent_transform)( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
-static const char *(* parent_getattrib)( AstObject *, const char *, int * );
-static int (* parent_testattrib)( AstObject *, const char *, int * );
-static void (* parent_clearattrib)( AstObject *, const char *, int * );
-static void (* parent_setattrib)( AstObject *, const char *, int * );
-static int (* parent_getobjsize)( AstObject *, int * );
-
-#if defined(THREAD_SAFE)
-static int (* parent_managelock)( AstObject *, int, int, AstObject **, int * );
-#endif
-
-
-#ifdef THREAD_SAFE
-/* Define how to initialise thread-specific globals. */
-#define GLOBAL_inits \
-   globals->Class_Init = 0; \
-   globals->GetAttrib_Buff[ 0 ] = 0;
-
-/* Create the function that initialises global data for this module. */
-astMAKE_INITGLOBALS(PolyMap)
-
-/* Define macros for accessing each item of thread specific global data. */
-#define class_init astGLOBAL(PolyMap,Class_Init)
-#define class_vtab astGLOBAL(PolyMap,Class_Vtab)
-#define getattrib_buff astGLOBAL(LutMap,GetAttrib_Buff)
-
-#include <pthread.h>
-
-
-#else
-
-static char getattrib_buff[ 101 ];
-
-/* Define the class virtual function table and its initialisation flag
-   as static variables. */
-static AstPolyMapVtab class_vtab;   /* Virtual function table */
-static int class_init = 0;       /* Virtual function table initialised? */
-
-#endif
-
-
-
-/* External Interface Function Prototypes. */
-/* ======================================= */
-/* The following functions have public prototypes only (i.e. no
-   protected prototypes), so we must provide local prototypes for use
-   within this module. */
-AstPolyMap *astPolyMapId_( int, int, int, const double[], int, const double[], const char *, ... );
-
-
-/* Prototypes for Private Member Functions. */
-/* ======================================== */
-static AstMapping *LinearGuess( AstPolyMap *, int * );
-static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
-static AstPolyMap **GetJacobian( AstPolyMap *, int * );
-static AstPolyMap *PolyTran( AstPolyMap *, int, double, double, int, const double *, const double *, int * );
-static double **SamplePoly1D( AstPolyMap *, int, double **, double, double, int, int *, double[2], int * );
-static double **SamplePoly2D( AstPolyMap *, int, double **, const double *, const double *, int, int *, double[4], int * );
-static double *FitPoly1D( AstPolyMap *, int, int, double, int, double **, double[2], int *, double *, int * );
-static double *FitPoly2D( AstPolyMap *, int, int, double, int, double **, double[4], int *, double *, int * );
-static int Equal( AstObject *, AstObject *, int * );
-static int GetObjSize( AstObject *, int * );
-static int GetTranForward( AstMapping *, int * );
-static int GetTranInverse( AstMapping *, int * );
-static int MPFunc1D( void *, int, int, const double *, double *, double *, int, int );
-static int MPFunc2D( void *, int, int, const double *, double *, double *, int, int );
-static int MapMerge( AstMapping *, int, int, int *, AstMapping ***, int **, int * );
-static int ReplaceTransformation( AstPolyMap *, int, double, double, int, const double *, const double *, int * );
-static void Copy( const AstObject *, AstObject *, int * );
-static void Delete( AstObject *obj, int * );
-static void Dump( AstObject *, AstChannel *, int * );
-static void FitPoly1DInit( AstPolyMap *, int, double **, AstMinPackData *, double *, int *);
-static void FitPoly2DInit( AstPolyMap *, int, double **, AstMinPackData *, double *, int *);
-static void FreeArrays( AstPolyMap *, int, int * );
-static void IterInverse( AstPolyMap *, AstPointSet *, AstPointSet *, int * );
-static void LMFunc1D(  const double *, double *, int, int, void * );
-static void LMFunc2D(  const double *, double *, int, int, void * );
-static void LMJacob1D( const double *, double *, int, int, void * );
-static void LMJacob2D( const double *, double *, int, int, void * );
-static void PolyCoeffs( AstPolyMap *, int, int, double *, int *, int * );
-static void PolyPowers( AstPolyMap *, double **, int, const int *, double **, int, int, int * );
-static void StoreArrays( AstPolyMap *, int, int, const double *, int * );
-
-#if defined(THREAD_SAFE)
-static int ManageLock( AstObject *, int, int, AstObject **, int * );
-#endif
-
-static const char *GetAttrib( AstObject *, const char *, int * );
-static int TestAttrib( AstObject *, const char *, int * );
-static void ClearAttrib( AstObject *, const char *, int * );
-static void SetAttrib( AstObject *, const char *, int * );
-
-static int GetIterInverse( AstPolyMap *, int * );
-static int TestIterInverse( AstPolyMap *, int * );
-static void ClearIterInverse( AstPolyMap *, int * );
-static void SetIterInverse( AstPolyMap *, int, int * );
-
-static int GetNiterInverse( AstPolyMap *, int * );
-static int TestNiterInverse( AstPolyMap *, int * );
-static void ClearNiterInverse( AstPolyMap *, int * );
-static void SetNiterInverse( AstPolyMap *, int, int * );
-
-static double GetTolInverse( AstPolyMap *, int * );
-static int TestTolInverse( AstPolyMap *, int * );
-static void ClearTolInverse( AstPolyMap *, int * );
-static void SetTolInverse( AstPolyMap *, double, int * );
-
-
-/* Member functions. */
-/* ================= */
-
-static void ClearAttrib( AstObject *this_object, const char *attrib, int *status ) {
-/*
-*  Name:
-*     ClearAttrib
-
-*  Purpose:
-*     Clear an attribute value for a PolyMap.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     void ClearAttrib( AstObject *this, const char *attrib, int *status )
-
-*  Class Membership:
-*     PolyMap member function (over-rides the astClearAttrib protected
-*     method inherited from the Mapping class).
-
-*  Description:
-*     This function clears the value of a specified attribute for a
-*     PolyMap, so that the default value will subsequently be used.
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     attrib
-*        Pointer to a null-terminated string specifying the attribute
-*        name.  This should be in lower case with no surrounding white
-*        space.
-*     status
-*        Pointer to the inherited status variable.
-*/
-
-/* Local Variables: */
-   AstPolyMap *this;             /* Pointer to the PolyMap structure */
-
-/* Check the global error status. */
-   if ( !astOK ) return;
-
-/* Obtain a pointer to the PolyMap structure. */
-   this = (AstPolyMap *) this_object;
-
-/* Check the attribute name and clear the appropriate attribute. */
-
-/* IterInverse. */
-/* ------------ */
-   if ( !strcmp( attrib, "iterinverse" ) ) {
-      astClearIterInverse( this );
-
-/* NiterInverse. */
-/* ------------- */
-   } else if ( !strcmp( attrib, "niterinverse" ) ) {
-      astClearNiterInverse( this );
-
-/* TolInverse. */
-/* ----------- */
-   } else if ( !strcmp( attrib, "tolinverse" ) ) {
-      astClearTolInverse( this );
-
-/* If the attribute is still not recognised, pass it on to the parent
-   method for further interpretation. */
-   } else {
-      (*parent_clearattrib)( this_object, attrib, status );
-   }
-}
-
-static int Equal( AstObject *this_object, AstObject *that_object, int *status ) {
-/*
-*  Name:
-*     Equal
-
-*  Purpose:
-*     Test if two PolyMaps are equivalent.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     int Equal( AstObject *this, AstObject *that, int *status )
-
-*  Class Membership:
-*     PolyMap member function (over-rides the astEqual protected
-*     method inherited from the astMapping class).
-
-*  Description:
-*     This function returns a boolean result (0 or 1) to indicate whether
-*     two PolyMaps are equivalent.
-
-*  Parameters:
-*     this
-*        Pointer to the first Object (a PolyMap).
-*     that
-*        Pointer to the second Object.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     One if the PolyMaps are equivalent, zero otherwise.
-
-*  Notes:
-*     - A value of zero will be returned if this function is invoked
-*     with the global status set, or if it should fail for any reason.
-*/
-
-/* Local Variables: */
-   AstPolyMap *that;
-   AstPolyMap *this;
-   int i, j, k;
-   int nin;
-   int nout;
-   int result;
-   int tmp;
-
-/* Initialise. */
-   result = 0;
-
-/* Check the global error status. */
-   if ( !astOK ) return result;
-
-/* Obtain pointers to the two PolyMap structures. */
-   this = (AstPolyMap *) this_object;
-   that = (AstPolyMap *) that_object;
-
-/* Check the second object is a PolyMap. We know the first is a
-   PolyMap since we have arrived at this implementation of the virtual
-   function. */
-   if( astIsAPolyMap( that ) ) {
-
-/* Get the number of inputs and outputs and check they are the same for both. */
-      nin = astGetNin( this );
-      nout = astGetNout( this );
-      if( astGetNin( that ) == nin && astGetNout( that ) == nout ) {
-
-/* If the Invert flags for the two PolyMaps differ, it may still be possible
-   for them to be equivalent. First compare the PolyMaps if their Invert
-   flags are the same. In this case all the attributes of the two PolyMaps
-   must be identical. */
-         if( astGetInvert( this ) == astGetInvert( that ) ) {
-
-/* Assume they are the same until we find a difference. */
-            result = 1;
-
-/* The "_f" and "_i" suffixes in the PolyMap structure array names refer
-   to the forward and inverse transformations of the original uninverted
-   PolyMap. So we need to ensure the "nout" and "nin" values also refer
-   to the uninverted values. So if the PolyMaps are inverted, swap nout
-   and nin. */
-            if(  astGetInvert( this ) ) {
-               tmp = nin;
-               nin = nout;
-               nout = tmp;
-            }
-
-/* Check properties of the forward transformation. */
-            if( this->ncoeff_f && that->ncoeff_f ) {
-               for( i = 0; i < nout && result; i++ ) {
-                  if( this->ncoeff_f[ i ] != that->ncoeff_f[ i ] ){
-                     result = 0;
-                  }
-               }
-            } else if( this->ncoeff_f || that->ncoeff_f ) {
-               result = 0;
-            }
-
-            if( this->mxpow_f && that->mxpow_f ) {
-               for( i = 0; i < nout && result; i++ ) {
-                  if( this->mxpow_f[ i ] != that->mxpow_f[ i ] ){
-                     result = 0;
-                  }
-               }
-            } else if( this->mxpow_f || that->mxpow_f ) {
-               result = 0;
-            }
-
-            if( this->coeff_f && that->coeff_f ) {
-               for( i = 0; i < nout && result; i++ ) {
-                  for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) {
-                     if( !astEQUAL( this->coeff_f[ i ][ j ],
-                                    that->coeff_f[ i ][ j ] ) ) {
-                        result = 0;
-                     }
-                  }
-               }
-            } else if( this->coeff_f || that->coeff_f ) {
-               result = 0;
-            }
-
-            if( this->power_f && that->power_f ) {
-               for( i = 0; i < nout && result; i++ ) {
-                  for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) {
-                     for( k = 0; k < nin && result; k++ ) {
-                        if( this->power_f[ i ][ j ][ k ] !=
-                            that->power_f[ i ][ j ][ k ] ) {
-                           result = 0;
-                        }
-                     }
-                  }
-               }
-            } else if( this->power_f || that->power_f ) {
-               result = 0;
-            }
-
-/* Check properties of the inverse transformation. */
-            if( this->ncoeff_i && that->ncoeff_i ) {
-               for( i = 0; i < nout && result; i++ ) {
-                  if( this->ncoeff_i[ i ] != that->ncoeff_i[ i ] ){
-                     result = 0;
-                  }
-               }
-            } else if( this->ncoeff_i || that->ncoeff_i ) {
-               result = 0;
-            }
-
-            if( this->mxpow_i && that->mxpow_i ) {
-               for( i = 0; i < nout && result; i++ ) {
-                  if( this->mxpow_i[ i ] != that->mxpow_i[ i ] ){
-                     result = 0;
-                  }
-               }
-            } else if( this->mxpow_i || that->mxpow_i ) {
-               result = 0;
-            }
-
-            if( this->coeff_f && that->coeff_f ) {
-               for( i = 0; i < nout && result; i++ ) {
-                  for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) {
-                     if( !astEQUAL( this->coeff_f[ i ][ j ],
-                                    that->coeff_f[ i ][ j ] ) ) {
-                        result = 0;
-                     }
-                  }
-               }
-            } else if( this->coeff_f || that->coeff_f ) {
-               result = 0;
-            }
-
-            if( this->power_f && that->power_f ) {
-               for( i = 0; i < nout && result; i++ ) {
-                  for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) {
-                     for( k = 0; k < nin && result; k++ ) {
-                        if( this->power_f[ i ][ j ][ k ] !=
-                            that->power_f[ i ][ j ][ k ] ) {
-                           result = 0;
-                        }
-                     }
-                  }
-               }
-            } else if( this->power_f || that->power_f ) {
-               result = 0;
-            }
-
-/* If the Invert flags for the two PolyMaps differ, the attributes of the two
-   PolyMaps must be inversely related to each other. */
-         } else {
-
-/* In the specific case of a PolyMap, Invert flags must be equal. */
-            result = 0;
-
-         }
-      }
-   }
-
-/* If an error occurred, clear the result value. */
-   if ( !astOK ) result = 0;
-
-/* Return the result, */
-   return result;
-}
-
-static double *FitPoly1D( AstPolyMap *this, int forward, int nsamp, double acc,
-                          int order, double **table, double scales[2], int *ncoeff,
-                          double *racc, int *status ){
-/*
-*  Name:
-*     FitPoly1D
-
-*  Purpose:
-*     Fit a (1-in,1-out) polynomial to a supplied set of data.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     double *FitPoly1D( AstPolyMap *this, int forward, int nsamp, double acc,
-*                        int order, double **table, double scales[2], int *ncoeff,
-*                        double *racc, int *status )
-
-*  Description:
-*     This function fits a least squares 1D polynomial curve to the
-*     positions in a supplied table, and returns the coefficients of a
-*     PolyMap to describe the fit. For the purposes of this function,
-*     the input to the fit is refered to as x1 and the output as y1. So
-*     the returned coefficients describe a PolyMap with forward
-*     transformation:
-*
-*     y1 = P1( x1 )
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     forward
-*        Non-zero if the forward transformation of "this" is being
-*        replaced. Zero if the inverse transformation of "this" is being
-*        replaced.
-*     nsamp
-*        The number of (x1,y1) positions in the supplied table.
-*     acc
-*        The required accuracy, expressed as a geodesic distance within
-*        the polynomials output space (not the normalised tabular space).
-*     order
-*        The maximum power (plus one) of x1 within P1. So for instance, a
-*        value of 3 refers to a quadratic polynomial.
-*     table
-*        Pointer to an array of 2 pointers. Each of these pointers points
-*        to an array of "nsamp" doubles, being the normalised and sampled
-*        values for x1 and y1 in that order.
-*     scales
-*        Array holding the scaling factors used to produced the normalised
-*        values in the two columns of the table. Multiplying the normalised
-*        table values by the scale factor produces input or output axis
-*        values for the returned PolyMap
-*     ncoeff
-*        Pointer to an int in which to return the number of coefficients
-*        described by the returned array.
-*     racc
-*        Pointer to a double in which to return the achieved accuracy
-*        (which may be greater than "acc").
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     A pointer to an array of doubles defining the polynomial in the
-*     form required by the PolyMap contructor. The number of coefficients
-*     is returned via "ncoeff". If the polynomial could not be found,
-*     then NULL is returned. The returned pointer should be freed using
-*     astFree when no longer needed.
-
-*/
-
-/* Local Variables: */
-   AstMinPackData data;
-   double *coeffs;
-   double *pc;
-   double *pr;
-   double *pxp1;
-   double *result;
-   double *work1;
-   double *work2;
-   double *work4;
-   double f1;
-   double f2;
-   double maxterm;
-   double term;
-   double tv;
-   int *work3;
-   int info;
-   int k;
-   int ncof;
-   int w1;
-
-/* Initialise returned value */
-   result = NULL;
-   *ncoeff = 0;
-   *racc = 10*acc;
-
-/* Check inherited status */
-   if( !astOK ) return result;
-
-/* Number of coefficients per poly. */
-   ncof = order;
-
-/* Initialise the elements of the structure. */
-   data.order = order;
-   data.nsamp = nsamp;
-   data.init_jac = 1;
-   data.xp1 = astMalloc( nsamp*order*sizeof( double ) );
-   data.xp2 = NULL;
-   data.y[ 0 ] = table[ 1 ];
-   data.y[ 1 ] = NULL;
-
-/* Work space to hold coefficients. */
-   coeffs = astMalloc( ncof*sizeof( double ) );
-
-/* Other work space. */
-   work1 = astMalloc( nsamp*sizeof( double ) );
-   work2 = astMalloc( ncof*nsamp*sizeof( double ) );
-   work3 = astMalloc( ncof*sizeof( int ) );
-   work4 = astMalloc( (5*ncof+nsamp)*sizeof( double ) );
-   if( astOK ) {
-
-/* Find all the required powers of x1 and store them in the "xp1"
-   component of the data structure. The required initialisation is done
-   differently for different subclasses of PolyMap, so we need to wrap
-   it up in a virtual function. */
-      astFitPoly1DInit( this, forward, table, &data, scales );
-
-/* The initial guess at the coefficient values represents a unit
-   transformation in (normalised) tabulated (x,y) values. Using normalised
-   values means that we are, effectively, including a guess at the linear
-   scaling factor between input and output of the PolyMap (e.g. the PolyMap
-   may have inputs in mm and outputs in radians). */
-      for( k = 0; k < ncof; k++ ) coeffs[ k ] = 0.0;
-      coeffs[ 1 ] = 1.0;
-
-/* Find the best coefficients */
-      info = lmder1( MPFunc1D, &data, nsamp, ncof, coeffs, work1, work2, nsamp,
-                     sqrt(DBL_EPSILON), work3, work4, (5*ncof+nsamp) );
-      if( info == 0 ) astError( AST__MNPCK, "astPolyMap(PolyTran): Minpack error "
-                                "detected (possible programming error).", status );
-
-/* Return the achieved accuracy. The "work1" array holds the normalised Y
-   residuals at each tabulated point. */
-      pr = work1;
-      tv = 0.0;
-      for( k = 0; k < nsamp; k++,pr++ ) tv += (*pr)*(*pr);
-      *racc = scales[ 1 ]*sqrt( tv/nsamp );
-
-/* The best fitting polynomial coefficient found above relate to the
-   polynomial between the scaled positions stored in "table". These
-   scaled positions are related to PolyMap input/output axis values via
-   the scale factors supplied in "scales". Find the initial factor for the
-   current output. */
-      f1 = scales[ 1 ];
-      f2 = 1.0;
-
-/* Look at each coefficient. */
-      pc = coeffs;
-      for( w1 = 0; w1 < order; w1++,pc++ ) {
-
-/* Get a pointer to the powers of X1 appropriate for the current coefficient,
-    at the first sample. */
-         pxp1 = data.xp1 + w1;
-
-/* We find the contribution which this coefficient makes to the total
-   polynomial value. Find the maximum contribution made at any sample
-   points. */
-         maxterm = 0.0;
-         for( k = 0; k < nsamp; k++ ) {
-
-/* Get the absolute value of the polynomial term that uses the current
-   coefficient. */
-            term = fabs( ( *pc )*( *pxp1 ) );
-
-/* Update the maximum term found at any sample. */
-            if( term > maxterm ) maxterm = term;
-
-/* Increment the pointers to refer to the next sample. */
-            pxp1 += order;
-         }
-
-/* If the maximum contribution made by this term is less than the
-   required accuracy, set the coefficient value to zero. */
-         if( maxterm*f1 < acc ) {
-            *pc = 0.0;
-
-/* Scale the best fitting polynomial coefficient found above to take
-   account of the fact that the tabulated input and output positions in
-   "table" were are not actual PolyMap input and output axis values, but
-   are scaled by the factors stored in "scales". */
-         } else {
-            *pc *= f1/f2;
-         }
-
-         f2 *= scales[ 0 ];
-      }
-
-/* Convert the array of coefficients into PolyMap form. */
-      result = astMalloc( ncof*3*sizeof( double ) );
-
-      *ncoeff = 0;
-      pr = result;
-      pc = coeffs;
-      for( w1 = 0; w1 < order; w1++,pc++ ) {
-         if( *pc != 0.0 ) {
-            *(pr++) = *pc;
-            *(pr++) = 1;
-            *(pr++) = w1;
-            (*ncoeff)++;
-         }
-      }
-
-/* Truncate the returned array. */
-      result = astRealloc( result, (*ncoeff)*3*sizeof( double ) );
-   }
-
-/* Free resources. */
-   coeffs = astFree( coeffs );
-   data.xp1 = astFree( data.xp1 );
-   work1 = astFree( work1 );
-   work2 = astFree( work2 );
-   work3 = astFree( work3 );
-   work4 = astFree( work4 );
-
-/* Return the coefficient array. */
-   return result;
-
-}
-
-static void FitPoly1DInit( AstPolyMap *this, int forward, double **table,
-                           AstMinPackData *data, double *scales, int *status ){
-/*
-*+
-*  Name:
-*     astFitPoly1DInit
-
-*  Purpose:
-*     Perform initialisation needed for FitPoly1D
-
-*  Type:
-*     Protected function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     void astFitPoly1DInit( AstPolyMap *this, int forward, double **table,
-*                            AstMinPackData *data, double *scales,
-*                            int *status )
-
-*  Class Membership:
-*     PolyMap virtual function.
-
-*  Description:
-*     This function performs initialisation needed for FitPoly1D.
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     forward
-*        Non-zero if the forward transformation of "this" is being
-*        replaced. Zero if the inverse transformation of "this" is being
-*        replaced.
-*     table
-*        Pointer to an array of 2 pointers. Each of these pointers points
-*        to an array of "nsamp" doubles, being the scaled and sampled values
-*        for x1 and y1 in that order.
-*     data
-*        Pointer to a structure holding information to pass the the
-*        service function invoked by the minimisation function.
-*     scales
-*        Array holding the scaling factors for the two columns of the table.
-*        Multiplying the table values by the scale factor produces PolyMap
-*        input or output axis values.
-*-
-*/
-
-/* Local Variables; */
-   double *px1;
-   double *pxp1;
-   double tv;
-   double x1;
-   int k;
-   int w1;
-
-/* Check the local error status. */
-   if ( !astOK ) return;
-
-/* Get pointers to the supplied x1 values. */
-   px1 = table[ 0 ];
-
-/* Get pointers to the location for the next power of x1. */
-   pxp1 = data->xp1;
-
-/* Loop round all samples. */
-   for( k = 0; k < data->nsamp; k++ ) {
-
-/* Get the current x1 value. */
-      x1 = *(px1++);
-
-/* Find all the required powers of x1 and store them in the "xp1"
-   component of the data structure. */
-      tv = 1.0;
-      for( w1 = 0; w1 < data->order; w1++ ) {
-         *(pxp1++) = tv;
-         tv *= x1;
-      }
-   }
-}
-
-static double *FitPoly2D( AstPolyMap *this, int forward, int nsamp, double acc,
-                          int order, double **table, double scales[4],
-                          int *ncoeff, double *racc, int *status ){
-/*
-*  Name:
-*     FitPoly2D
-
-*  Purpose:
-*     Fit a (2-in,2-out) polynomial to a supplied set of data.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     double *FitPoly2D( AstPolyMap *this, int forward, int nsamp, double acc,
-*                        int order, double **table, double scales[4],
-*                        int *ncoeff, double *racc, int *status )
-
-*  Description:
-*     This function fits a pair of least squares 2D polynomial surfaces
-*     to the positions in a supplied table, and returns the coefficients
-*     of a PolyMap to describe the fit. For the purposes of this function,
-*     the inputs to the fit is refered to as (x1,x2) and the output as
-*     (y1,y2). So the returned coefficients describe a PolyMap with forward
-*     transformations:
-*
-*     y1 = P1( x1, x2 )
-*     y2 = P2( x1, x2 )
-*
-*     P1 and P2 have the same maximum power on each input (specified by
-*     the "order" parameter).
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     forward
-*        Non-zero if the forward transformation of "this" is being
-*        replaced. Zero if the inverse transformation of "this" is being
-*        replaced.
-*     nsamp
-*        The number of (x1,x2,y1,y2) positions in the supplied table.
-*     acc
-*        The required accuracy, expressed as a geodesic distance within
-*        the polynomials output space (not the normalised tabular space).
-*     order
-*        The maximum power (plus one) of x1 or x2 within P1 and P2. So for
-*        instance, a value of 3 refers to a quadratic polynomial. Note, cross
-*        terms with total powers greater than or equal to "order" are not
-*        included in the fit. So the maximum number of terms in the fitted
-*        polynomial is order*(order+1)/2.
-*     table
-*        Pointer to an array of 4 pointers. Each of these pointers points
-*        to an array of "nsamp" doubles, being the normalised and sampled
-*        values for x1, x2, y1 or y2 in that order.
-*     scales
-*        Array holding the scaling factors used to produced the normalised
-*        values in the four columns of the table. Multiplying the normalised
-*        table values by the scale factor produces input or output axis
-*        values for the returned PolyMap
-*     ncoeff
-*        Pointer to an ant in which to return the number of coefficients
-*        described by the returned array.
-*     racc
-*        Pointer to a double in which to return the achieved accuracy
-*        (which may be greater than "acc").
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     A pointer to an array of doubles defining the polynomial in the
-*     form required by the PolyMap contructor. The number of coefficients
-*     is returned via "ncoeff". If the polynomial could not be found with
-*     sufficient accuracy , then NULL is returned. The returned pointer
-*     should be freed using astFree when no longer needed.
-
-*/
-
-/* Local Variables: */
-   AstMinPackData data;
-   double *coeffs;
-   double *pc;
-   double *pr;
-   double *pxp1;
-   double *pxp2;
-   double *result;
-   double *work1;
-   double *work2;
-   double *work4;
-   double f1;
-   double f20;
-   double f2;
-   double f3;
-   double facc;
-   double maxterm;
-   double term;
-   double tv;
-   int *work3;
-   int info;
-   int iout;
-   int k;
-   int ncof;
-   int w12;
-   int w1;
-   int w2;
-
-/* Initialise returned value */
-   result = NULL;
-   *ncoeff = 0;
-   *racc = 10*acc;
-
-/* Check inherited status */
-   if( !astOK ) return result;
-
-/* Number of coefficients per poly. */
-   ncof = order*( order + 1 )/2;
-
-/* Initialise the elements of the structure. */
-   data.order = order;
-   data.nsamp = nsamp;
-   data.init_jac = 1;
-   data.xp1 = astMalloc( nsamp*order*sizeof( double ) );
-   data.xp2 = astMalloc( nsamp*order*sizeof( double ) );
-   data.y[ 0 ] = table[ 2 ];
-   data.y[ 1 ] = table[ 3 ];
-
-/* Work space to hold coefficients. */
-   coeffs = astMalloc( 2*ncof*sizeof( double ) );
-
-/* Other work space. */
-   work1 = astMalloc( 2*nsamp*sizeof( double ) );
-   work2 = astMalloc( 4*ncof*nsamp*sizeof( double ) );
-   work3 = astMalloc( 2*ncof*sizeof( int ) );
-   work4 = astMalloc( 2*(5*ncof+nsamp)*sizeof( double ) );
-   if( astOK ) {
-
-/* Find all the required powers of x1 and x2 and store them in the "xp1"
-   and "xp2" components of the data structure. The required initialisation
-   is done differently for different subclasses of PolyMap, so we need to
-   wrap it up in a virtual function. */
-      astFitPoly2DInit( this, forward, table, &data, scales );
-
-/* The initial guess at the coefficient values represents a unit
-   transformation in (normalised) tabulated (x,y) values. Using normalised
-   values means that we are, effectively, including a guess at the linear
-   scaling factor between input and output of the PolyMap (e.g. the PolyMap
-   may have inputs in mm and outputs in radians). */
-      for( k = 0; k < 2*ncof; k++ ) coeffs[ k ] = 0.0;
-      coeffs[ 1 ] = 1.0;
-      coeffs[ 5 ] = 1.0;
-
-/* Find the best coefficients */
-      info = lmder1( MPFunc2D, &data, 2*nsamp, 2*ncof, coeffs, work1, work2,
-                     2*nsamp, sqrt(DBL_EPSILON), work3, work4, 2*(5*ncof+nsamp) );
-      if( info == 0 ) astError( AST__MNPCK, "astPolyMap(PolyTran): Minpack error "
-                                "detected (possible programming error).", status );
-
-/* Return the achieved accuracy. */
-      pr = work1;
-      tv = 0.0;
-      for( k = 0; k < 2*nsamp; k++,pr++ ) tv += (*pr)*(*pr);
-      facc = 1.0/(scales[2]*scales[2]) + 1.0/(scales[3]*scales[3]);
-      *racc = sqrt( tv/(2*nsamp*facc) );
-
-/* Pointer to the first coefficient. */
-      pc = coeffs;
-
-/* Look at coefficients for each output in turn. */
-      for( iout = 0; iout < 2 && astOK; iout++ ) {
-
-/* The best fitting polynomial coefficient found above relate to the
-   polynomial between the scaled positions stored in "table". These
-   scaled positions are related to PolyMap input/output axis values via
-   the scale factors supplied in "scales". Find the initial factor for the
-   current output. */
-         f1 = scales[ 2 + iout ];
-
-/* Look at each coefficient for the current output. */
-         f20 = 1.0;
-         for( w12 = 0; w12 < order; w12++ ) {
-            f3 = 1.0;
-            f2 = f20;
-            for( w2 = 0; w2 <= w12; w2++,pc++ ) {
-               w1 = w12 - w2;
-
-/* Get pointers to the powers of X1 and X2 appropriate for the current
-   coefficient, at the first sample. */
-               pxp1 = data.xp1 + w1;
-               pxp2 = data.xp2 + w2;
-
-/* We find the contribution which this coefficient makes to the total
-   polynomial value. Find the maximum contribution made at any sample
-   points. */
-               maxterm = 0.0;
-               for( k = 0; k < nsamp; k++ ) {
-
-/* Get the absolute value of the polynomial term that uses the current
-   coefficient. */
-                  term = fabs( ( *pc )*( *pxp1 )*( *pxp2 ) );
-
-/* Update the maximum term found at any sample. */
-                  if( term > maxterm ) maxterm = term;
-
-/* Increment the pointers to refer to the next sample. */
-                  pxp1 += order;
-                  pxp2 += order;
-               }
-
-/* If the maximum contribution made by this term is less than the
-   required accuracy, set the coefficient value to zero. */
-               if( maxterm*f1 < acc ) {
-                  *pc = 0.0;
-
-/* Scale the best fitting polynomial coefficient found above to take
-   account of the fact that the tabulated input and output positions in
-   "table" were are not actual PolyMap input and output axis values, but
-   are scaled by the factors stored in "scales". */
-               } else {
-                  *pc *= f1/( f2*f3 );
-               }
-
-               f2 /= scales[ 0 ];
-               f3 *= scales[ 1 ];
-            }
-
-            f20 *= scales[ 0 ];
-         }
-      }
-
-/* Convert the array of coefficients into PolyMap form. */
-      result = astMalloc( 2*ncof*4*sizeof( double ) );
-
-      *ncoeff = 0;
-      pr = result;
-      pc = coeffs;
-      for( iout = 0; iout < 2 && astOK; iout++ ) {
-         for( w12 = 0; w12 < order; w12++ ) {
-            for( w2 = 0; w2 <= w12; w2++,pc++ ) {
-               w1 = w12 - w2;
-               if( *pc != 0.0 ) {
-                  *(pr++) = *pc;
-                  *(pr++) = iout + 1;
-                  *(pr++) = w1;
-                  *(pr++) = w2;
-                  (*ncoeff)++;
-               }
-            }
-         }
-      }
-
-/* Truncate the returned array. */
-      result = astRealloc( result, (*ncoeff)*4*sizeof( double ) );
-   }
-
-/* Free resources. */
-   coeffs = astFree( coeffs );
-   data.xp1 = astFree( data.xp1 );
-   data.xp2 = astFree( data.xp2 );
-   work1 = astFree( work1 );
-   work2 = astFree( work2 );
-   work3 = astFree( work3 );
-   work4 = astFree( work4 );
-
-/* Return the coefficient array. */
-   return result;
-
-}
-
-static void FitPoly2DInit( AstPolyMap *this, int forward, double **table,
-                           AstMinPackData *data, double *scales, int *status ){
-/*
-*+
-*  Name:
-*     astFitPoly2DInit
-
-*  Purpose:
-*     Perform initialisation needed for FitPoly2D
-
-*  Type:
-*     Protected function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     void astFitPoly2DInit( AstPolyMap *this, int forward, double **table,
-*                            AstMinPackData *data, double *scales,
-*                            int *status )
-
-*  Class Membership:
-*     PolyMap virtual function.
-
-*  Description:
-*     This function performs initialisation needed for FitPoly2D.
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     forward
-*        Non-zero if the forward transformation of "this" is being
-*        replaced. Zero if the inverse transformation of "this" is being
-*        replaced.
-*     table
-*        Pointer to an array of 4 pointers. Each of these pointers points
-*        to an array of "nsamp" doubles, being the scaled and sampled values
-*        for x1, x2, y1 or y2 in that order.
-*     data
-*        Pointer to a structure holding information to pass the the
-*        service function invoked by the minimisation function.
-*     scales
-*        Array holding the scaling factors for the four columns of the table.
-*        Multiplying the table values by the scale factor produces PolyMap
-*        input or output axis values.
-*-
-*/
-
-/* Local Variables; */
-   double *px1;
-   double *px2;
-   double *pxp1;
-   double *pxp2;
-   double tv;
-   double x1;
-   double x2;
-   int k;
-   int w1;
-   int w2;
-
-/* Check the local error status. */
-   if ( !astOK ) return;
-
-/* Get pointers to the supplied x1 and x2 values. */
-   px1 = table[ 0 ];
-   px2 = table[ 1 ];
-
-/* Get pointers to the location for the next power of x1 and x2. */
-   pxp1 = data->xp1;
-   pxp2 = data->xp2;
-
-/* Loop round all samples. */
-   for( k = 0; k < data->nsamp; k++ ) {
-
-/* Get the current x1 and x2 values. */
-      x1 = *(px1++);
-      x2 = *(px2++);
-
-/* Find all the required powers of x1 and store them in the "xp1"
-   component of the data structure. */
-      tv = 1.0;
-      for( w1 = 0; w1 < data->order; w1++ ) {
-         *(pxp1++) = tv;
-         tv *= x1;
-      }
-
-/* Find all the required powers of x2 and store them in the "xp2"
-   comonent of the data structure. */
-      tv = 1.0;
-      for( w2 = 0; w2 < data->order; w2++ ) {
-         *(pxp2++) = tv;
-         tv *= x2;
-      }
-   }
-}
-
-static void FreeArrays( AstPolyMap *this, int forward, int *status ) {
-/*
-*  Name:
-*     FreeArrays
-
-*  Purpose:
-*     Free the dynamic arrays contained within a PolyMap
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     void FreeArrays( AstPolyMap *this, int forward, int *status )
-
-*  Description:
-*     This function frees all the dynamic arrays allocated as part of a
-*     PolyMap.
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     forward
-*        If non-zero, the arrays for the forward transformation are freed.
-*        Otherwise, the arrays for the inverse transformation are freed.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     void
-
-*  Notes:
-*     This function attempts to execute even if the global error status is
-*     set.
-*/
-
-/* Local Variables: */
-   int nin;                     /* Number of inputs */
-   int nout;                    /* Number of outputs */
-   int i;                       /* Loop count */
-   int j;                       /* Loop count */
-
-/* Get the number of inputs and outputs of the uninverted Mapping. */
-   nin = ( (AstMapping *) this )->nin;
-   nout = ( (AstMapping *) this )->nout;
-
-/* Free the dynamic arrays for the forward transformation. */
-   if( forward != astGetInvert( this ) ) {
-
-      if( this->coeff_f ) {
-         for( i = 0; i < nout; i++ ) {
-            this->coeff_f[ i ] = astFree( this->coeff_f[ i ] );
-         }
-         this->coeff_f = astFree( this->coeff_f );
-      }
-
-      if( this->power_f ) {
-         for( i = 0; i < nout; i++ ) {
-            if( this->ncoeff_f && this->power_f[ i ] ) {
-               for( j = 0; j < this->ncoeff_f[ i ]; j++ ) {
-                  this->power_f[ i ][ j ] = astFree( this->power_f[ i ][ j ] );
-               }
-            }
-            this->power_f[ i ] = astFree( this->power_f[ i ] );
-         }
-         this->power_f = astFree( this->power_f );
-      }
-
-      this->ncoeff_f = astFree( this->ncoeff_f );
-      this->mxpow_f = astFree( this->mxpow_f );
-
-/* Free the dynamic arrays for the inverse transformation. */
-   } else {
-
-      if( this->coeff_i ) {
-         for( i = 0; i < nin; i++ ) {
-            this->coeff_i[ i ] = astFree( this->coeff_i[ i ] );
-         }
-         this->coeff_i = astFree( this->coeff_i );
-      }
-
-      if(this->power_i ) {
-         for( i = 0; i < nin; i++ ) {
-            if( this->ncoeff_i && this->power_i[ i ] ) {
-               for( j = 0; j < this->ncoeff_i[ i ]; j++ ) {
-                  this->power_i[ i ][ j ] = astFree( this->power_i[ i ][ j ] );
-               }
-            }
-            this->power_i[ i ] = astFree( this->power_i[ i ] );
-         }
-         this->power_i = astFree( this->power_i );
-      }
-
-      this->ncoeff_i = astFree( this->ncoeff_i );
-      this->mxpow_i = astFree( this->mxpow_i );
-   }
-}
-
-static const char *GetAttrib( AstObject *this_object, const char *attrib, int *status ) {
-/*
-*  Name:
-*     GetAttrib
-
-*  Purpose:
-*     Get the value of a specified attribute for a PolyMap.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     const char *GetAttrib( AstObject *this, const char *attrib, int *status )
-
-*  Class Membership:
-*     PolyMap member function (over-rides the protected astGetAttrib
-*     method inherited from the Mapping class).
-
-*  Description:
-*     This function returns a pointer to the value of a specified
-*     attribute for a PolyMap, formatted as a character string.
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     attrib
-*        Pointer to a null-terminated string containing the name of
-*        the attribute whose value is required. This name should be in
-*        lower case, with all white space removed.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     - Pointer to a null-terminated string containing the attribute
-*     value.
-
-*  Notes:
-*     - The returned string pointer may point at memory allocated
-*     within the PolyMap, or at static memory. The contents of the
-*     string may be over-written or the pointer may become invalid
-*     following a further invocation of the same function or any
-*     modification of the PolyMap. A copy of the string should
-*     therefore be made if necessary.
-*     - A NULL pointer will be returned if this function is invoked
-*     with the global error status set, or if it should fail for any
-*     reason.
-*/
-
-/* Local Variables: */
-   astDECLARE_GLOBALS           /* Pointer to thread-specific global data */
-   AstPolyMap *this;            /* Pointer to the PolyMap structure */
-   const char *result;          /* Pointer value to return */
-   double dval;                 /* Floating point attribute value */
-   int ival;                    /* Integer attribute value */
-
-/* Initialise. */
-   result = NULL;
-
-/* Check the global error status. */
-   if ( !astOK ) return result;
-
-/* Get a pointer to the thread specific global data structure. */
-   astGET_GLOBALS(this_object);
-
-/* Obtain a pointer to the PolyMap structure. */
-   this = (AstPolyMap *) this_object;
-
-/* Compare "attrib" with each recognised attribute name in turn,
-   obtaining the value of the required attribute. If necessary, write
-   the value into "getattrib_buff" as a null-terminated string in an appropriate
-   format.  Set "result" to point at the result string. */
-
-/* IterInverse. */
-/* ------------ */
-   if ( !strcmp( attrib, "iterinverse" ) ) {
-      ival = astGetIterInverse( this );
-      if ( astOK ) {
-         (void) sprintf( getattrib_buff, "%d", ival );
-         result = getattrib_buff;
-      }
-
-/* NiterInverse. */
-/* ------------- */
-   } else if ( !strcmp( attrib, "niterinverse" ) ) {
-      ival = astGetNiterInverse( this );
-      if ( astOK ) {
-         (void) sprintf( getattrib_buff, "%d", ival );
-         result = getattrib_buff;
-      }
-
-/* TolInverse. */
-/* ----------- */
-   } else if ( !strcmp( attrib, "tolinverse" ) ) {
-      dval = astGetTolInverse( this );
-      if ( astOK ) {
-         (void) sprintf( getattrib_buff, "%.*g", AST__DBL_DIG, dval );
-         result = getattrib_buff;
-      }
-
-/* If the attribute name was not recognised, pass it on to the parent
-   method for further interpretation. */
-   } else {
-      result = (*parent_getattrib)( this_object, attrib, status );
-   }
-
-/* Return the result. */
-   return result;
-
-}
-
-static AstPolyMap **GetJacobian( AstPolyMap *this, int *status ){
-/*
-*  Name:
-*     GetJacobian
-
-*  Purpose:
-*     Get a description of a Jacobian matrix for the original
-*     fwd transformation of a PolyMap.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     AstPolyMap **GetJacobian( AstPolyMap *this, int *status )
-
-*  Description:
-*     This function returns a set of PolyMaps which define the Jacobian
-*     matrix of the original forward transformation of the supplied PolyMap
-*     (i.e. the Negated attribute is assumed to be zero).
-*
-*     The Jacobian matrix has "nout" rows and "nin" columns, where "nin"
-*     and "nout" are the number of inputs and outputs of the original PolyMap.
-*     Row "i", column "j" of the matrix holds the rate of change of the
-*     i'th PolyMap output with respect to the j'th PolyMap input.
-*
-*     Since the values in the Jacobian matrix vary across the input space
-*     of the PolyMap, the matrix is returned in the form of a set of new
-*     PolyMaps which generate the elements of the Jacobian for any given
-*     position in the input space. The "nout" values in a single column of
-*     the Jacobian matrix are generated by the "nout" outputs from a single
-*     new PolyMap. The whole matrix is described by "nin" PolyMaps.
-*
-*     The returned PolyMaps are cached in the supplied PolyMap object in
-*     order to speed up subsequent calls to this function.
-
-*  Parameters:
-*     this
-*        The PolyMap for which the Jacbian is required.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     A pointer to an array of "nin" PolyMap pointers, where "nin" is the number
-*     of inputs for the sipplied PolyMap. The returned array should not be changed
-*     in any way, and the PolyMaps should not be freed (they will be freed when
-*     the supplied PolyMap is deleted).
-
-*/
-
-/* Local Variables: */
-   double *coeffs;
-   double *pc;
-   int icof;
-   int icol;
-   int iin;
-   int irow;
-   int ncof;
-   int ncof_row;
-   int ncof_total;
-   int nin;
-   int nout;
-   int power;
-
-/* Check inherited status */
-   if( !astOK ) return NULL;
-
-/* Ensure there is a Jacobian stored in the PolyMap. */
-   if( !this->jacobian ) {
-
-/* Get the number of inputs and outputs of the original forward
-   transformation. */
-      if( astGetInvert( this ) ){
-         nout = astGetNin( this );
-         nin = astGetNout( this );
-      } else {
-         nin = astGetNin( this );
-         nout = astGetNout( this );
-      }
-
-/* Allocate memory to hold pointers to the PolyMaps used to describe the
-   Jacobian matrix. */
-      this->jacobian = astCalloc( nin, sizeof(AstPolyMap *) );
-
-/* Find the total number of coefficients used to describe the supplied
-   PolyMap (forward transformation) and allocate work space to hold the
-   coefficients for a single new PolyMap forward transformation. */
-      ncof = 0;
-      for( irow = 0; irow <  nout; irow++ ) {
-         ncof += this->ncoeff_f[ irow ];
-      }
-      coeffs = astMalloc( ncof*( 2 + nin )*sizeof( double ) );
-
-/* Check pointers can be used safely. */
-      if( astOK ) {
-
-/* The Jacobian matrix has "nout" rows and "nin" columns. The "nout" values
-   in a single column of the Jacobian matrix corresponds to the "nout" outputs
-   from a single PolyMap. The whole matrix is described by "nin" PolyMaps.
-   Loop over each column of the Matrix, creating the corresponding PolyMap
-   for each. */
-         for( icol = 0; icol <  nin; icol++ ) {
-
-/* Initialise the total number of coefficients used to describe the
-   element of the PolyMap. */
-            ncof_total = 0;
-
-/* Loop over each row of the Jacobian matrix (i.e. each PolyMap output). */
-            pc = coeffs;
-            for( irow = 0; irow <  nout; irow++ ) {
-
-/* Loop over each coefficient used in the polynomial that generates the
-   current PolyMap output. */
-               ncof_row = this->ncoeff_f[ irow ];
-               for( icof = 0; icof <  ncof_row; icof++ ) {
-
-/* Get the power of input "icol" associated with the current coefficient. */
-                  power = (int)( this->power_f[ irow ][ icof ][ icol ] + 0.5 );
-
-/* We can skip the coefficient if the power is zero. */
-                  if( power > 0 ) {
-                     ncof_total++;
-
-/* Store the coefficient value, modified so that it describes a
-   polynomial that has been differentiated with respect to input "icol". */
-                     *(pc++) = this->coeff_f[ irow ][ icof ]*power;
-
-/* Store the output PolyMap to which the coeff relates. */
-                     *(pc++) = irow + 1;
-
-/* Store the powers of the inputs associated with the coeff. These are
-   the same as the original powers, except that the power of "icol"
-   (the input with respect to which the output has been differentiated)
-   is reduced by one. */
-                     for( iin = 0; iin <  nin; iin++ ) {
-                        if( iin != icol ) {
-                           *(pc++) = this->power_f[ irow ][ icof ][ iin ];
-                        } else {
-                           *(pc++) = this->power_f[ irow ][ icof ][ iin ] - 1;
-                        }
-                     }
-                  }
-               }
-            }
-
-/* Create the PolyMap and store a pointer to it in the jacobian array in
-   the supplied PolyMap. */
-            (this->jacobian)[ icol ] = astPolyMap( nin, nout, ncof_total, coeffs,
-                                                   0, NULL, " ", status );
-         }
-      }
-
-/* Free resources */
-      coeffs = astFree( coeffs );
-   }
-
-/* Return the Jacobian. */
-   return this->jacobian;
-}
-
-static int GetObjSize( AstObject *this_object, int *status ) {
-/*
-*  Name:
-*     GetObjSize
-
-*  Purpose:
-*     Return the in-memory size of an Object.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     int GetObjSize( AstObject *this, int *status )
-
-*  Class Membership:
-*     PolyMap member function (over-rides the astGetObjSize protected
-*     method inherited from the parent class).
-
-*  Description:
-*     This function returns the in-memory size of the supplied PolyMap,
-*     in bytes.
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     The Object size, in bytes.
-
-*  Notes:
-*     - A value of zero will be returned if this function is invoked
-*     with the global status set, or if it should fail for any reason.
-*/
-
-/* Local Variables: */
-   AstPolyMap *this;
-   int ic;
-   int nc;
-   int result;
-
-/* Initialise. */
-   result = 0;
-
-/* Check the global error status. */
-   if ( !astOK ) return result;
-
-/* Obtain a pointers to the PolyMap structure. */
-   this = (AstPolyMap *) this_object;
-
-/* Invoke the GetObjSize method inherited from the parent class, and then
-   add on any components of the class structure defined by this class
-   which are stored in dynamically allocated memory. */
-   result = (*parent_getobjsize)( this_object, status );
-
-   if( this->jacobian ) {
-      nc = astGetNin( this );
-      for( ic = 0; ic < nc; ic++ ) {
-         result +=  astGetObjSize( (this->jacobian)[ ic ] );
-      }
-      result += sizeof( AstPolyMap * )*nc;
-   }
-
-/* If an error occurred, clear the result value. */
-   if ( !astOK ) result = 0;
-
-/* Return the result, */
-   return result;
-}
-
-static int GetTranForward( AstMapping *this, int *status ) {
-/*
-*
-*  Name:
-*     GetTranForward
-
-*  Purpose:
-*     Determine if a PolyMap defines a forward coordinate transformation.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     #include "mapping.h"
-*     int GetTranForward( AstMapping *this, int *status )
-
-*  Class Membership:
-*     PolyMap member function (over-rides the astGetTranForward method
-*     inherited from the Mapping class).
-
-*  Description:
-*     This function returns a value indicating if the PolyMap is able
-*     to perform a forward coordinate transformation.
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     Zero if the forward coordinate transformation is not defined, or 1 if it
-*     is.
-
-*  Notes:
-*     -  A value of zero will be returned if this function is invoked with the
-*     global error status set, or if it should fail for any reason.
-*/
-
-/* Local Variables: */
-   AstPolyMap *map;            /* Pointer to PolyMap to be queried */
-   int result;                 /* The returned value */
-
-/* Check the global error status. */
-   if ( !astOK ) return 0;
-
-/* Obtain a pointer to the PolyMap. */
-   map = (AstPolyMap *) this;
-
-/* First deal with cases where the PolyMap has not been inverted. */
-   if( !astGetInvert( this ) ) {
-
-/* The PolyMap has a defined forward transformation if one or more
-   coefficients values were supplied for the original forward
-   transformation. It is not possible to replace the original forward
-   transformation with an iterative algorithm. */
-      result = map->ncoeff_f ? 1 : 0;
-
-/* Now deal with cases where the PolyMap has been inverted. */
-   } else {
-
-/* The PolyMap has a defined forward transformation if one or more
-   coefficients values were supplied for the original inverse
-   transformation, or if the original inverse transformation is being
-   approximated using an iterative algorithm. */
-      result = ( map->ncoeff_i || astGetIterInverse( map ) ) ? 1 : 0;
-   }
-
-/* Return the result. */
-   return result;
-}
-
-static int GetTranInverse( AstMapping *this, int *status ) {
-/*
-*
-*  Name:
-*     GetTranInverse
-
-*  Purpose:
-*     Determine if a PolyMap defines an inverse coordinate transformation.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     #include "matrixmap.h"
-*     int GetTranInverse( AstMapping *this, int *status )
-
-*  Class Membership:
-*     PolyMap member function (over-rides the astGetTranInverse method
-*     inherited from the Mapping class).
-
-*  Description:
-*     This function returns a value indicating if the PolyMap is able
-*     to perform an inverse coordinate transformation.
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     Zero if the inverse coordinate transformation is not defined, or 1 if it
-*     is.
-
-*  Notes:
-*     -  A value of zero will be returned if this function is invoked with the
-*     global error status set, or if it should fail for any reason.
-*/
-
-/* Local Variables: */
-   AstPolyMap *map;            /* Pointer to PolyMap to be queried */
-   int result;                 /* The returned value */
-
-/* Check the global error status. */
-   if ( !astOK ) return 0;
-
-/* Obtain a pointer to the PolyMap. */
-   map = (AstPolyMap *) this;
-
-/* First deal with cases where the PolyMap has not been inverted. */
-   if( !astGetInvert( this ) ) {
-
-/* The PolyMap has a defined inverse transformation if one or more
-   coefficients values were supplied for the original inverse
-   transformation, or if the original inverse transformation is being
-   approximated using an iterative algorithm. */
-      result = ( map->ncoeff_i || astGetIterInverse( map ) ) ? 1 : 0;
-
-/* Now deal with cases where the PolyMap has been inverted. */
-   } else {
-
-/* The PolyMap has a defined inverse transformation if one or more
-   coefficients values were supplied for the original forward
-   transformation. It is not possible to replace the original forward
-   transformation with an iterative algorithm. */
-      result = map->ncoeff_f ? 1 : 0;
-   }
-
-/* Return the result. */
-   return result;
-}
-
-void astInitPolyMapVtab_(  AstPolyMapVtab *vtab, const char *name, int *status ) {
-/*
-*+
-*  Name:
-*     astInitPolyMapVtab
-
-*  Purpose:
-*     Initialise a virtual function table for a PolyMap.
-
-*  Type:
-*     Protected function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     void astInitPolyMapVtab( AstPolyMapVtab *vtab, const char *name )
-
-*  Class Membership:
-*     PolyMap vtab initialiser.
-
-*  Description:
-*     This function initialises the component of a virtual function
-*     table which is used by the PolyMap class.
-
-*  Parameters:
-*     vtab
-*        Pointer to the virtual function table. The components used by
-*        all ancestral classes will be initialised if they have not already
-*        been initialised.
-*     name
-*        Pointer to a constant null-terminated character string which contains
-*        the name of the class to which the virtual function table belongs (it
-*        is this pointer value that will subsequently be returned by the Object
-*        astClass function).
-*-
-*/
-
-/* Local Variables: */
-   astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
-   AstObjectVtab *object;        /* Pointer to Object component of Vtab */
-   AstMappingVtab *mapping;      /* Pointer to Mapping component of Vtab */
-
-/* Check the local error status. */
-   if ( !astOK ) return;
-
-
-/* Get a pointer to the thread specific global data structure. */
-   astGET_GLOBALS(NULL);
-
-/* Initialize the component of the virtual function table used by the
-   parent class. */
-   astInitMappingVtab( (AstMappingVtab *) vtab, name );
-
-/* Store a unique "magic" value in the virtual function table. This
-   will be used (by astIsAPolyMap) to determine if an object belongs
-   to this class.  We can conveniently use the address of the (static)
-   class_check variable to generate this unique value. */
-   vtab->id.check = &class_check;
-   vtab->id.parent = &(((AstMappingVtab *) vtab)->id);
-
-/* Initialise member function pointers. */
-/* ------------------------------------ */
-/* Store pointers to the member functions (implemented here) that provide
-   virtual methods for this class. */
-   vtab->PolyPowers = PolyPowers;
-   vtab->FitPoly1DInit = FitPoly1DInit;
-   vtab->FitPoly2DInit = FitPoly2DInit;
-   vtab->PolyTran = PolyTran;
-   vtab->PolyCoeffs = PolyCoeffs;
-
-   vtab->ClearIterInverse = ClearIterInverse;
-   vtab->GetIterInverse = GetIterInverse;
-   vtab->SetIterInverse = SetIterInverse;
-   vtab->TestIterInverse = TestIterInverse;
-
-   vtab->ClearNiterInverse = ClearNiterInverse;
-   vtab->GetNiterInverse = GetNiterInverse;
-   vtab->SetNiterInverse = SetNiterInverse;
-   vtab->TestNiterInverse = TestNiterInverse;
-
-   vtab->ClearTolInverse = ClearTolInverse;
-   vtab->GetTolInverse = GetTolInverse;
-   vtab->SetTolInverse = SetTolInverse;
-   vtab->TestTolInverse = TestTolInverse;
-
-/* Save the inherited pointers to methods that will be extended, and
-   replace them with pointers to the new member functions. */
-   object = (AstObjectVtab *) vtab;
-   mapping = (AstMappingVtab *) vtab;
-
-   parent_getobjsize = object->GetObjSize;
-   object->GetObjSize = GetObjSize;
-
-#if defined(THREAD_SAFE)
-   parent_managelock = object->ManageLock;
-   object->ManageLock = ManageLock;
-#endif
-
-   parent_clearattrib = object->ClearAttrib;
-   object->ClearAttrib = ClearAttrib;
-   parent_getattrib = object->GetAttrib;
-   object->GetAttrib = GetAttrib;
-   parent_setattrib = object->SetAttrib;
-   object->SetAttrib = SetAttrib;
-   parent_testattrib = object->TestAttrib;
-   object->TestAttrib = TestAttrib;
-
-   parent_transform = mapping->Transform;
-   mapping->Transform = Transform;
-   mapping->GetTranForward = GetTranForward;
-   mapping->GetTranInverse = GetTranInverse;
-
-/* Store replacement pointers for methods which will be over-ridden by
-   new member functions implemented here. */
-   object->Equal = Equal;
-   mapping->MapMerge = MapMerge;
-
-/* Declare the destructor and copy constructor. */
-   astSetDelete( (AstObjectVtab *) vtab, Delete );
-   astSetCopy( (AstObjectVtab *) vtab, Copy );
-
-/* Declare the class dump function. */
-   astSetDump( vtab, Dump, "PolyMap", "Polynomial transformation" );
-
-/* If we have just initialised the vtab for the current class, indicate
-   that the vtab is now initialised, and store a pointer to the class
-   identifier in the base "object" level of the vtab. */
-   if( vtab == &class_vtab ) {
-      class_init = 1;
-      astSetVtabClassIdentifier( vtab, &(vtab->id) );
-   }
-}
-
-static void IterInverse( AstPolyMap *this, AstPointSet *out, AstPointSet *result,
-                         int *status ){
-/*
-*  Name:
-*     IterInverse
-
-*  Purpose:
-*     Use an iterative method to evaluate the original inverse transformation
-*     of a PolyMap at a set of (original) output positions.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     void IterInverse( AstPolyMap *this, AstPointSet *out, AstPointSet *result,
-*                       int *status )
-
-*  Description:
-*     This function transforms a set of PolyMap positions using the original
-*     inverse transformation of the PolyMap (i.e. the Negated attribute
-*     is assumed to be zero). An iterative Newton-Raphson method is used
-*     which only required the original forward transformation of the PolyMap
-*     to be defined.
-
-*  Parameters:
-*     this
-*        The PolyMap.
-*     out
-*        A PointSet holding the positions that are to be transformed using
-*        the original inverse transformation. These corresponds to
-*        outputs of the original (i.e. uninverted) PolyMap
-*     result
-*        A PointSet into which the transformed positions are to be stored.
-*        These corresponds to inputs of the original (i.e. uninverted) PolyMap
-
-*     status
-*        Pointer to the inherited status variable.
-
-*/
-
-/* Local Variables: */
-   AstMapping *lintrunc;
-   AstPointSet *work;
-   AstPointSet **ps_jac;
-   AstPolyMap **jacob;
-   double *vec;
-   double *pb;
-   double **ptr_work;
-   double ***ptr_jac;
-   double *mat;
-   double **ptr_out;
-   double **ptr_in;
-   double *pa;
-   double det;
-   double maxerr;
-   double vlensq;
-   double xlensq;
-   double xx;
-   int *flags;
-   int *iw;
-   int fwd;
-   int icol;
-   int icoord;
-   int ipoint;
-   int irow;
-   int iter;
-   int maxiter;
-   int nconv;
-   int ncoord;
-   int npoint;
-   int sing;
-
-/* Check inherited status */
-   if( !astOK ) return;
-
-/* Check the PolyMap has equal numbers of inputs and outputs. */
-   ncoord = astGetNin( this );
-   if( ncoord != astGetNout( this ) ) {
-      astError( AST__INTER, "astTransform(%s): Supplied %s has unequal numbers"
-                " of inputs and outputs and therefore an iterative inverse "
-                "cannot be used (internal AST Programming errpr).", status,
-                astGetClass(this), astGetClass(this) );
-   }
-
-/* Get information about the Jacobian matrix for the forward polynomial
-   transformation. This matrix is a ncoord X ncoord matrix, in which
-   element (row=I,col=J) is the rate of change of output coord I with
-   respect to input coord J, of the supplied PolyMap's forward transformation.
-   The numerical values of the matrix vary depending on where it is
-   evaluated within the input space of the PolyMap. For this reason, the
-   "jacob" variable holds a vector of "ncoord" PolyMaps. The outputs of
-   each of these PolyMaps corresponds to a single column in the Jacobian
-   matrix. */
-   jacob = GetJacobian( this, status );
-
-/* Get the number of points to be transformed. */
-   npoint = astGetNpoint( out );
-
-/* Get another PointSet to hold intermediate results. */
-   work = astPointSet( npoint, ncoord, " ", status );
-
-/* See if the PolyMap has been inverted.*/
-   fwd = !astGetInvert( this );
-
-/* Get pointers to the data arrays for all PointSets. Note, here "in" and
-   "out" refer to inputs and outputs of the PolyMap (i.e. the forward
-   transformation). These are respectively *outputs* and *inputs* of the
-   inverse transformation. */
-   ptr_in = astGetPoints( result );  /* Returned input positions */
-   ptr_out = astGetPoints( out );    /* Supplied output positions */
-   ptr_work = astGetPoints( work );  /* Work space */
-
-/* Allocate an array of PointSets to hold the elements of the Jacobian
-   matrix. */
-   ptr_jac = astMalloc( sizeof( double ** )*ncoord );
-   ps_jac = astCalloc( ncoord, sizeof( AstPointSet * ) );
-   if( astOK ) {
-      for( icoord = 0; icoord < ncoord; icoord++ ) {
-         ps_jac[ icoord ] = astPointSet( npoint, ncoord, " ", status );
-         ptr_jac[ icoord ] = astGetPoints( ps_jac[ icoord ] );
-      }
-   }
-
-/* Allocate an array to hold flags indicating if each position has
-   converged. Initialise it to hold zero at every element. */
-   flags = astCalloc( npoint, sizeof( int ) );
-
-/* Allocate memory to hold the Jacobian matrix at a single point. */
-   mat = astMalloc( sizeof( double )*ncoord*ncoord );
-
-/* Allocate memory to hold the offset vector. */
-   vec = astMalloc( sizeof( double )*ncoord );
-
-/* Allocate memory to hold work space for palDmat. */
-   iw = astMalloc( sizeof( int )*ncoord );
-
-/* Check pointers can be used safely. */
-   if( astOK ) {
-
-/* Store the initial guess at the required input positions. These are
-   determined by transforming the supplied output positions using the
-   inverse of a linear truncation of the PolyMap's forward
-   transformation. */
-      lintrunc = LinearGuess( this, status );
-      (void) astTransform( lintrunc, out, 0, result );
-      lintrunc = astAnnul( lintrunc );
-
-/* Get the maximum number of iterations to perform. */
-      maxiter = astGetNiterInverse( this );
-
-/* Get the target relative error for the returned input axis values, and
-   square it. */
-      maxerr = astGetTolInverse( this );
-      maxerr *= maxerr;
-
-/* Initialise the number of positions which have reached the required
-   accuracy. */
-      nconv = 0;
-
-/* Loop round doing iterations of a Newton-Raphson algorithm, until
-   all points have achieved the required relative error, or the
-   maximum number of iterations have been performed. */
-      for( iter = 0; iter < maxiter && nconv < npoint && astOK; iter++ ) {
-
-/* Use the original forward transformation of the supplied PolyMap to
-   transform the current guesses at the required input positions into
-   the corresponding output positions. Store the results in the "work"
-   PointSet. */
-         (void) astTransform( this, result, fwd, work );
-
-/* Modify the work PointSet so that it holds the offsets from the output
-   positions produced by the current input position guesses, and the
-   required output positions. */
-         for( icoord = 0; icoord < ncoord; icoord++ ) {
-            pa = ptr_out[ icoord ];
-            pb = ptr_work[ icoord ];
-            for( ipoint = 0; ipoint< npoint; ipoint++,pa++,pb++ ) {
-               if( *pa != AST__BAD && *pb != AST__BAD ){
-                  *pb = *pa - *pb;
-               } else {
-                  *pb = AST__BAD;
-               }
-            }
-         }
-
-/* Evaluate the elements of the Jacobian matrix at the current input
-   position guesses. */
-         for( icoord = 0; icoord < ncoord; icoord++ ) {
-            (void) astTransform( jacob[ icoord ], result, 1, ps_jac[ icoord ] );
-         }
-
-/* For each position, we now invert the matrix equation
-
-    Dy = Jacobian.Dx
-
-   to find a guess at the vector (dx) holding the offsets from the
-   current input positions guesses to their required values. Loop over all
-   points. */
-         for( ipoint = 0; ipoint < npoint; ipoint++ ) {
-
-/* Do not change positions that have already converged. */
-            if( !flags[ ipoint ] ) {
-
-/* Get the numerical values for the elements of the Jacobian matrix at
-   the current point. */
-               pa = mat;
-               for( irow = 0; irow < ncoord; irow++ ) {
-                  for( icol = 0; icol < ncoord; icol++ ) {
-                     *(pa++) = ptr_jac[ icol ][ irow ][ ipoint ];
-                  }
-
-/* Store the offset from the current output position to the required
-   output position. */
-                  vec[ irow ] = ptr_work[ irow ][ ipoint ];
-               }
-
-/* Find the corresponding offset from the current input position to the required
-   input position. */
-               palDmat( ncoord, mat, vec, &det, &sing, iw );
-
-/* If the matrix was singular, the input position cannot be evaluated so
-   store a bad value for it and indicate it has converged. */
-               if( sing ) {
-                  for( icoord = 0; icoord < ncoord; icoord++ ) {
-                     ptr_in[ icoord ][ ipoint ] = AST__BAD;
-                  }
-                  flags[ ipoint ] = 1;
-                  nconv++;
-
-/* Otherwise, update the input position guess. */
-               } else {
-                  vlensq = 0.0;
-                  xlensq = 0.0;
-                  pa = vec;
-                  for( icoord = 0; icoord < ncoord; icoord++,pa++ ) {
-                     xx = ptr_in[ icoord ][ ipoint ] + (*pa);
-                     ptr_in[ icoord ][ ipoint ] = xx;
-                     xlensq += xx*xx;
-                     vlensq += (*pa)*(*pa);
-                  }
-
-/* Check for convergence. */
-                  if( vlensq <= maxerr*xlensq ) {
-                     flags[ ipoint ] = 1;
-                     nconv++;
-                  }
-               }
-            }
-         }
-      }
-   }
-
-/* Free resources. */
-   vec = astFree( vec );
-   iw = astFree( iw );
-   mat = astFree( mat );
-   flags = astFree( flags );
-   work = astAnnul( work );
-
-   if( ps_jac ) {
-      for( icoord = 0; icoord < ncoord; icoord++ ) {
-         ps_jac[ icoord ] = astAnnul( ps_jac[ icoord ] );
-      }
-      ps_jac = astFree( ps_jac );
-   }
-
-   ptr_jac = astFree( ptr_jac );
-
-}
-
-static AstMapping *LinearGuess( AstPolyMap *this, int *status ){
-/*
-*  Name:
-*     LinearTruncation
-
-*  Purpose:
-*     Get a Mapping representing a linear approximation of a PolyMap
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     AstMapping *LinearGuess( AstPolyMap *this, int *status )
-
-*  Description:
-*     This function returns a linear Mapping approximating the original
-*     forward transformation of the supplied PolyMap (i.e. the Invert
-*     flag is assume to be zero). The linear and constant terms in the
-*     supplied PolyMap are used for all outputs that have such terms. Any
-*     other outputs are assumed to be equal to the corresponding inputs.
-*     If the forward transformation of the PolyMap is defined, then it is
-*     used to determine the returned Mapping. Otherwise, the inverse
-*     transformation of the PolyMap is used. The forward transformation of
-*     the returned Mapping always represents the forward transformation of
-*     the original (i.e. uninverted) PolyMap.
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     The returned linear Mapping, or NULL if an error occurs.
-
-*/
-
-/* Local Variables: */
-   AstMapping *result;           /* Pointer to returned Mapping */
-   AstMapping *mm;               /* MatrixMap representing linear terms */
-   AstMapping *sm;               /* ShiftMap representing offset terms */
-   double **coeff;               /* Pointer to coefficient value arrays */
-   double *matrix;               /* Linear scaling matrix */
-   double *vector;               /* Offset vector */
-   int ***power;                 /* Pointer to coefficient power arrays */
-   int *ncoeff;                  /* Pointer to no. of coefficients */
-   int flag;                     /* Indicates nature of current coeff */
-   int ico;                      /* Coefficient index */
-   int iin;                      /* Index of transformation input */
-   int inverse;                  /* PolyMap inverse transformation selected? */
-   int iout;                     /* Index of transformation output */
-   int jin;                      /* Transformation input for current coeff */
-   int nin;                      /* No. of i/ps for selected transformation */
-   int nout;                     /* No. of o/ps for selected transformation */
-   int ok;                       /* Current output is defined? */
-   int pow;                      /* Power for current coeff and input */
-
-/* Initialise returned value */
-   result = NULL;
-
-/* Check inherited status */
-   if( !astOK ) return result;
-
-/* If the required Mapping has already been determined, return a pointer
-   to it. */
-   if( this->lintrunc ) return astClone( this->lintrunc );
-
-/* Get pointers to the arrays holding the required coefficient
-   values and powers. Use the forward transformation of the original
-   (uninverted) PolyMap if it is defined and the inverse transformation
-   otherwise. */
-   if ( this->ncoeff_f ){
-      ncoeff = this->ncoeff_f;
-      coeff = this->coeff_f;
-      power = this->power_f;
-      nin = astGetNin( this );
-      nout = astGetNout( this );
-      inverse = 0;
-
-   } else {
-      ncoeff = this->ncoeff_i;
-      coeff = this->coeff_i;
-      power = this->power_i;
-      nout = astGetNin( this );
-      nin = astGetNout( this );
-      inverse = 1;
-   }
-
-/* Allocate memory to hold a matrix holding the linear terms of the
-   PolyMap transformation selected above. Initialise it to hold zeros. */
-   matrix = astCalloc( nin*nout, sizeof( double * ) );
-
-/* Allocate memory to hold a vector holding the offset terms of the
-   PolyMap transformation selected above. Initialise it to hold zeros. */
-   vector = astCalloc( nout, sizeof( double * ) );
-   if( astOK ){
-
-/* Loop round all the outputs of the transformation selected above. */
-
-      for( iout = 0; iout < nout; iout++ ){
-
-/* Loop round all the coefficients for the current transformation output. */
-         for( ico = 0; ico < ncoeff[iout]; ico++ ){
-
-/* Initialise a tristate flag that indicates if the current coefficient is
-   constant (0), linear (1) or something else (2). */
-            flag = 0;
-
-/* Look for the first input axis for which the current coeff uses a
-   non-zero power. If the power is not one, the current coeff is neither
-   linear nor constant, and so we can pass on to the next coeff. If any
-   subsequent non-zero power is found, the coeff is not linear. */
-            jin = 0;
-            for( iin = 0; iin < nin; iin++ ){
-               pow = power[iout][ico][iin];
-               if( pow ) {
-                  if( pow != 1 ) {
-                     flag = 2;
-                     break;
-                  } else if( flag == 1 ) {
-                     flag = 2;
-                     break;
-                  } else {
-                     flag = 1;
-                     jin = iin;
-                  }
-               }
-            }
-
-/* If the coeff is a constant term, insert it into the offset vector. */
-            if( flag == 0 ) {
-               vector[ iout ] = coeff[ iout ][ ico ];
-
-/* If the coeff is a linear term, insert it into the matrix. */
-            } else if( flag == 1 ) {
-               matrix[ jin + iout*nin ] = coeff[ iout ][ ico ];
-            }
-         }
-      }
-
-/* If any PolyMap output has no linear terms, the matrix will not be
-   invertable. So insert a 1.0 value on the diagonal of such rows. */
-      for( iout = 0; iout < nout; iout++ ){
-         ok = 0;
-         for( iin = 0; iin < nin; iin++ ) {
-            if( matrix[ iin + iout*nin ] != 0.0 ) {
-               ok = 1;
-               break;
-            }
-         }
-         if( ! ok ) matrix[ iout + iout*nin ] = 1.0;
-      }
-
-/* Create a MatrixMap to represent the matrix. */
-      mm = (AstMapping *) astMatrixMap( nin, nout, 0, matrix, " ", status );
-
-/* If the MatrixMap cannot be inverted, use a unit map instead. */
-      if( !astGetTranInverse( mm ) ) {
-         mm = astAnnul( mm );
-         mm = (AstMapping *) astUnitMap( nin, " ", status );
-      }
-
-/* Create a ShiftMap to represent the offset. */
-      sm = (AstMapping *) astShiftMap( nout, vector, " ", status );
-
-/* Combine them in series into a single compound Mapping. */
-      result = (AstMapping *) astCmpMap( mm, sm, 1, " ", status );
-
-/* If this CmpMap represents the inverse transformation of the original
-   (i.e. uninverted) PolyMap, invert it so that it represents the forward
-   transformation. */
-      if( inverse ) astInvert( result );
-
-/* Free resources */
-      mm = astAnnul( mm );
-      sm = astAnnul( sm );
-
-/* Store the Mapping in the PolyMap structure so we do not need to
-   recalculate it every time it is needed. */
-      if( astOK ) this->lintrunc = astClone( result );
-   }
-
-/* Free resources */
-   matrix = astFree( matrix );
-   vector = astFree( vector );
-
-/* Return the required Mapping. */
-   return result;
-
-}
-
-static void LMFunc1D( const double *p, double *hx, int m, int n, void *adata ){
-/*
-*  Name:
-*     LMFunc1D
-
-*  Purpose:
-*     Evaluate a test 1D polynomal.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     void LMFunc1D( const double *p, double *hx, int m, int n, void *adata )
-
-*  Description:
-*     This function finds the residuals implied by a supplied set of
-*     candidate polynomial coefficients. Each residual is a candidate
-*     polynomial evaluated at one of the sample points, minus the
-*     supplied target value for the polynomial at that test point.
-*
-*     The minimisation process minimises the sum of the squared residuals.
-
-*  Parameters:
-*     p
-*        An array of "m" coefficients for the candidate polynomial. The
-*        coefficients are ordered C0, C1, C2, etc.
-*     hx
-*        An array in which to return the "n" residuals. The residual at
-*        sample "k" is returned in element (k).
-*     m
-*        The length of the "p" array. This should be equal to order.
-*     n
-*        The length of the "hx" array. This should be equal to nsamp.
-*     adata
-*        Pointer to a structure holding the sample positions and values,
-*        and other information.
-
-*/
-
-/* Local Variables: */
-   AstMinPackData *data;
-   double *px1;
-   double *py;
-   const double *vp;
-   double *vr;
-   double res;
-   int k;
-   int w1;
-
-/* Get a pointer to the data structure. */
-   data = (AstMinPackData *) adata;
-
-/* Initialise a pointer to the current returned residual value. */
-   vr = hx;
-
-/* Initialise a pointer to the sampled Y values for the polynomial output. */
-   py = data->y[ 0 ];
-
-/* Initialise a pointer to the powers of the input X values at the curent
-   (i.e. first) sample. */
-   px1 = data->xp1;
-
-/* Loop over the index of the sample to which this residual refers. */
-   for( k = 0; k < data->nsamp; k++ ) {
-
-/* Initialise a pointer to the first coefficient (the  constant term) for the
-   polynomial output coordinate. */
-      vp = p;
-
-/* Initialise this residual to hold the sampled Y value. Increment the
-   pointer to the next sampled value for the current polynomial output. */
-      res = -( *(py++) );
-
-/* Loop over the coefficients. */
-      for( w1 = 0; w1 < data->order; w1++ ) {
-
-/* Increment the residual by the value of the current term Cw1*(x1^w1).
-   Increment the pointer to the next coefficient (C). Also increment the
-   pointer to the next higher power of X1. */
-         res += ( *(vp++) )*( *(px1++) );
-      }
-
-/* Store the complete residual in the returned array, and increment the
-   pointer to the next residual. */
-      *(vr++) = res;
-   }
-}
-
-static void LMFunc2D( const double *p, double *hx, int m, int n, void *adata ){
-/*
-*  Name:
-*     LMFunc2D
-
-*  Purpose:
-*     Evaluate a test 2D polynomal.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     void LMFunc2D( const double *p, double *hx, int m, int n, void *adata )
-
-*  Description:
-*     This function finds the residuals implied by a supplied set of
-*     candidate polynomial coefficients. Each residual is a candidate
-*     polynomial (either P1 or P2) evaluated at one of the sample points
-*     (x1,x2), minus the supplied target value for the polynomial at that
-*     test point.
-*
-*     The minimisation process minimises the sum of the squared residuals.
-
-*  Parameters:
-*     p
-*        An array of "m" coefficients for the candidate polynomial. All the
-*        coefficients for polynomial P1 come first, followed by those for P2.
-*        Within each each polynomial, the coefficients are order C00, C10,
-*        C01, C20, C11, C02, C30, C21, C12, C03, etc. So the coefficient
-*        of (x1^j*x2^k) (=Cjk) for polynomial Pi is stored in element
-*        [k + (j + k)*(j + k + 1)/2 + i*order*(order+1)/2] of the "p" array.
-*     hx
-*        An array in which to return the "n" residuals. The residual at
-*        sample "k" for polynomial "i" is returned in element (k + nsamp*i).
-*     m
-*        The length of the "p" array. This should be equal to order*(order+1).
-*     n
-*        The length of the "hx" array. This should be equal to 2*nsamp.
-*     adata
-*        Pointer to a structure holding the sample positions and values,
-*        and other information.
-
-*/
-
-/* Local Variables: */
-   AstMinPackData *data;
-   const double *vp0;
-   const double *vp;
-   double *py;
-   double *px1;
-   double *px10;
-   double *px20;
-   double *vr;
-   double res;
-   double *px2;
-   int iout;
-   int k;
-   int w12;
-   int w2;
-
-/* Get a pointer to the data structure. */
-   data = (AstMinPackData *) adata;
-
-/* Initialise a pointer to the current returned residual value. */
-   vr = hx;
-
-/* Initialise a pointer to the first coefficient (the  constant term) for the
-   current (i.e. first) polynomial output  coordinate. */
-   vp0 = p;
-
-/* Loop over each polynomial output coordinate. */
-   for( iout = 0; iout < 2; iout++ ) {
-
-/* Initialise a pointer to the sampled Y values for the first polynomial
-   output. */
-      py = data->y[ iout ];
-
-/* Initialise pointers to the powers of the input X values at the curent
-   (i.e. first) sample. */
-      px10 = data->xp1;
-      px20 = data->xp2;
-
-/* Loop over the index of the sample to which this residual refers. */
-      for( k = 0; k < data->nsamp; k++ ) {
-
-/* Reset the pointer to the first coefficient (the  constant term)
-   for the current polynomial output  coordinate. */
-         vp = vp0;
-
-/* Initialise this residual to hold the sampled Y value. Increment the
-   pointer to the next sampled value for the current polynomial output. */
-         res = -( *(py++) );
-
-/* The w12 value is the sum of the powers of X1 and X2. So w12=0
-   corresponds to the constant term in the polynomial, and (e.g.) w12=6
-   corresponds to all the terms for which the sum of the powerss (w1+w2)
-   is 6. Loop over all possible w12 values. */
-         for( w12 = 0; w12 < data->order; w12++ ) {
-
-/* The next coeff refers to (x1^w12)*(x2^0). Get pointers to the values
-   holding x1^w12 and x2^0. */
-            px1 = px10++;
-            px2 = px20;
-
-/* Loop over powers of x2. The corresponding power of x1 is "w12-w2", but
-   is not explicitly needed here. So x1 moves down from w12 to zero, as
-   w2 moves up from zero to w12. */
-            for( w2 = 0; w2 <= w12; w2++ ) {
-
-/* Increment the residual by the value of the current term Cw1w2*(x1^w1)*(x2^w2).
-   Increment the pointer tio the next coefficient (C). Also decrement the
-   pointer to the next lower power of X1, and increment the pointer to the next
-   higher power of X2. */
-               res += ( *(vp++) )*( *(px1--) )*( *(px2++) );
-            }
-         }
-
-/* Move on to the x2 powers for the next sample. Don't need to do this
-   for x1 since px10 is incremented within the above loop. */
-         px20 += data->order;
-
-/* Store the complete residual in the returned array, and increment the
-   pointer to the next residual. */
-         *(vr++) = res;
-      }
-
-/* Get a pointer to the first coefficient (the  constant term) for the
-   next polynomial output  coordinate. */
-      vp0 += data->order*( 1 + data->order )/2;
-   }
-}
-
-static void LMJacob1D( const double *p, double *jac, int m, int n, void *adata ){
-/*
-*  Name:
-*     LMJacob1D
-
-*  Purpose:
-*     Evaluate the Jacobian matrix of a test 1D polynomal.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     void LMJacob1D( const double *p, double *jac, int m, int n, void *adata )
-
-*  Description:
-*     This function finds the Jacobian matrix that describes the rate of
-*     change of every residual with respect to every polynomial coefficient.
-*     Each residual is a candidate polynomial evaluated at one of the sample
-*     points minus the supplied target value for the polynomial at that test
-*     point.
-*
-*     For a polynomial the Jacobian matrix is constant (i.e. does not
-*     depend on the values of the polynomial coefficients). So we only
-*     evaluate it on the first call.
-
-*  Parameters:
-*     p
-*        An array of "m" coefficients for the candidate polynomial.
-*     jac
-*        An array in which to return the "m*n" elements of the Jacobian
-*        matrix. The rate of change of residual "r" with respect to
-*        coefficient "c" is returned in element "r + c*n". The residual
-*        at sample "k" of polynomial Pi has an "r" index of (k + nsamp*i).
-*        The coefficient of (x1^j) for polynomial Pi has a "c" index
-*        of j.
-*     m
-*        The length of the "p" array. This should be equal to order.
-*     n
-*        The number of residuals. This should be equal to nsamp.
-*     adata
-*        Pointer to a structure holding the sample positions and values,
-*        and other information.
-
-*/
-
-/* Local Variables: */
-   AstMinPackData *data;
-   int k;
-   int w1;
-
-/* Get a pointer to the data structure. */
-   data = (AstMinPackData *) adata;
-
-/* The Jacobian of the residuals with respect to the polynomial
-   coefficients is constant (i.e. does not depend on the values of the
-   polynomial coefficients). So we only need to calculate it once. If
-   this is the first call, calculate the Jacobian and return it in "jac".
-   otherwise, just return immediately retaining the supplied "jac" values
-   (which will be the values returned by the previous call to this
-   function). */
-   if( data->init_jac ) {
-      data->init_jac = 0;
-
-/* Loop over all residuals. */
-      for( k = 0; k < n; k++ ) {
-
-/* Loop over all parameters (i.e. polynomial coefficients). */
-         for( w1 = 0; w1 < m; w1++ ) {
-
-/* Store the Jacobian. */
-            jac[ k + w1*n ] = (data->xp1)[ w1 + k*data->order ];
-         }
-      }
-   }
-}
-
-static void LMJacob2D( const double *p, double *jac, int m, int n, void *adata ){
-/*
-*  Name:
-*     LMJacob2D
-
-*  Purpose:
-*     Evaluate the Jacobian matrix of a test 2D polynomal.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     void LMJacob2D( const double *p, double *jac, int m, int n, void *adata )
-
-*  Description:
-*     This function finds the Jacobian matrix that describes the rate of
-*     change of every residual with respect to every polynomial coefficient.
-*     Each residual is a candidate polynomial (either P1 or P2) evaluated
-*     at one of the sample points (x1,x2), minus the supplied target value
-*     for the polynomial at that test point.
-*
-*     For a polynomial the Jacobian matrix is constant (i.e. does not
-*     depend on the values of the polynomial coefficients). So we only
-*     evaluate it on the first call.
-
-*  Parameters:
-*     p
-*        An array of "m" coefficients for the candidate polynomial. All the
-*        coefficients for polynomial P1 come first, followed by those for P2.
-*        Within each each polynomial, the coefficients are order C00, C10,
-*        C01, C20, C11, C02, C30, C21, C12, C03, etc. So the coefficient
-*        of (x1^j*x2^k) (=Cjk) for polynomial Pi is stored in element
-*        [k + (j + k)*(j + k + 1)/2 + i*order*(order+1)/2] of the "p" array.
-*     jac
-*        An array in which to return the "m*n" elements of the Jacobian
-*        matrix. The rate of change of residual "r" with respect to
-*        coefficient "c" is returned in element "r + c*n". The residual
-*        at sample "k" of polynomial Pi has an "r" index of (k + nsamp*i).
-*        The coefficient of (x1^j*x2^k) for polynomial Pi has a "c" index
-*        of [k + (j + k)*(j + k + 1)/2 + i*order*(order+1)/2].
-*     m
-*        The length of the "p" array. This should be equal to order*(order+1).
-*     n
-*        The number of residuals. This should be equal to 2*nsamp.
-*     adata
-*        Pointer to a structure holdin gthe sample positions and values,
-*        and other information.
-
-*/
-
-/* Local Variables: */
-   AstMinPackData *data;
-   int iout;
-   int k;
-   int ncof;
-   int vp;
-   int vr;
-   int w1;
-   int w12;
-   int w2;
-
-/* Get a pointer to the data structure. */
-   data = (AstMinPackData *) adata;
-
-/* The Jacobian of the residuals with respect to the polynomial
-   coefficients is constant (i.e. does not depend on the values of the
-   polynomial coefficients). So we only need to calculate it once. If
-   this is the first call, calculate the Jacobian and return it in "jac".
-   otherwise, just return immediately retaining the supplied "jac" values
-   (which will be the values returned by the previous call to this
-   function). */
-   if( data->init_jac ) {
-      data->init_jac = 0;
-
-/* Store the number of coefficients in one polynomial. */
-      ncof = data->order*( 1 + data->order )/2;
-
-/* Loop over all residuals. */
-      for( vr = 0; vr < n; vr++ ) {
-
-/* Calculate the polynomial output index, and sample index, that creates
-   the current residual. */
-         iout = vr/data->nsamp;
-         k = vr - iout*data->nsamp;
-
-/* Loop over all parameters (i.e. polynomial coefficients). */
-         for( vp = 0; vp < m; vp++ ) {
-
-/* If this coefficient is not used in the creation of the current
-   polynomial output value, then the Jacobian value is zero. */
-            if( vp/ncof != iout ) {
-               jac[ vr + vp*n ] = 0.0;
-
-/* Otherwise, get the powers of the two polynomial inputs, to which
-   the current coefficient relates. */
-            } else {
-               w12 = ( -1.0 + sqrt( 1.0 + 8.0*(vp - iout*ncof) ) )/2.0;
-               w2 = vp - iout*ncof - w12*( w12 + 1 )/2;
-               w1 = w12 - w2;
-
-/* Store the Jacobian. */
-               jac[ vr + vp*n ] = (data->xp1)[ w1 + k*data->order ]*
-                                  (data->xp2)[ w2 + k*data->order ];
-            }
-         }
-      }
-   }
-}
-
-#if defined(THREAD_SAFE)
-static int ManageLock( AstObject *this_object, int mode, int extra,
-                       AstObject **fail, int *status ) {
-/*
-*  Name:
-*     ManageLock
-
-*  Purpose:
-*     Manage the thread lock on an Object.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     #include "object.h"
-*     AstObject *ManageLock( AstObject *this, int mode, int extra,
-*                            AstObject **fail, int *status )
-
-*  Class Membership:
-*     PolyMap member function (over-rides the astManageLock protected
-*     method inherited from the parent class).
-
-*  Description:
-*     This function manages the thread lock on the supplied Object. The
-*     lock can be locked, unlocked or checked by this function as
-*     deteremined by parameter "mode". See astLock for details of the way
-*     these locks are used.
-
-*  Parameters:
-*     this
-*        Pointer to the Object.
-*     mode
-*        An integer flag indicating what the function should do:
-*
-*        AST__LOCK: Lock the Object for exclusive use by the calling
-*        thread. The "extra" value indicates what should be done if the
-*        Object is already locked (wait or report an error - see astLock).
-*
-*        AST__UNLOCK: Unlock the Object for use by other threads.
-*
-*        AST__CHECKLOCK: Check that the object is locked for use by the
-*        calling thread (report an error if not).
-*     extra
-*        Extra mode-specific information.
-*     fail
-*        If a non-zero function value is returned, a pointer to the
-*        Object that caused the failure is returned at "*fail". This may
-*        be "this" or it may be an Object contained within "this". Note,
-*        the Object's reference count is not incremented, and so the
-*        returned pointer should not be annulled. A NULL pointer is
-*        returned if this function returns a value of zero.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*    A local status value:
-*        0 - Success
-*        1 - Could not lock or unlock the object because it was already
-*            locked by another thread.
-*        2 - Failed to lock a POSIX mutex
-*        3 - Failed to unlock a POSIX mutex
-*        4 - Bad "mode" value supplied.
-
-*  Notes:
-*     - This function attempts to execute even if an error has already
-*     occurred.
-*/
-
-/* Local Variables: */
-   AstPolyMap *this;
-   int ic;
-   int result;
-   int nc;
-
-/* Initialise */
-   result = 0;
-
-/* Check the supplied pointer is not NULL. */
-   if( !this_object ) return result;
-
-/* Obtain a pointer to the PolyMap structure. */
-   this = (AstPolyMap *) this_object;
-
-/* Invoke the ManageLock method inherited from the parent class. */
-   if( !result ) result = (*parent_managelock)( this_object, mode, extra,
-                                                fail, status );
-
-/* Invoke the astManageLock method on any Objects contained within
-   the supplied Object. */
-   if( this->jacobian ) {
-      nc = astGetNin( this );
-      for( ic = 0; ic < nc && result; ic++ ) {
-         result = astManageLock( (this->jacobian)[ ic ], mode,
-                                 extra, fail );
-      }
-   }
-
-   return result;
-
-}
-#endif
-
-static int MapMerge( AstMapping *this, int where, int series, int *nmap,
-                     AstMapping ***map_list, int **invert_list, int *status ) {
-/*
-*  Name:
-*     MapMerge
-
-*  Purpose:
-*     Simplify a sequence of Mappings containing a PolyMap.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     #include "mapping.h"
-*     int MapMerge( AstMapping *this, int where, int series, int *nmap,
-*                   AstMapping ***map_list, int **invert_list, int *status )
-
-*  Class Membership:
-*     PolyMap method (over-rides the protected astMapMerge method
-*     inherited from the Mapping class).
-
-*  Description:
-*     This function attempts to simplify a sequence of Mappings by
-*     merging a nominated PolyMap in the sequence with its neighbours,
-*     so as to shorten the sequence if possible.
-*
-*     In many cases, simplification will not be possible and the
-*     function will return -1 to indicate this, without further
-*     action.
-*
-*     In most cases of interest, however, this function will either
-*     attempt to replace the nominated PolyMap with a Mapping which it
-*     considers simpler, or to merge it with the Mappings which
-*     immediately precede it or follow it in the sequence (both will
-*     normally be considered). This is sufficient to ensure the
-*     eventual simplification of most Mapping sequences by repeated
-*     application of this function.
-*
-*     In some cases, the function may attempt more elaborate
-*     simplification, involving any number of other Mappings in the
-*     sequence. It is not restricted in the type or scope of
-*     simplification it may perform, but will normally only attempt
-*     elaborate simplification in cases where a more straightforward
-*     approach is not adequate.
-
-*  Parameters:
-*     this
-*        Pointer to the nominated PolyMap which is to be merged with
-*        its neighbours. This should be a cloned copy of the PolyMap
-*        pointer contained in the array element "(*map_list)[where]"
-*        (see below). This pointer will not be annulled, and the
-*        PolyMap it identifies will not be modified by this function.
-*     where
-*        Index in the "*map_list" array (below) at which the pointer
-*        to the nominated PolyMap resides.
-*     series
-*        A non-zero value indicates that the sequence of Mappings to
-*        be simplified will be applied in series (i.e. one after the
-*        other), whereas a zero value indicates that they will be
-*        applied in parallel (i.e. on successive sub-sets of the
-*        input/output coordinates).
-*     nmap
-*        Address of an int which counts the number of Mappings in the
-*        sequence. On entry this should be set to the initial number
-*        of Mappings. On exit it will be updated to record the number
-*        of Mappings remaining after simplification.
-*     map_list
-*        Address of a pointer to a dynamically allocated array of
-*        Mapping pointers (produced, for example, by the astMapList
-*        method) which identifies the sequence of Mappings. On entry,
-*        the initial sequence of Mappings to be simplified should be
-*        supplied.
-*
-*        On exit, the contents of this array will be modified to
-*        reflect any simplification carried out. Any form of
-*        simplification may be performed. This may involve any of: (a)
-*        removing Mappings by annulling any of the pointers supplied,
-*        (b) replacing them with pointers to new Mappings, (c)
-*        inserting additional Mappings and (d) changing their order.
-*
-*        The intention is to reduce the number of Mappings in the
-*        sequence, if possible, and any reduction will be reflected in
-*        the value of "*nmap" returned. However, simplifications which
-*        do not reduce the length of the sequence (but improve its
-*        execution time, for example) may also be performed, and the
-*        sequence might conceivably increase in length (but normally
-*        only in order to split up a Mapping into pieces that can be
-*        more easily merged with their neighbours on subsequent
-*        invocations of this function).
-*
-*        If Mappings are removed from the sequence, any gaps that
-*        remain will be closed up, by moving subsequent Mapping
-*        pointers along in the array, so that vacated elements occur
-*        at the end. If the sequence increases in length, the array
-*        will be extended (and its pointer updated) if necessary to
-*        accommodate any new elements.
-*
-*        Note that any (or all) of the Mapping pointers supplied in
-*        this array may be annulled by this function, but the Mappings
-*        to which they refer are not modified in any way (although
-*        they may, of course, be deleted if the annulled pointer is
-*        the final one).
-*     invert_list
-*        Address of a pointer to a dynamically allocated array which,
-*        on entry, should contain values to be assigned to the Invert
-*        attributes of the Mappings identified in the "*map_list"
-*        array before they are applied (this array might have been
-*        produced, for example, by the astMapList method). These
-*        values will be used by this function instead of the actual
-*        Invert attributes of the Mappings supplied, which are
-*        ignored.
-*
-*        On exit, the contents of this array will be updated to
-*        correspond with the possibly modified contents of the
-*        "*map_list" array.  If the Mapping sequence increases in
-*        length, the "*invert_list" array will be extended (and its
-*        pointer updated) if necessary to accommodate any new
-*        elements.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     If simplification was possible, the function returns the index
-*     in the "map_list" array of the first element which was
-*     modified. Otherwise, it returns -1 (and makes no changes to the
-*     arrays supplied).
-
-*  Notes:
-*     - A value of -1 will be returned if this function is invoked
-*     with the global error status set, or if it should fail for any
-*     reason.
-*/
-
-/* Local Variables: */
-   AstPolyMap *pmap0;    /* Nominated PolyMap */
-   AstPolyMap *pmap1;    /* Neighbouring PolyMap */
-   int i;                /* Index of neighbour */
-   int nin;              /* Number of input coordinates for nominated PolyMap */
-   int nout;             /* Number of output coordinates for nominated PolyMap */
-   int ok;               /* Are PolyMaps equivalent? */
-   int oldinv0;          /* Original Invert value in pmap0 */
-   int oldinv1;          /* Original Invert value in pmap1 */
-   int result;           /* Result value to return */
-
-/* Initialise. */
-   result = -1;
-
-/* Check the global error status. */
-   if ( !astOK ) return result;
-
-/* Save a pointer to the nominated PolyMap. */
-   pmap0 = (AstPolyMap *) ( *map_list )[ where ];
-
-/* The only simplification which can currently be performed is to merge a PolyMap
-   with its own inverse. This can only be done in series. Obviously,
-   there are potentially other simplications which could be performed, but
-   time does not currently allow these to be coded. */
-   if( series ) {
-
-/* Temporarily set the Invert flag of the nominated PolyMap to the
-   required value, first saving the original value so that it can be
-   re-instated later. */
-      oldinv0 = astGetInvert( pmap0 );
-      astSetInvert( pmap0, ( *invert_list )[ where ] );
-
-/* Get the number of inputs and outputs to the nominated PolyMap. */
-      nin = astGetNin( pmap0 );
-      nout = astGetNout( pmap0 );
-
-/* Check each neighbour. */
-      for( i = where - 1; i <= where + 1; i += 2 ) {
-
-/* Continue with the next pass if the neighbour does not exist. */
-         if( i < 0 || i >= *nmap ) continue;
-
-/* Continue with the next pass if this neighbour is not a PermMap. */
-         if( ! astIsAPolyMap( ( *map_list )[ i ] ) ) continue;
-
-/* Get a pointer to it. */
-         pmap1 = (AstPolyMap *) ( *map_list )[ i ];
-
-/* Check it is used in the opposite direction to the nominated PolyMap. */
-         if( ( *invert_list )[ i ] == ( *invert_list )[ where ] ) continue;
-
-/* We use the astEqual method to check that the two PolyMaps are equal.
-   But at the moment they may not be equal because they may have
-   different Invert flags. Therefore, temporarily set the invert flag
-   of the neighbour so that it is the same as the nominated PolyMap,
-   first saving the original value so that it can be re-instated later.
-   Note, we have already checked that the two PolyMaps are used in opposite
-   directions within the CmpMap. */
-         oldinv1 = astGetInvert( pmap1 );
-         astSetInvert( pmap1, ( *invert_list )[ where ] );
-
-/* Use astEqual to check that the coefficients etc are equal in the two
-   PolyMaps. */
-         ok = astEqual( pmap0, pmap1 );
-
-/* Re-instate the original value of the Invert flag in the neighbour. */
-         astSetInvert( pmap1, oldinv1 );
-
-/* Pass on to the next neighbour if the current neighbour differs from
-   the nominated PolyMap. */
-         if( !ok ) continue;
-
-/* If we get this far, then the nominated PolyMap and the current
-   neighbour cancel each other out, so replace each by a UnitMap. */
-         pmap0 = astAnnul( pmap0 );
-         pmap1 = astAnnul( pmap1 );
-         if( i < where ) {
-            ( *map_list )[ where ] = (AstMapping *) astUnitMap( nout, "", status );
-            ( *map_list )[ i ] = (AstMapping *) astUnitMap( nout, "", status );
-            ( *invert_list )[ where ] = 0;
-            ( *invert_list )[ i ] = 0;
-            result = i;
-         } else {
-            ( *map_list )[ where ] = (AstMapping *) astUnitMap( nin, "", status );
-            ( *map_list )[ i ] = (AstMapping *) astUnitMap( nin, "", status );
-            ( *invert_list )[ where ] = 0;
-            ( *invert_list )[ i ] = 0;
-            result = where;
-         }
-
-/* Leave the loop. */
-         break;
-      }
-
-/* If the nominated PolyMap was not replaced by a UnitMap, then
-   re-instate its original value for the Invert flag. */
-      if( pmap0 ) astSetInvert( pmap0, oldinv0 );
-   }
-
-/* Return the result. */
-   return result;
-}
-
-static int MPFunc1D( void *p, int m, int n, const double *x, double *fvec,
-                     double *fjac, int ldfjac, int iflag ) {
-/*
-*  Name:
-*     MPFunc1D
-
-*  Purpose:
-*     Evaluate a test 1D polynomal or its Jacobian.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     int MPFunc1D( void *p, int m, int n, const double *x, double *fvec,
-*                   double *fjac, int ldfjac, int iflag )
-
-*  Description:
-*     This function finds the residuals or Jacobian implied by a supplied
-*     set of candidate polynomial coefficients. Each residual is a candidate
-*     polynomial evaluated at one of the sample points, minus the
-*     supplied target value for the polynomial at that test point.
-*
-*     The minimisation process minimises the sum of the squared residuals.
-
-*  Parameters:
-*     p
-*        A pointer to data needed to evaluate the required residuals or
-*        Jacobians.
-*     m
-*        The number of residuals.
-*     n
-*        The number of coefficients.
-*     x
-*        An array of "n" coefficients for the candidate polynomial.
-*     fvec
-*        An array in which to return the "m" residuals.
-*     fjac
-*        An array in which to return the "mxn" Jacobian values.
-*     ldflac
-*        Unused (should always be equal to "m").
-*     iflag
-*        If 1 calculate the residuals, otherwise calculate the Jacobians.
-
-*/
-
-/* If required, calculate the function values at X, and return this
-   vector in fvec. Do not alter fjac. */
-   if( iflag == 1 ) {
-      LMFunc1D( x, fvec, n, m, p );
-
-/* Otherwise, calculate the Jacobian values at X, and return this
-   matrix in fjac. Do not alter fvec. */
-   } else {
-      LMJacob1D( x, fjac, n, m, p );
-   }
-
-   return 0;
-}
-
-static int MPFunc2D( void *p, int m, int n, const double *x, double *fvec,
-                     double *fjac, int ldfjac, int iflag ) {
-/*
-*  Name:
-*     MPFunc1D
-
-*  Purpose:
-*     Evaluate a test 2D polynomal or its Jacobian.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     int MPFunc2D( void *p, int m, int n, const double *x, double *fvec,
-*                   double *fjac, int ldfjac, int iflag )
-
-*  Description:
-*     This function finds the residuals or Jacobian implied by a supplied
-*     set of candidate polynomial coefficients. Each residual is a candidate
-*     polynomial evaluated at one of the sample points, minus the
-*     supplied target value for the polynomial at that test point.
-*
-*     The minimisation process minimises the sum of the squared residuals.
-
-*  Parameters:
-*     p
-*        A pointer to data needed to evaluate the required residuals or
-*        Jacobians.
-*     m
-*        The number of residuals.
-*     n
-*        The number of coefficients.
-*     x
-*        An array of "n" coefficients for the candidate polynomial.
-*     fvec
-*        An array in which to return the "m" residuals.
-*     fjac
-*        An array in which to return the "mxn" Jacobian values.
-*     ldflac
-*        Unused (should always be equal to "m").
-*     iflag
-*        If 1 calculate the residuals, otherwise calculate the Jacobians.
-
-*/
-
-/* If required, calculate the function values at X, and return this
-   vector in fvec. Do not alter fjac. */
-   if( iflag == 1 ) {
-      LMFunc2D( x, fvec, n, m, p );
-
-/* Otherwise, calculate the Jacobian values at X, and return this
-   matrix in fjac. Do not alter fvec. */
-   } else {
-      LMJacob2D( x, fjac, n, m, p );
-
-   }
-
-   return 0;
-}
-
-static void PolyCoeffs( AstPolyMap *this, int forward, int nel, double *coeffs,
-                        int *ncoeff, int *status ){
-/*
-*++
-*  Name:
-c     astPolyCoeffs
-f     AST_POLYCOEFFS
-
-*  Purpose:
-*     Retrieve the coefficient values used by a PolyMap.
-
-*  Type:
-*     Public function.
-
-*  Synopsis:
-c     #include "polymap.h"
-c     void astPolyCoeffs( AstPolyMap *this, int forward, int nel, double *coeffs,
-c                         int *ncoeff )
-f     CALL AST_POLYCOEFFS( THIS, FORWARD, NEL, COEFFS, NCOEFF, STATUS )
-
-*  Class Membership:
-*     PolyMap method.
-
-*  Description:
-*     This function returns the coefficient values used by either the
-*     forward or inverse transformation of a PolyMap, in the same form
-*     that they are supplied to the PolyMap constructor.
-*
-*     Usually, you should call this method first with
-c     "nel"
-f     NEL
-*     set to zero to determine the number of coefficients used by the
-*     PolyMap. This allows you to allocate an array of the correct size to
-*     hold all coefficient data. You should then call this method a
-*     second time to get the coefficient data.
-
-*  Parameters:
-c     this
-f     THIS = INTEGER (Given)
-*        Pointer to the original Mapping.
-c     forward
-f     FORWARD = LOGICAL (Given)
-c        If non-zero,
-f        If .TRUE.,
-*        the coefficients of the forward PolyMap transformation are
-*        returned. Otherwise the inverse transformation coefficients are
-*        returned.
-c     nel
-f     NEL = INTEGER (Given)
-*        The length of the supplied
-c        "coeffs"
-f        COEFFS
-*        array. It should be at least "ncoeff*( nin + 2 )" if
-c        "foward" is non-zero,
-f        FORWARD is .TRUE.,
-*        and "ncoeff*( nout + 2 )" otherwise, where "ncoeff" is the
-*        number of coefficients to be returned. If a value of zero
-*        is supplied, no coefficient values are returned, but the
-*        number of coefficients used by the transformation is still
-*        returned in
-c        "ncoeff".
-f        NCOEFF.
-c     coeffs
-f     COEFFS( NEL ) = DOUBLE PRECISION (Returned)
-*        An array in which to return the coefficients used by the
-*        requested transformation of the PolyMap. Ignored if
-c        "nel" is zero.
-f        NEL is zero.
-*        The coefficient data is returned in the form in which it is
-*        supplied to the PolyMap constructor. That is, each group of
-*        "2 + nin" or "2 + nout" adjacent elements describe a single
-*        coefficient of the forward or inverse transformation. See the
-*        PolyMap constructor documentation for further details.
-*
-*        If the supplied array is too short to hold all the coefficients,
-*        trailing coefficients are excluded. If the supplied array is
-*        longer than needed to hold all the coefficients, trailing
-*        elements are filled with zeros.
-c     ncoeff
-f     NCOEFF = INTEGER (Returned)
-*        The number of coefficients used by the requested transformation.
-*        A value of zero is returned if the transformation does not
-*        have any defining polynomials. A value is returned for this argument
-*        even if
-c        "nel" is zero.
-f        NEL is zero.
-f     STATUS = INTEGER (Given and Returned)
-f        The global status.
-
-*--
-*/
-
-/* Local Variables: */
-   double **coeff;
-   int ***power;
-   int *nco;
-   int *pp;
-   int iax;
-   int icoeff;
-   int iel;
-   int ipoly;
-   int nax;
-   int npoly;
-
-/* Initialise */
-   *ncoeff = 0;
-
-/* Check the inherited status. */
-   if ( !astOK ) return;
-
-/* Fill any supplied array with zeros. */
-   if( nel ) memset( coeffs, 0, nel*sizeof( *coeffs ) );
-
-/* Get the values to use, taking account of whether the PolyMap has been
-   inverted or not. */
-   if( forward != astGetInvert( this ) ){
-      nco = this->ncoeff_f;
-      power = this->power_f;
-      coeff = this->coeff_f;
-      npoly = astGetNout( this );
-      nax = astGetNin( this );
-   } else {
-      nco = this->ncoeff_i;
-      power = this->power_i;
-      coeff = this->coeff_i;
-      npoly = astGetNin( this );
-      nax = astGetNout( this );
-   }
-
-/* Notheg to do if there are no coeffs. */
-   if( nco && power && coeff ) {
-
-/* Initialise index of next value to store in returned array. */
-      iel = 0;
-
-/* Loop round each 1D polynomial. */
-      for( ipoly = 0; ipoly < npoly; ipoly++ ){
-
-/* Loop round coefficients. */
-         for( icoeff = 0; icoeff < nco[ ipoly ]; icoeff++ ) {
-
-/* Store the coefficient value in the next element of the returned array,
-   if the array is not already full. Increment the pointer to the next
-   element. */
-            if( iel < nel ) coeffs[ iel++ ] = coeff[ipoly][icoeff];
-
-/* Next store the integer index of the PolyMap input or output which uses
-   the coefficient within its defining polynomial (the first axis has
-   index 1). */
-            if( iel < nel ) coeffs[ iel++ ] = ipoly + 1;
-
-/* The remaining elements of the group give the integer powers to use
-   with each input or output coordinate value. */
-            pp = power[ipoly][icoeff];
-            for( iax = 0; iax < nax; iax++,pp++ ) {
-               if( iel < nel ) coeffs[ iel++ ] = *pp;
-            }
-         }
-
-/* Increment the total number of coefficients used by the transformation. */
-         *ncoeff += nco[ ipoly ];
-      }
-   }
-}
-
-static void PolyPowers( AstPolyMap *this, double **work, int ncoord,
-                        const int *mxpow, double **ptr, int point,
-                        int fwd, int *status ){
-/*
-*+
-*  Name:
-*     astPolyPowers
-
-*  Purpose:
-*     Find the required powers of the input axis values.
-
-*  Type:
-*     Protected function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     void astPolyPowers( AstPolyMap *this, double **work, int ncoord,
-*                         const int *mxpow, double **ptr, int point,
-*                         int fwd )
-
-*  Class Membership:
-*     PolyMap virtual function.
-
-*  Description:
-*     This function is used by astTransform to calculate the powers of
-*     the axis values for a single input position. In the case of
-*     sub-classes, the powers may not be simply powers of the supplied
-*     axis values but may be more complex quantities such as a Chebyshev
-*     polynomial of the required degree evaluated at the input axis values.
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     work
-*        An array of "ncoord" pointers, each pointing to an array of
-*        length "max(2,mxpow)". The required values are placed in this
-*        array on exit.
-*     ncoord
-*        The number of axes.
-*     mxpow
-*        Pointer to an array holding the maximum power required of each
-*        axis value. Should have "ncoord" elements.
-*     ptr
-*        An array of "ncoord" pointers, each pointing to an array holding
-*        the axis values. Each of these arrays of axis values must have
-*        at least "point+1" elements.
-*     point
-*        The zero based index of the point within "ptr" that holds the
-*        axis values to be exponentiated.
-*     fwd
-*        Do the supplied coefficients define the foward transformation of
-*        the PolyMap?
-*-
-*/
-
-/* Local Variables; */
-   double *pwork;
-   double x;
-   int coord;
-   int ip;
-
-/* Check the local error status. */
-   if ( !astOK ) return;
-
-/* For the base PolyMap class, this method simply raises each input axis
-   value to the required power. Loop over all input axes. */
-   for( coord = 0; coord < ncoord; coord++ ) {
-
-/* Get a pointer to the array in which the powers of the current axis
-   value are to be returned. */
-      pwork = work[ coord ];
-
-/* Anything to the power zero is 1.0. */
-      pwork[ 0 ] = 1.0;
-
-/* Get the input axis value. If it is bad, store bad values for all
-   remaining powers. */
-      x = ptr[ coord ][ point ];
-      if( x == AST__BAD ) {
-         for( ip = 1; ip <= mxpow[ coord ]; ip++ ) pwork[ ip ] = AST__BAD;
-
-/* Otherwise, form and store the required powers of the input axis value. */
-      } else {
-         for( ip = 1; ip <= mxpow[ coord ]; ip++ ) {
-            pwork[ ip ] = pwork[ ip - 1 ]*x;
-         }
-      }
-   }
-}
-
-static AstPolyMap *PolyTran( AstPolyMap *this, int forward, double acc,
-                             double maxacc, int maxorder, const double *lbnd,
-                             const double *ubnd, int *status ){
-/*
-*++
-*  Name:
-c     astPolyTran
-f     AST_POLYTRAN
-
-*  Purpose:
-*     Fit a PolyMap inverse or forward transformation.
-
-*  Type:
-*     Public function.
-
-*  Synopsis:
-c     #include "polymap.h"
-c     AstPolyMap *astPolyTran( AstPolyMap *this, int forward, double acc,
-c                              double maxacc, int maxorder, const double *lbnd,
-c                              const double *ubnd )
-f     RESULT = AST_POLYTRAN( THIS, FORWARD, ACC, MAXACC, MAXORDER, LBND,
-f                            UBND, STATUS )
-
-*  Class Membership:
-*     PolyMap method.
-
-*  Description:
-*     This function creates a new PolyMap which is a copy of the supplied
-*     PolyMap, in which a specified transformation (forward or inverse)
-*     has been replaced by a new polynomial transformation. The
-*     coefficients of the new transformation are estimated by sampling
-*     the other transformation and performing a least squares polynomial
-*     fit in the opposite direction to the sampled positions and values.
-*
-*     This method can only be used on (1-input,1-output) or (2-input,2-output)
-*     PolyMaps.
-*
-*     The transformation to create is specified by the
-c     "forward" parameter.
-f     FORWARD parameter.
-*     In what follows "X" refers to the inputs of the PolyMap, and "Y" to
-*     the outputs of the PolyMap. The forward transformation transforms
-*     input values (X) into output values (Y), and the inverse transformation
-*     transforms output values (Y) into input values (X). Within a PolyMap,
-*     each transformation is represented by an independent set of
-*     polynomials, P_f or P_i: Y=P_f(X) for the forward transformation and
-*     X=P_i(Y) for the inverse transformation.
-*
-c     The "forward"
-f     The FORWARD
-*     parameter specifies the transformation to be replaced. If it is
-c     non-zero,
-f     is .TRUE.,
-*     a new forward transformation is created
-*     by first finding the input values (X) using the inverse transformation
-*     (which must be available) at a regular grid of points (Y) covering a
-*     rectangular region of the PolyMap's output space. The coefficients of
-*     the required forward polynomial, Y=P_f(X), are chosen in order to
-*     minimise the sum of the squared residuals between the sampled values
-*     of Y and P_f(X).
-*
-c     If "forward" is zero (probably the most likely case),
-f     If FORWARD is .FALSE. (probably the most likely case),
-*     a new inverse transformation is created by
-*     first finding the output values (Y) using the forward transformation
-*     (which must be available) at a regular grid of points (X) covering a
-*     rectangular region of the PolyMap's input space. The coefficients of
-*     the required inverse polynomial, X=P_i(Y), are chosen in order to
-*     minimise the sum of the squared residuals between the sampled values
-*     of X and P_i(Y).
-*
-*     This fitting process is performed repeatedly with increasing
-*     polynomial orders (starting with linear) until the target
-*     accuracy is achieved, or a specified maximum order is reached. If
-*     the target accuracy cannot be achieved even with this maximum-order
-*     polynomial, the best fitting maximum-order polynomial is returned so
-*     long as its accuracy is better than
-c     "maxacc".
-f     MAXACC.
-*     If it is not, a NULL pointer is returned but no error is reported.
-
-*  Parameters:
-c     this
-f     THIS = INTEGER (Given)
-*        Pointer to the original Mapping.
-c     forward
-f     FORWARD = LOGICAL (Given)
-c        If non-zero,
-f        If .TRUE.,
-*        the forward PolyMap transformation is replaced. Otherwise the
-*        inverse transformation is replaced.
-c     acc
-f     ACC = DOUBLE (Given)
-*        The target accuracy, expressed as a geodesic distance within
-*        the PolyMap's input space (if
-c        "forward" is zero)  or output space (if "forward" is non-zero).
-f        FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.).
-c     maxacc
-f     MAXACC = DOUBLE (Given)
-*        The maximum allowed accuracy for an acceptable polynomial,
-*        expressed as a geodesic distance within the PolyMap's input
-*        space (if
-c        "forward" is zero)  or output space (if "forward" is non-zero).
-f        FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.).
-c     maxorder
-f     MAXORDER = INTEGER (Given)
-*        The maximum allowed polynomial order. This is one more than the
-*        maximum power of either input axis. So for instance, a value of
-*        3 refers to a quadratic polynomial. Note, cross terms with total
-*        powers greater than or equal to
-c        maxorder
-f        MAXORDER
-*        are not inlcuded in the fit. So the maximum number of terms in
-*        each of the fitted polynomials is
-c        maxorder*(maxorder+1)/2.
-f        MAXORDER*(MAXORDER+1)/2.
-c     lbnd
-f     LBND( * ) = DOUBLE PRECISION (Given)
-c        Pointer to an
-f        An
-*        array holding the lower bounds of a rectangular region within
-*        the PolyMap's input space (if
-c        "forward" is zero)  or output space (if "forward" is non-zero).
-f        FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.).
-*        The new polynomial will be evaluated over this rectangle. The
-*        length of this array should equal the value of the PolyMap's Nin
-*        or Nout attribute, depending on
-c        "forward".
-f        FORWARD.
-c     ubnd
-f     UBND( * ) = DOUBLE PRECISION (Given)
-c        Pointer to an
-f        An
-*        array holding the upper bounds of a rectangular region within
-*        the PolyMap's input space (if
-c        "forward" is zero)  or output space (if "forward" is non-zero).
-f        FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.).
-*        The new polynomial will be evaluated over this rectangle.  The
-*        length of this array should equal the value of the PolyMap's Nin
-*        or Nout attribute, depending on
-c        "forward".
-f        FORWARD.
-f     STATUS = INTEGER (Given and Returned)
-f        The global status.
-
-*  Returned Value:
-c     astPolyTran()
-f     AST_POLYTRAN = INTEGER
-*        A pointer to the new PolyMap.
-c        A NULL pointer
-f        AST__NULL
-*        will be returned if the fit fails to achieve the accuracy
-*        specified by
-c        "maxacc",
-f        MAXACC,
-*        but no error will be reported.
-
-*  Applicability:
-*     PolyMap
-*        All PolyMaps have this method.
-*     ChebyMap
-c        The ChebyMap implementation of this method allows
-c        NULL pointers to be supplied for "lbnd" and/or "ubnd",
-c        in which case the corresponding bounds supplied when the ChebyMap
-c        was created are used.
-*        The returned PolyMap will be a ChebyMap, and the new transformation
-*        will be defined as a weighted sum of Chebyshev functions of the
-*        first kind.
-
-*  Notes:
-*     - The IterInverse attribute is always cleared in the returned PolyMap.
-*     This means that the returned PolyMap will always use the new fit by
-*     default, rather than the iterative inverse, regardless of the setting
-*     of IterInverse in the supplied PolyMap.
-*     - This function can only be used on 1D or 2D PolyMaps which have
-*     the same number of inputs and outputs.
-*     - A null Object pointer (AST__NULL) will be returned if this
-c     function is invoked with the AST error status set, or if it
-f     function is invoked with STATUS set to an error value, or if it
-*     should fail for any reason.
-*--
-*/
-
-/* Local Variables: */
-   AstPolyMap *result;
-   int ok;
-
-/* Initialise. */
-   result = NULL;
-
-/* Check the inherited status. */
-   if ( !astOK ) return result;
-
-/* Take a copy of the supplied PolyMap. */
-   result = astCopy( this );
-
-/* Replace the required transformation. */
-   ok = ReplaceTransformation( result, forward, acc, maxacc, maxorder, lbnd,
-                               ubnd, status );
-
-/* If an error occurred, or the fit was not good enough, annul the returned
-   PolyMap. */
-   if ( !ok || !astOK ) {
-      result = astAnnul( result );
-
-/* Otherwise, ensure the new fit is used in preference to the iterative
-   inverse by clearing the IterInverse attribute. */
-   } else {
-      astClearIterInverse( result );
-   }
-
-/* Return the result. */
-   return result;
-}
-
-static int ReplaceTransformation( AstPolyMap *this, int forward, double acc,
-                                  double maxacc, int maxorder, const double *lbnd,
-                                  const double *ubnd, int *status ){
-/*
-*  Name:
-*     ReplaceTransformation
-
-*  Purpose:
-*     Create a new inverse or forward transformation for a PolyMap.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     int ReplaceTransformation( AstPolyMap *this, int forward, double acc,
-*                                double maxacc, int maxorder, const double *lbnd,
-*                                const double *ubnd, int *status )
-
-*  Description:
-*     This function creates a new forward or inverse transformation for
-*     the supplied PolyMap (replacing any existing transformation), by
-*     sampling the other transformation and performing a least squares
-*     polynomial fit to the sample positions and values.
-*
-*     The transformation to create is specified by the "forward" parameter.
-*     In what follows "X" refers to the inputs of the PolyMap, and "Y" to
-*     the outputs of the PolyMap. The forward transformation transforms
-*     input values (X) into output values (Y), and the inverse transformation
-*     transforms output values (Y) into input values (X). Within a PolyMap,
-*     each transformation is represented by an independent set of
-*     polynomials: Y=P_f(X) for the forward transformation and X=P_i(Y)
-*     for the inverse transformation.
-*
-*     If "forward" is zero then a new inverse transformation is created by
-*     first finding the output values (Y) using the forward transformation
-*     (which must be available) at a regular grid of points (X) covering a
-*     rectangular region of the PolyMap's input space. The coefficients of
-*     the required inverse polynomial, X=P_i(Y), are chosen in order to
-*     minimise the sum of the squared residuals between the sampled values
-*     of X and P_i(Y).
-*
-*     If "forward" is non-zero then a new forward transformation is created
-*     by first finding the input values (X) using the inverse transformation
-*     (which must be available) at a regular grid of points (Y) covering a
-*     rectangular region of the PolyMap's output space. The coefficients of
-*     the required forward polynomial, Y=P_f(X), are chosen in order to
-*     minimise the sum of the squared residuals between the sampled values
-*     of Y and P_f(X).
-*
-*     This fitting process is performed repeatedly with increasing
-*     polynomial orders (starting with linear) until the target
-*     accuracy is achieved, or a specified maximum order is reached. If
-*     the target accuracy cannot be achieved even with this maximum-order
-*     polynomial, the best fitting maximum-order polynomial is returned so
-*     long as its accuracy is better than "maxacc".
-
-*  Parameters:
-*     this
-*        The PolyMap.
-*     forward
-*        If non-zero, then the forward PolyMap transformation is
-*        replaced. Otherwise the inverse transformation is replaced.
-*     acc
-*        The target accuracy, expressed as a geodesic distance within
-*        the PolyMap's input space (if "forward" is zero) or output
-*        space (if "forward" is non-zero).
-*     maxacc
-*        The maximum allowed accuracy for an acceptable polynomial,
-*        expressed as a geodesic distance within the PolyMap's input
-*        space (if "forward" is zero) or output space (if "forward" is
-*        non-zero).
-*     maxorder
-*        The maximum allowed polynomial order. This is one more than the
-*        maximum power of either input axis. So for instance, a value of
-*        3 refers to a quadratic polynomial. Note, cross terms with total
-*        powers greater than or equal to maxorder are not inlcuded in the
-*        fit. So the maximum number of terms in each of the fitted polynomials
-*        is maxorder*(maxorder+1)/2.
-*     lbnd
-*        An array holding the lower bounds of a rectangular region within
-*        the PolyMap's input space (if "forward" is zero) or output
-*        space (if "forward" is non-zero). The new polynomial will be
-*        evaluated over this rectangle.
-*     ubnd
-*        An array holding the upper bounds of a rectangular region within
-*        the PolyMap's input space (if "forward" is zero) or output
-*        space (if "forward" is non-zero). The new polynomial will be
-*        evaluated over this rectangle.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     - Non-zero if a fit was performed succesfully, to at least the
-*     maximum allowed accuracy. Zero if the fit failed, or an error
-*     occurred.
-
-*  Notes:
-*     - No error is reported if the fit fails to achieve the required
-*     maximum accuracy.
-*     - An error is reported if the transformation that is not being
-*     replaced is not defined.
-*     - An error is reported if the PolyMap does not have equal numbers
-*     of inputs and outputs.
-*     - An error is reported if the PolyMap has more than 2 inputs or outputs.
-
-*/
-
-/* Local Variables: */
-   double **table;
-   double *cofs;
-   double racc;
-   double scales[ 4 ];
-   int ndim;
-   int ncof;
-   int nsamp;
-   int order;
-   int result;
-
-/* Check inherited status */
-   if( !astOK ) return 0;
-
-/* Check the PolyMap can be used. */
-   ndim = astGetNin( this );
-   if( astGetNout( this ) != ndim && astOK ) {
-      astError( AST__BADNI, "astPolyTran(%s): Supplied %s has "
-                "different number of inputs (%d) and outputs (%d).",
-                status, astGetClass( this ), astGetClass( this ), ndim,
-                astGetNout( this ) );
-
-   } else if( ndim > 2 && astOK ) {
-      astError( AST__BADNI, "astPolyTran(%s): Supplied %s has "
-                "too many inputs and outputs (%d) - must be 1 or 2.",
-                status, astGetClass( this ), astGetClass( this ), ndim );
-   }
-
-   if( forward != astGetInvert( this ) ){
-      if( ! this->ncoeff_i && astOK ) {
-         astError( AST__NODEF, "astPolyTran(%s): Supplied %s has "
-                   "no inverse transformation.", status, astGetClass( this ),
-                   astGetClass( this ) );
-      }
-   } else {
-      if( ! this->ncoeff_f && astOK ) {
-         astError( AST__NODEF, "astPolyTran(%s): Supplied %s has "
-                   "no forward transformation.", status, astGetClass( this ),
-                   astGetClass( this ) );
-      }
-   }
-
-/* Check the bounds can be used. */
-   if( !lbnd || !ubnd ) {
-      astError( AST__NODEF, "astPolyTran(%s): No upper and/or lower bounds "
-                "supplied.", status, astGetClass( this ) );
-   } else if( lbnd[ 0 ] >= ubnd[ 0 ] && astOK ) {
-      astError( AST__NODEF, "astPolyTran(%s): Supplied upper "
-                "bound for the first axis (%g) is less than or equal to the "
-                "supplied lower bound (%g).", status, astGetClass( this ),
-                lbnd[ 0 ], ubnd[ 0 ] );
-   } else if( ndim == 2 && lbnd[ 1 ] >= ubnd[ 1 ] && astOK ) {
-      astError( AST__NODEF, "astPolyTran(%s): Supplied upper "
-                "bound for the second axis (%g) is less than or equal to the "
-                "supplied lower bound (%g).", status, astGetClass( this ),
-                lbnd[ 1 ], ubnd[ 1 ] );
-   }
-
-/* Initialise pointer to work space. */
-   table = NULL;
-
-/* Loop over increasing polynomial orders until the required accuracy is
-   achieved, up to the specified maximum order. The "order" value is one more
-   than the maximum power in the polynomial (so a quadratic has "order" 3). */
-   if( maxorder < 2 ) maxorder = 2;
-   for( order = 2; order <= maxorder; order++ ) {
-
-/* First do 2D PolyMaps. */
-      if( ndim == 2 ) {
-
-/* Sample the requested polynomial transformation at a grid of points. This
-   grid covers the user-supplied region, using 2*order points on each
-   axis. */
-         table = SamplePoly2D( this, !forward, table, lbnd, ubnd, 2*order,
-                               &nsamp, scales, status );
-
-/* Fit the polynomial. Always fit a linear polynomial ("order" 2) to any
-   dummy second axis. If successfull, replace the PolyMap transformation
-   and break out of the order loop. */
-         cofs = FitPoly2D( this, forward,  nsamp, acc, order, table, scales,
-                           &ncof, &racc, status );
-
-/* Now do 1D PolyMaps. */
-      } else {
-         table = SamplePoly1D( this, !forward, table, lbnd[ 0 ], ubnd[ 0 ],
-                               2*order, &nsamp, scales, status );
-         cofs = FitPoly1D( this, forward, nsamp, acc, order, table, scales,
-                           &ncof, &racc, status );
-      }
-
-/* If the fit was succesful, replace the PolyMap transformation and break
-   out of the order loop. */
-      if( cofs && ( racc < acc || ( racc < maxacc && order == maxorder ) ) ) {
-         StoreArrays( this, forward, ncof, cofs, status );
-         break;
-      } else {
-        cofs = astFree( cofs );
-      }
-   }
-
-/* If no fit was produced, return zero. */
-   result = cofs ? 1 : 0;
-
-/* Free resources. */
-   cofs = astFree( cofs );
-   table = astFreeDouble( table );
-
-/* Return the result. */
-   return result;
-}
-
-static double **SamplePoly1D( AstPolyMap *this, int forward, double **table,
-                              double lbnd, double ubnd, int npoint, int *nsamp,
-                              double scales[2], int *status ){
-/*
-*  Name:
-*     SamplePoly1D
-
-*  Purpose:
-*     Create a table of input and output positions for a 1D PolyMap.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     double **SamplePoly1D( AstPolyMap *this, int forward, double **table,
-*                            double lbnd, double ubnd, int npoint, int *nsamp,
-*                            double scales[2], int *status )
-
-*  Description:
-*     This function creates a table containing samples of the requested
-*     polynomial transformation at a grid of input points. This grid covers
-*     the user-supplied region, using "npoint" points.
-
-*  Parameters:
-*     this
-*        The PolyMap.
-*     forward
-*        If non-zero, then the forward PolyMap transformation is sampled.
-*        Otherwise the inverse transformation is sampled.
-*     table
-*        Pointer to a previous table created by this function, which is
-*        to be re-used, or NULL.
-*     lbnd
-*        The lower bounds of the region within the PolyMap's 1D input space
-*        (if "forward" is non-zero) or output space (if "forward" is zero).
-*        The new polynomial will be evaluated over this region.
-*     ubnd
-*        The upper bounds of the region within the PolyMap's 1D input space
-*        (if "forward" is non-zero) or output space (if "forward" is zero).
-*        The new polynomial will be evaluated over this region.
-*     npoint
-*        The number of points to use.
-*     nsamp
-*        Address of an int in which to return the total number of samples
-*        in the returned table.
-*     scales
-*        Array in which to return the scaling factors for the two columns
-*        of the returned table. Multiplying the returned table values by
-*        the scale factor produces PolyMap input or output axis values.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*        Pointer to an array of 2 pointers. Each of these pointers points
-*        to an array of "nsamp" doubles, being the sampled values for y1
-*        and x1 in that order. Here x1 is the input value for the sampled
-*        transformation (these are spaced on the regular grid specified
-*        by lbnd, ubnd and npoint), and y1 is the output position produced
-*        by the sampled transformation. The returned values are scaled so
-*        that each column has an RMS value of 1.0. The scaling factors that
-*        convert scaled values into original values are returned in "scales".
-*        The returned pointer should be freed using astFreeDouble when no
-*        longer needed.
-
-*/
-
-/* Local Variables: */
-   AstPointSet *ps1;
-   AstPointSet *ps2;
-   double **result;
-   double *p0;
-   double *p1;
-   double *ptr1[ 1 ];
-   double *ptr2[ 1 ];
-   double delta0;
-   double rms;
-   double sum;
-   double val0;
-   int i;
-   int icol;
-
-/* Initialise returned value */
-   result = table;
-   *nsamp = 0;
-
-/* Check inherited status */
-   if( !astOK ) return result;
-
-/* Ensure we have a table of the correct size. */
-   *nsamp = npoint;
-   if( !result ) result = astCalloc( 2, sizeof( double * ) );
-   if( result ) {
-      for( i = 0; i < 2; i++ ) {
-         result[ i ] = astRealloc( result[ i ] , (*nsamp)*sizeof( double ) );
-      }
-   }
-
-/* Work out the step sizes for the grid. */
-   delta0 = ( ubnd - lbnd )/( npoint - 1 );
-
-/* Create a PointSet to hold the grid of input positions. Use column 1
-   of the table to hold the PointSet values. */
-   ps1 = astPointSet( *nsamp, 1, " ", status );
-   ptr1[ 0 ] = result[ 1 ];
-   astSetPoints( ps1, ptr1 );
-
-/* Create a PointSet to hold the grid of output positions. Use column 0
-   of the table to hold the PointSet values. */
-   ps2 = astPointSet( *nsamp, 1, " ", status );
-   ptr2[ 0 ] = result[ 0 ];
-   astSetPoints( ps2, ptr2 );
-   if( astOK ) {
-
-/* Calculate the grid of input positions and store in the PointSet and
-   therefore also in the returned table. Rounding error may cause the
-   final point to be just over the upper bound, so we leave the final
-   point out of the loop and set it explicitly to the upper bound
-   afterwards. */
-      val0 = lbnd;
-      p0 = ptr1[ 0 ];
-      for( i = 0; i < npoint-1; i++ ) {
-         *(p0++) = val0;
-         val0 += delta0;
-      }
-      *p0 = ubnd;
-
-/* Transform the input grid to get the output grid. */
-      (void) astTransform( this, ps1, forward, ps2 );
-
-/* Scale each column in turn. */
-      for( icol = 0; icol < 2; icol++ ) {
-
-/* Find the RMS of the values in the column. */
-         sum = 0.0;
-         p0 = result[ icol ];
-         p1 = p0 + (*nsamp);
-         for( ; p0 < p1; p0++ ) sum += ( *p0 )*( *p0 );
-         rms = sqrt( sum/(*nsamp) );
-
-/* Divide the table values by the RMS. */
-         p0 = result[ icol ];
-         p1 = p0 + (*nsamp);
-         for( ; p0 < p1; p0++ ) *p0 /= rms;
-
-/* Return the RMS as the scale factor. */
-         scales[ icol ] = rms;
-      }
-   }
-
-/* Free resources */
-   ps1 = astAnnul( ps1 );
-   ps2 = astAnnul( ps2 );
-
-/* If an error occurred, free the returned array. */
-   if( !astOK ) result = astFreeDouble( result );
-
-/* Return a pointer to the table. */
-   return result;
-}
-
-static double **SamplePoly2D( AstPolyMap *this, int forward, double **table,
-                              const double *lbnd, const double *ubnd, int npoint,
-                              int *nsamp, double scales[4], int *status ){
-/*
-*  Name:
-*     SamplePoly2D
-
-*  Purpose:
-*     Create a table of input and output positions for a 2D PolyMap.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     double **SamplePoly2D( AstPolyMap *this, int forward, double **table,
-*                            const double *lbnd, const double *ubnd, int npoint,
-*                            int *nsamp, double scales[4], int *status )
-
-*  Description:
-*     This function creates a table containing samples of the requested
-*     polynomial transformation at a grid of input points. This grid covers
-*     the user-supplied region, using "npoint" points on each axis.
-
-*  Parameters:
-*     this
-*        The PolyMap.
-*     forward
-*        If non-zero, then the forward PolyMap transformation is sampled.
-*        Otherwise the inverse transformation is sampled.
-*     table
-*        Pointer to a previous table created by this function, which is
-*        to be re-used, or NULL.
-*     lbnd
-*        An array holding the lower bounds of a rectangular region within
-*        the PolyMap's input space (if "forward" is non-zero) or output
-*        space (if "forward" is zero). The new polynomial will be
-*        evaluated over this rectangle.
-*     ubnd
-*        An array holding the upper bounds of a rectangular region within
-*        the PolyMap's input space (if "forward" is non-zero) or output
-*        space (if "forward" is zero). The new polynomial will be
-*        evaluated over this rectangle.
-*     npoint
-*        The number of points along each edge of the grid.
-*     nsamp
-*        Address of an int in which to return the total number of samples
-*        in the returned table.
-*     scales
-*        Array in which to return the scaling factors for the four
-*        columns of the returned table. Multiplying the returned table
-*        values by the scale factor produces PolyMap input or output axis
-*        values.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*        Pointer to an array of 4 pointers. Each of these pointers points
-*        to an array of "nsamp" doubles, being the sampled values for y1,
-*        y2, x1 or x2 in that order. Here (x1,x2) are the input values
-*        for the sampled transformation (these are spaced on the regular
-*        grid specified by lbnd, ubnd and npoint), and (y1,y2) are the
-*        output positions produced by the sampled transformation. The
-*        returned values are scaled so that each column has an RMS value
-*        of 1.0. The scaling factors that convert scaled values into
-*        original values are returned in "scales". The returned pointer
-*        should be freed using astFreeDouble when no longer needed.
-
-*/
-
-/* Local Variables: */
-   AstPointSet *ps1;
-   AstPointSet *ps2;
-   double **result;
-   double *p0;
-   double *p1;
-   double *ptr1[ 2 ];
-   double *ptr2[ 2 ];
-   double delta1;
-   double delta0;
-   double rms;
-   double sum;
-   double val0;
-   double val1;
-   int i;
-   int icol;
-   int j;
-
-/* Initialise returned value */
-   result = table;
-   *nsamp = 0;
-
-/* Check inherited status */
-   if( !astOK ) return result;
-
-/* Ensure we have a table of the correct size. */
-   *nsamp = npoint*npoint;
-   if( !result ) result = astCalloc( 4, sizeof( double * ) );
-   if( result ) {
-      for( i = 0; i < 4; i++ ) {
-         result[ i ] = astRealloc( result[ i ] , (*nsamp)*sizeof( double ) );
-      }
-   }
-
-/* Work out the step sizes for the grid. */
-   delta0 = ( ubnd[ 0 ] - lbnd[ 0 ] )/( npoint - 1 );
-   delta1 = ( ubnd[ 1 ] - lbnd[ 1 ] )/( npoint - 1 );
-
-/* Create a PointSet to hold the grid of input positions. Use columns 2
-   and 3 of the table to hold the PointSet values. */
-   ps1 = astPointSet( *nsamp, 2, " ", status );
-   ptr1[ 0 ] = result[ 2 ];
-   ptr1[ 1 ] = result[ 3 ];
-   astSetPoints( ps1, ptr1 );
-
-/* Create a PointSet to hold the grid of output positions. Use columns 0
-   and 1 of the table to hold the PointSet values. */
-   ps2 = astPointSet( *nsamp, 2, " ", status );
-   ptr2[ 0 ] = result[ 0 ];
-   ptr2[ 1 ] = result[ 1 ];
-   astSetPoints( ps2, ptr2 );
-   if( astOK ) {
-
-/* Calculate the grid of input positions and store in the PointSet and
-   therefore also in the returned table. Rounding error may cause the
-   final point on each axis to be just over the upper bound, so we leave
-   the final point out of the loop and set it explicitly to the upper bound
-   afterwards. */
-      val0 = lbnd[ 0 ];
-      p0 = ptr1[ 0 ];
-      p1 = ptr1[ 1 ];
-      for( i = 0; i < npoint - 1; i++ ) {
-         val1 = lbnd[ 1 ];
-         for( j = 0; j < npoint-1; j++ ) {
-             *(p0++) = val0;
-             *(p1++) = val1;
-             val1 += delta1;
-         }
-         *(p0++) = val0;
-         *(p1++) = ubnd[ 1 ];
-         val0 += delta0;
-      }
-
-
-      val1 = lbnd[ 1 ];
-      for( j = 0; j < npoint-1; j++ ) {
-          *(p0++) = ubnd[ 0 ];
-          *(p1++) = val1;
-          val1 += delta1;
-      }
-      *(p0++) = ubnd[ 0 ];
-      *(p1++) = ubnd[ 1 ];
-
-
-/* Transform the input grid to get the output grid. */
-      (void) astTransform( this, ps1, forward, ps2 );
-
-/* Scale each pair of columns in turn. Use the same scale factor for
-   each axis in order to ensure an isotropic metric. */
-      for( icol = 0; icol < 4; icol += 2 ) {
-
-/* Find the RMS of the values in the two columns. */
-         sum = 0.0;
-         p0 = result[ icol ];
-         p1 = p0 + (*nsamp);
-         for( ; p0 < p1; p0++ ) sum += ( *p0 )*( *p0 );
-
-         p0 = result[ icol + 1 ];
-         p1 = p0 + (*nsamp);
-         for( ; p0 < p1; p0++ ) sum += ( *p0 )*( *p0 );
-
-         rms = sqrt( sum/(2*(*nsamp)) );
-
-/* Divide the table values by the RMS. */
-         p0 = result[ icol ];
-         p1 = p0 + (*nsamp);
-         for( ; p0 < p1; p0++ ) *p0 /= rms;
-
-         p0 = result[ icol + 1 ];
-         p1 = p0 + (*nsamp);
-         for( ; p0 < p1; p0++ ) *p0 /= rms;
-
-/* Return the RMS as the scale factor. */
-         scales[ icol ] = rms;
-         scales[ icol + 1 ] = rms;
-      }
-   }
-
-/* Free resources */
-   ps1 = astAnnul( ps1 );
-   ps2 = astAnnul( ps2 );
-
-/* If an error occurred, free the returned array. */
-   if( !astOK ) result = astFreeDouble( result );
-
-/* Return a pointer to the table. */
-   return result;
-}
-
-static void SetAttrib( AstObject *this_object, const char *setting, int *status ) {
-/*
-*  Name:
-*     SetAttrib
-
-*  Purpose:
-*     Set an attribute value for a PolyMap.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     void SetAttrib( AstObject *this, const char *setting )
-
-*  Class Membership:
-*     PolyMap member function (over-rides the astSetAttrib protected
-*     method inherited from the Mapping class).
-
-*  Description:
-*     This function assigns an attribute value for a PolyMap, the
-*     attribute and its value being specified by means of a string of
-*     the form:
-*
-*        "attribute= value "
-*
-*     Here, "attribute" specifies the attribute name and should be in
-*     lower case with no white space present. The value to the right
-*     of the "=" should be a suitable textual representation of the
-*     value to be assigned and this will be interpreted according to
-*     the attribute's data type.  White space surrounding the value is
-*     only significant for string attributes.
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     setting
-*        Pointer to a null-terminated string specifying the new attribute
-*        value.
-*/
-
-/* Local Variables: */
-   AstPolyMap *this;             /* Pointer to the PolyMap structure */
-   double dval;                  /* Floating point attribute value */
-   int ival;                     /* Integer attribute value */
-   int len;                      /* Length of setting string */
-   int nc;                       /* Number of characters read by astSscanf */
-
-/* Check the global error status. */
-   if ( !astOK ) return;
-
-/* Obtain a pointer to the PolyMap structure. */
-   this = (AstPolyMap *) this_object;
-
-/* Obtain the length of the setting string. */
-   len = (int) strlen( setting );
-
-/* Test for each recognised attribute in turn, using "astSscanf" to parse
-   the setting string and extract the attribute value (or an offset to
-   it in the case of string values). In each case, use the value set
-   in "nc" to check that the entire string was matched. Once a value
-   has been obtained, use the appropriate method to set it. */
-
-/* IterInverse. */
-/* ------------ */
-   if ( nc = 0,
-        ( 1 == astSscanf( setting, "iterinverse= %d %n", &ival, &nc ) )
-        && ( nc >= len ) ) {
-      astSetIterInverse( this, ival );
-
-/* NiterInverse. */
-/* ------------- */
-   } else if ( nc = 0,
-        ( 1 == astSscanf( setting, "niterinverse= %d %n", &ival, &nc ) )
-        && ( nc >= len ) ) {
-      astSetNiterInverse( this, ival );
-
-/* TolInverse. */
-/* ----------- */
-   } else if ( nc = 0,
-        ( 1 == astSscanf( setting, "tolinverse= %lg %n", &dval, &nc ) )
-        && ( nc >= len ) ) {
-      astSetTolInverse( this, dval );
-
-/* If the attribute is still not recognised, pass it on to the parent
-   method for further interpretation. */
-   } else {
-      (*parent_setattrib)( this_object, setting, status );
-   }
-}
-
-static void StoreArrays( AstPolyMap *this, int forward, int ncoeff,
-                         const double *coeff, int *status ){
-/*
-*  Name:
-*     StoreArrays
-
-*  Purpose:
-*     Store the dynamic arrays for a single transformation within a PolyMap
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     void StoreArrays( AstPolyMap *this, int forward, int ncoeff,
-*                       const double *coeff, int *status )
-
-*  Class Membership:
-*     PolyMap initialiser.
-
-*  Description:
-*     This function sets up the arrays within a PolyMap structure that
-*     describes either the forward or inverse transformation.
-
-*  Parameters:
-*     this
-*        The PolyMap.
-*     forward
-*        If non-zero, replace the forward transformation. Otherwise,
-*        replace the inverse transformation.
-*     ncoeff
-*        The number of non-zero coefficients necessary to define the
-*        specified transformation of the PolyMap. If zero is supplied, the
-*        transformation will be undefined.
-*     coeff
-*        An array containing "ncof*( 2 + nin )" elements. Each group
-*	 of "2 + nin" adjacent elements describe a single coefficient of
-*	 the transformation. Within each such group, the first
-*	 element is the coefficient value; the next element is the
-*	 integer index of the PolyMap output which uses the coefficient
-*	 within its defining polynomial (the first output has index 1);
-*	 the remaining elements of the group give the integer powers to
-*	 use with each input coordinate value (powers must not be
-*	 negative).
-*     status
-*        Pointer to inherited status.
-*/
-
-/* Local Variables: */
-   const double *group;          /* Pointer to start of next coeff. description */
-   int *pows;                    /* Pointer to powers for current coeff. */
-   int gsize;                    /* Length of each coeff. description */
-   int i;                        /* Loop count */
-   int ico;                      /* Index of next coeff. for current input or output */
-   int iin;                      /* Input index extracted from coeff. description */
-   int iout;                     /* Output index extracted from coeff. description */
-   int j;                        /* Loop count */
-   int nin;                      /* Number of orignal inputs */
-   int nout;                     /* Number of original outputs */
-   int pow;                      /* Power extracted from coeff. description */
-
-/* Check the global status. */
-   if ( !astOK ) return;
-
-/* First Free any existing arrays. */
-   FreeArrays( this, forward, status );
-
-/* Get the number of inputs and outputs of the original uninverted PolyMap.
-   Swap the transformation to be set if the PolyMap has been inverted. */
-   if( astGetInvert( this ) ) {
-      nout = astGetNin( this );
-      nin = astGetNout( this );
-      forward = !forward;
-   } else {
-      nin = astGetNin( this );
-      nout = astGetNout( this );
-   }
-
-/* Now initialise the original forward transformation, if required. */
-   if( forward && ncoeff > 0 ) {
-
-/* Create the arrays decribing the forward transformation. */
-      this->ncoeff_f = astMalloc( sizeof( int )*(size_t) nout );
-      this->mxpow_f = astMalloc( sizeof( int )*(size_t) nin );
-      this->power_f = astMalloc( sizeof( int ** )*(size_t) nout );
-      this->coeff_f = astMalloc( sizeof( double * )*(size_t) nout );
-      if( astOK ) {
-
-/* Initialise the count of coefficients for each output coordinate to zero. */
-         for( i = 0; i < nout; i++ ) this->ncoeff_f[ i ] = 0;
-
-/* Initialise max power for each input coordinate to zero. */
-         for( j = 0; j < nin; j++ ) this->mxpow_f[ j ] = 0;
-
-/* Scan through the supplied forward coefficient array, counting the
-   number of coefficients which relate to each output. Also find the
-   highest power used for each input axis. Report errors if any unusable
-   values are found in the supplied array. */
-         group = coeff;
-         gsize = 2 + nin;
-         for( i = 0; i < ncoeff && astOK; i++, group += gsize ) {
-
-            iout = floor( group[ 1 ] + 0.5 );
-            if( iout < 1 || iout > nout ) {
-               astError( AST__BADCI, "astInitPolyMap(%s): Forward "
-                         "coefficient %d referred to an illegal output "
-                         "coordinate %d.", status, astGetClass( this ), i + 1,
-                         iout );
-               astError( AST__BADCI, "This number should be in the "
-                         "range 1 to %d.", status, nout );
-               break;
-            }
-
-            this->ncoeff_f[ iout - 1 ]++;
-
-            for( j = 0; j < nin; j++ ) {
-               pow = floor( group[ 2 + j ] + 0.5 );
-               if( pow < 0 ) {
-                  astError( AST__BADPW, "astInitPolyMap(%s): Forward "
-                            "coefficient %d has a negative power (%d) "
-                            "for input coordinate %d.", status,
-                            astGetClass( this ), i + 1, pow, j + 1 );
-                  astError( AST__BADPW, "All powers should be zero or "
-                            "positive." , status);
-                  break;
-               }
-               if( pow > this->mxpow_f[ j ] ) this->mxpow_f[ j ] = pow;
-            }
-         }
-
-/* Allocate the arrays to store the input powers associated with each
-   coefficient, and the coefficient values. Reset the coefficient count
-   for each axis to zero afterwards so that we can use the array as an index
-   to the next vacant slot within the following loop. */
-         for( i = 0; i < nout; i++ ) {
-            this->power_f[ i ] = astMalloc( sizeof( int * )*
-                                           (size_t) this->ncoeff_f[ i ] );
-            this->coeff_f[ i ] = astMalloc( sizeof( double )*
-                                           (size_t) this->ncoeff_f[ i ] );
-            this->ncoeff_f[ i ] = 0;
-         }
-
-         if( astOK ) {
-
-/* Extract the coefficient values and powers form the supplied array and
-   store them in the arrays created above. */
-            group = coeff;
-            for( i = 0; i < ncoeff && astOK; i++, group += gsize ) {
-               iout = floor( group[ 1 ] + 0.5 ) - 1;
-               ico = ( this->ncoeff_f[ iout ] )++;
-               this->coeff_f[ iout ][ ico ] = group[ 0 ];
-
-               pows = astMalloc( sizeof( int )*(size_t) nin );
-               this->power_f[ iout ][ ico ] = pows;
-               if( astOK ) {
-                  for( j = 0; j < nin; j++ ) {
-                     pows[ j ] = floor( group[ 2 + j ] + 0.5 );
-                  }
-               }
-            }
-         }
-      }
-   }
-
-/* Now initialise the original inverse transformation, if required. */
-   if( !forward && ncoeff > 0 ) {
-
-/* Create the arrays decribing the inverse transformation. */
-      this->ncoeff_i = astMalloc( sizeof( int )*(size_t) nin );
-      this->mxpow_i = astMalloc( sizeof( int )*(size_t) nout );
-      this->power_i = astMalloc( sizeof( int ** )*(size_t) nin );
-      this->coeff_i = astMalloc( sizeof( double * )*(size_t) nin );
-      if( astOK ) {
-
-/* Initialise the count of coefficients for each input coordinate to zero. */
-         for( i = 0; i < nin; i++ ) this->ncoeff_i[ i ] = 0;
-
-/* Initialise max power for each output coordinate to zero. */
-         for( j = 0; j < nout; j++ ) this->mxpow_i[ j ] = 0;
-
-/* Scan through the supplied inverse coefficient array, counting the
-   number of coefficients which relate to each input. Also find the
-   highest power used for each output axis. Report errors if any unusable
-   values are found in the supplied array. */
-         group = coeff;
-
-         gsize = 2 + nout;
-         for( i = 0; i < ncoeff && astOK; i++, group += gsize ) {
-
-            iin = floor( group[ 1 ] + 0.5 );
-            if( iin < 1 || iin > nin ) {
-               astError( AST__BADCI, "astInitPolyMap(%s): Inverse "
-                         "coefficient %d referred to an illegal input "
-                         "coordinate %d.", status, astGetClass( this ),
-                         i + 1, iin );
-               astError( AST__BADCI, "This number should be in the "
-                         "range 1 to %d.", status, nin );
-               break;
-            }
-
-            this->ncoeff_i[ iin - 1 ]++;
-
-            for( j = 0; j < nout; j++ ) {
-               pow = floor( group[ 2 + j ] + 0.5 );
-               if( pow < 0 ) {
-                  astError( AST__BADPW, "astInitPolyMap(%s): Inverse "
-                            "coefficient %d has a negative power (%d) "
-                            "for output coordinate %d.", status,
-                            astGetClass( this ), i + 1, pow, j + 1 );
-                  astError( AST__BADPW, "All powers should be zero or "
-                            "positive." , status);
-                  break;
-               }
-               if( pow > this->mxpow_i[ j ] ) this->mxpow_i[ j ] = pow;
-            }
-         }
-
-/* Allocate the arrays to store the output powers associated with each
-   coefficient, and the coefficient values. Reset the coefficient count
-   for each axis to zero afterwards so that we can use the array as an index
-   to the next vacant slot within the following loop. */
-         for( i = 0; i < nin; i++ ) {
-            this->power_i[ i ] = astMalloc( sizeof( int * )*
-                                           (size_t) this->ncoeff_i[ i ] );
-            this->coeff_i[ i ] = astMalloc( sizeof( double )*
-                                           (size_t) this->ncoeff_i[ i ] );
-            this->ncoeff_i[ i ] = 0;
-         }
-
-         if( astOK ) {
-
-/* Extract the coefficient values and powers form the supplied array and
-   store them in the arrays created above. */
-            group = coeff;
-            for( i = 0; i < ncoeff && astOK; i++, group += gsize ) {
-               iin = floor( group[ 1 ] + 0.5 ) - 1;
-               ico = ( this->ncoeff_i[ iin ] )++;
-               this->coeff_i[ iin ][ ico ] = group[ 0 ];
-
-               pows = astMalloc( sizeof( int )*(size_t) nout );
-               this->power_i[ iin ][ ico ] = pows;
-               if( astOK ) {
-                  for( j = 0; j < nout; j++ ) {
-                     pows[ j ] = floor( group[ 2 + j ] + 0.5 );
-                  }
-               }
-            }
-         }
-      }
-   }
-}
-
-static int TestAttrib( AstObject *this_object, const char *attrib, int *status ) {
-/*
-*  Name:
-*     TestAttrib
-
-*  Purpose:
-*     Test if a specified attribute value is set for a PolyMap.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     int TestAttrib( AstObject *this, const char *attrib, int *status )
-
-*  Class Membership:
-*     PolyMap member function (over-rides the astTestAttrib protected
-*     method inherited from the Mapping class).
-
-*  Description:
-*     This function returns a boolean result (0 or 1) to indicate whether
-*     a value has been set for one of a PolyMap's attributes.
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     attrib
-*        Pointer to a null-terminated string specifying the attribute
-*        name.  This should be in lower case with no surrounding white
-*        space.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     One if a value has been set, otherwise zero.
-
-*  Notes:
-*     - A value of zero will be returned if this function is invoked
-*     with the global status set, or if it should fail for any reason.
-*/
-
-/* Local Variables: */
-   AstPolyMap *this;             /* Pointer to the PolyMap structure */
-   int result;                   /* Result value to return */
-
-/* Initialise. */
-   result = 0;
-
-/* Check the global error status. */
-   if ( !astOK ) return result;
-
-/* Obtain a pointer to the PolyMap structure. */
-   this = (AstPolyMap *) this_object;
-
-/* Check the attribute name and test the appropriate attribute. */
-
-/* IterInverse. */
-/* ------------ */
-   if ( !strcmp( attrib, "iterinverse" ) ) {
-      result = astTestIterInverse( this );
-
-/* NiterInverse. */
-/* ------------- */
-   } else if ( !strcmp( attrib, "niterinverse" ) ) {
-      result = astTestNiterInverse( this );
-
-/* TolInverse. */
-/* ----------- */
-   } else if ( !strcmp( attrib, "tolinverse" ) ) {
-      result = astTestTolInverse( this );
-
-/* If the attribute is still not recognised, pass it on to the parent
-   method for further interpretation. */
-   } else {
-      result = (*parent_testattrib)( this_object, attrib, status );
-   }
-
-/* Return the result, */
-   return result;
-}
-
-static AstPointSet *Transform( AstMapping *this, AstPointSet *in,
-                               int forward, AstPointSet *out, int *status ) {
-/*
-*  Name:
-*     Transform
-
-*  Purpose:
-*     Apply a PolyMap to transform a set of points.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     AstPointSet *Transform( AstMapping *this, AstPointSet *in,
-*                             int forward, AstPointSet *out, int *status )
-
-*  Class Membership:
-*     PolyMap member function (over-rides the astTransform protected
-*     method inherited from the Mapping class).
-
-*  Description:
-*     This function takes a PolyMap and a set of points encapsulated in a
-*     PointSet and transforms the points.
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap.
-*     in
-*        Pointer to the PointSet holding the input coordinate data.
-*     forward
-*        A non-zero value indicates that the forward coordinate transformation
-*        should be applied, while a zero value requests the inverse
-*        transformation.
-*     out
-*        Pointer to a PointSet which will hold the transformed (output)
-*        coordinate values. A NULL value may also be given, in which case a
-*        new PointSet will be created by this function.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     Pointer to the output (possibly new) PointSet.
-
-*  Notes:
-*     -  A null pointer will be returned if this function is invoked with the
-*     global error status set, or if it should fail for any reason.
-*     -  The number of coordinate values per point in the input PointSet must
-*     match the number of columns in the PolyMap being applied.
-*     -  The number of coordinate values per point in the output PointSet will
-*     equal the number of rows in the PolyMap being applied.
-*     -  If an output PointSet is supplied, it must have space for sufficient
-*     number of points and coordinate values per point to accommodate the
-*     result. Any excess space will be ignored.
-*/
-
-/* Local Variables: */
-   AstPointSet *result;          /* Pointer to output PointSet */
-   AstPolyMap *map;              /* Pointer to PolyMap to be applied */
-   double **coeff;               /* Pointer to coefficient value arrays */
-   double **ptr_in;              /* Pointer to input coordinate data */
-   double **ptr_out;             /* Pointer to output coordinate data */
-   double **work;                /* Pointer to exponentiated axis values */
-   double *outcof;               /* Pointer to next coefficient value */
-   double outval;                /* Output axis value */
-   double term;                  /* Term to be added to output value */
-   double xp;                    /* Exponentiated input axis value */
-   int ***power;                 /* Pointer to coefficient power arrays */
-   int **outpow;                 /* Pointer to next set of axis powers */
-   int *mxpow;                   /* Pointer to max used power for each input */
-   int *ncoeff;                  /* Pointer to no. of coefficients */
-   int in_coord;                 /* Index of output coordinate */
-   int ico;                      /* Coefficient index */
-   int nc;                       /* No. of coefficients in polynomial */
-   int ncoord_in;                /* Number of coordinates per input point */
-   int ncoord_out;               /* Number of coordinates per output point */
-   int npoint;                   /* Number of points */
-   int out_coord;                /* Index of output coordinate */
-   int point;                    /* Loop counter for points */
-   int pow;                      /* Next axis power */
-
-/* Check the global error status. */
-   if ( !astOK ) return NULL;
-
-/* Obtain a pointer to the PolyMap. */
-   map = (AstPolyMap *) this;
-
-/* Apply the parent mapping using the stored pointer to the Transform member
-   function inherited from the parent Mapping class. This function validates
-   all arguments and generates an output PointSet if necessary, but does not
-   actually transform any coordinate values. */
-   result = (*parent_transform)( this, in, forward, out, status );
-
-/* Determine whether to apply the original forward or inverse mapping,
-   according to the direction specified and whether the mapping has been
-   inverted. */
-   if ( astGetInvert( map ) ) forward = !forward;
-
-/* We will now extend the parent astTransform method by performing the
-   calculations needed to generate the output coordinate values. */
-
-/* If we are using the original inverse transformatiom, and the IterInverse
-   attribute is non-zero, use an iterative inverse algorithm rather than any
-   inverse transformation defined within the PolyMap. */
-   if( !forward && astGetIterInverse(map) ) {
-      IterInverse( map, in, result, status );
-
-/* Otherwise, determine the numbers of points and coordinates per point from
-   the input and output PointSets and obtain pointers for accessing the input
-   and output coordinate values. */
-   } else {
-      ncoord_in = astGetNcoord( in );
-      ncoord_out = astGetNcoord( result );
-      npoint = astGetNpoint( in );
-      ptr_in = astGetPoints( in );
-      ptr_out = astGetPoints( result );
-
-/* Get a pointer to the arrays holding the required coefficient
-   values and powers, according to the direction of mapping required. */
-      if ( forward ) {
-         ncoeff = map->ncoeff_f;
-         coeff = map->coeff_f;
-         power = map->power_f;
-         mxpow = map->mxpow_f;
-      } else {
-         ncoeff = map->ncoeff_i;
-         coeff = map->coeff_i;
-         power = map->power_i;
-         mxpow = map->mxpow_i;
-      }
-
-/* Allocate memory to hold the required powers of the input axis values. */
-      work = astMalloc( sizeof( double * )*(size_t) ncoord_in );
-      for( in_coord = 0; in_coord < ncoord_in; in_coord++ ) {
-         work[ in_coord ] = astMalloc( sizeof( double )*
-                           (size_t) ( astMAX( 2, mxpow[in_coord]+1 ) ) );
-      }
-
-/* Perform coordinate arithmetic. */
-/* ------------------------------ */
-      if ( astOK ) {
-
-/* Loop to apply the polynomial to each point in turn.*/
-         for ( point = 0; point < npoint; point++ ) {
-
-/* Find the required powers of the input axis values and store them
-   in the work array. Note, using a virtual method here slows the PolyMap
-   Transform function down by about 5%, compared to doing the equivalent
-   calculations in-line. But we need some way to allow the ChebyMap class
-   to over-ride the calculation of the powers, so we must do something
-   like this. If the 5% slow-down is too much, it can be reduced down to
-   about 2% by replacing the invocation of the astPolyPowers_ interface
-   function with a direct call to the implementation function itself.
-   This involves replacing the astPolyPowers call below with this:
-
-   (**astMEMBER(this,PolyMap,PolyPowers))( (AstPolyMap *) this, work, ncoord_in,
-                                           mxpow, ptr_in, point, forward, status );
-
-   The above could be wrapped up in an alternative implementation of the
-   astPolyPowers macro, so that it looks the same as the existing code.
-   In fact, this scheme could be more widely used to speed up invocation
-   of virtual functions within AST. The disadvantage is that the interface
-   functions for some virtual methods includes some extra processing,
-   over and above simply invoking the implementation function.
-
-   Of course the other way to get rid of the 5% slow down, is to
-   revert to using in-line code below, and then replicate this entire
-   function in the ChebyMap class, making suitable changes to use
-   Chebyshev functions in place of simple powers. But that is bad
-   structuring... */
-            astPolyPowers( this, work, ncoord_in, mxpow, ptr_in, point,
-                           forward );
-
-/* Loop round each output. */
-            for( out_coord = 0; out_coord < ncoord_out; out_coord++ ) {
-
-/* Initialise the output value. */
-               outval = 0.0;
-
-/* Get pointers to the coefficients and powers for this output. */
-               outcof = coeff[ out_coord ];
-               outpow = power[ out_coord ];
-
-/* Loop round all polynomial coefficients.*/
-               nc = ncoeff[ out_coord ];
-               for ( ico = 0; ico < nc && outval != AST__BAD;
-                     ico++, outcof++, outpow++ ) {
-
-/* Initialise the current term to be equal to the value of the coefficient.
-   If it is bad, store a bad output value. */
-                  term = *outcof;
-                  if( term == AST__BAD ) {
-                     outval = AST__BAD;
-
-/* Otherwise, loop round all inputs */
-                  } else {
-                     for( in_coord = 0; in_coord < ncoord_in; in_coord++ ) {
-
-/* Get the power of the current input axis value used by the current
-   coefficient. If it is zero, pass on. */
-                        pow = (*outpow)[ in_coord ];
-                        if( pow > 0 ) {
-
-/* Get the axis value raised to the appropriate power. */
-                           xp = work[ in_coord ][ pow ];
-
-/* If bad, set the output value bad and break. */
-                           if( xp == AST__BAD ) {
-                              outval = AST__BAD;
-                              break;
-
-/* Otherwise multiply the current term by the exponentiated axis value. */
-                           } else {
-                              term *= xp;
-                           }
-                        }
-                     }
-                  }
-
-/* Increment the output value by the current term of the polynomial. */
-                  if( outval != AST__BAD ) outval += term;
-
-               }
-
-/* Store the output value. */
-               ptr_out[ out_coord ][ point ] = outval;
-
-            }
-         }
-      }
-
-/* Free resources. */
-      for( in_coord = 0; in_coord < ncoord_in; in_coord++ ) {
-         work[ in_coord ] = astFree( work[ in_coord ] );
-      }
-      work = astFree( work );
-   }
-
-/* Return a pointer to the output PointSet. */
-   return result;
-}
-
-/* Functions which access class attributes. */
-/* ---------------------------------------- */
-/* Implement member functions to access the attributes associated with
-   this class using the macros defined for this purpose in the
-   "object.h" file. For a description of each attribute, see the class
-   interface (in the associated .h file). */
-
-/*
-*att++
-*  Name:
-*     IterInverse
-
-*  Purpose:
-*     Provide an iterative inverse transformation?
-
-*  Type:
-*     Public attribute.
-
-*  Synopsis:
-*     Integer (boolean).
-
-*  Description:
-*     This attribute indicates whether the original inverse transformation
-*     of the PolyMap should be implemented via an iterative Newton-Raphson
-*     approximation that uses the forward transformation to transform
-*     candidate input positions until an output position is found which
-*     is close to the required output position. By default, an iterative
-*     inverse is provided if, and only if, no inverse polynomial was supplied
-*     when the PolyMap was constructed.
-*
-*     Note, the term "inverse transformation" here refers to the inverse
-*     transformation of the original PolyMap, ignoring any subsequent
-*     inversions. Also, "input" and "output" refer to the inputs and
-*     outputs of the original PolyMap.
-
-*     The NiterInverse and TolInverse attributes provide parameters that
-*     control the behaviour of the inverse approximation method.
-*
-*     The iterative inverse returns AST__BAD axis values at positions
-*     for which the inverse transformation is undefined. For instance,
-*     if the forward transformation is y = x*x, the iterative inverse
-*     will return x = AST__BAD at y = -1. If the inverse transformation
-*     is multiply defined, the position returned by the iterative inverse
-*     will be the position of the solution that is closest to the
-*     supplied position. For instance, using the above example, y = x*x,
-*     the iterative inverse will return x = +2 at y = 4, because x = +2
-*     is the closest solution to 4 (the other solution is x = -2).
-
-*  Applicability:
-*     PolyMap
-*        All PolyMaps have this attribute.
-*     ChebyMap
-*        The ChebyMap class does not currently provide an option for an
-*        iterative inverse, and so the IterInverse value is always zero.
-*        Setting or clearing the IterInverse attribute of a ChebyMap has
-*        no effect.
-
-*  Notes:
-*     - The transformation replaced by the iterative algorithm is the
-*     transformation from the original PolyMap output space to the
-*     original PolyMap input space (i.e. the input and output spaces as
-*     defined by the arguments of the PolyMap constructor). This is still
-*     the case even if the PolyMap has subsequently been inverted. In
-*     other words if a PolyMap is created and then inverted, setting
-*     the IterInverse to a non-zero value will replace the forward
-*     transformation of the inverted PolyMap (i.e. the inverse
-*     transformation of the original PolyMap). It is not possible to
-*     replace the other transformation (i.e. from the original PolyMap
-*     input space to the original PolyMap output space) with an iterative
-*     algorithm.
-*     - If a PolyMap that has an iterative inverse transformation is
-*     subsequently inverted, the inverted PolyMap will have an iterative
-*     forward transformation.
-*     - An iterative inverse can only be used if the PolyMap has equal
-*     numbers of inputs and outputs, as given by the Nin and Nout
-*     attributes. An error will be reported if IterInverse is set non-zero
-*     for a PolyMap that does not meet this requirement.
-
-*att--
-*/
-astMAKE_CLEAR(PolyMap,IterInverse,iterinverse,-INT_MAX)
-astMAKE_GET(PolyMap,IterInverse,int,0,( ( this->iterinverse == -INT_MAX ) ?
-                                          (this->ncoeff_i == 0) : this->iterinverse ))
-astMAKE_SET(PolyMap,IterInverse,int,iterinverse,
-  (((astGetNin(this)==astGetNout(this))||!value)?((value?1:0)):(astError(AST__ATTIN,"astSetIterInverse(%s):"
-  "Cannot use an iterative inverse because the %s has unequal numbers of "
-  "inputs and outputs.", status, astGetClass(this),astGetClass(this)),this->iterinverse)))
-astMAKE_TEST(PolyMap,IterInverse,( this->iterinverse != -INT_MAX ))
-
-/* NiterInverse. */
-/* --------- */
-/*
-*att++
-*  Name:
-*     NiterInverse
-
-*  Purpose:
-*     Maximum number of iterations for the iterative inverse transformation.
-
-*  Type:
-*     Public attribute.
-
-*  Synopsis:
-*     Integer.
-
-*  Description:
-*     This attribute controls the iterative inverse transformation
-*     used if the IterInverse attribute is non-zero.
-*
-*     Its value gives the maximum number of iterations of the
-*     Newton-Raphson algorithm to be used for each transformed position.
-*     The default value is 4. See also attribute TolInverse.
-
-*  Applicability:
-*     PolyMap
-*        All PolyMaps have this attribute.
-
-*att--
-*/
-astMAKE_CLEAR(PolyMap,NiterInverse,niterinverse,-INT_MAX)
-astMAKE_GET(PolyMap,NiterInverse,int,0,( this->niterinverse == -INT_MAX ? 4 : this->niterinverse))
-astMAKE_SET(PolyMap,NiterInverse,int,niterinverse,value)
-astMAKE_TEST(PolyMap,NiterInverse,( this->niterinverse != -INT_MAX ))
-
-/* TolInverse. */
-/* ----------- */
-/*
-*att++
-*  Name:
-*     TolInverse
-
-*  Purpose:
-*     Target relative error for the iterative inverse transformation.
-
-*  Type:
-*     Public attribute.
-
-*  Synopsis:
-*     Floating point.
-
-*  Description:
-*     This attribute controls the iterative inverse transformation
-*     used if the IterInverse attribute is non-zero.
-*
-*     Its value gives the target relative error in the axis values of
-*     each transformed position. Further iterations will be performed
-*     until the target relative error is reached, or the maximum number
-*     of iterations given by attribute NiterInverse is reached.
-
-*     The default value is 1.0E-6.
-
-*  Applicability:
-*     PolyMap
-*        All PolyMaps have this attribute.
-*att--
-*/
-astMAKE_CLEAR(PolyMap,TolInverse,tolinverse,AST__BAD)
-astMAKE_GET(PolyMap,TolInverse,double,0.0,( this->tolinverse == AST__BAD ? 1.0E-6 : this->tolinverse))
-astMAKE_SET(PolyMap,TolInverse,double,tolinverse,value)
-astMAKE_TEST(PolyMap,TolInverse,( this->tolinverse != AST__BAD ))
-
-/* Copy constructor. */
-/* ----------------- */
-static void Copy( const AstObject *objin, AstObject *objout, int *status ) {
-/*
-*  Name:
-*     Copy
-
-*  Purpose:
-*     Copy constructor for PolyMap objects.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     void Copy( const AstObject *objin, AstObject *objout, int *status )
-
-*  Description:
-*     This function implements the copy constructor for PolyMap objects.
-
-*  Parameters:
-*     objin
-*        Pointer to the object to be copied.
-*     objout
-*        Pointer to the object being constructed.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     void
-
-*  Notes:
-*     -  This constructor makes a deep copy, including a copy of the
-*     coefficients associated with the input PolyMap.
-*/
-
-
-/* Local Variables: */
-   AstPolyMap *in;               /* Pointer to input PolyMap */
-   AstPolyMap *out;              /* Pointer to output PolyMap */
-   int nin;                      /* No. of input coordinates */
-   int nout;                     /* No. of output coordinates */
-   int i;                        /* Loop count */
-   int j;                        /* Loop count */
-
-/* Check the global error status. */
-   if ( !astOK ) return;
-
-/* Obtain pointers to the input and output PolyMaps. */
-   in = (AstPolyMap *) objin;
-   out = (AstPolyMap *) objout;
-
-/* Nullify the pointers stored in the output object since these will
-   currently be pointing at the input data (since the output is a simple
-   byte-for-byte copy of the input). Otherwise, the input data could be
-   freed by accidient if the output object is deleted due to an error
-   occuring in this function. */
-   out->ncoeff_f = NULL;
-   out->power_f = NULL;
-   out->coeff_f = NULL;
-   out->mxpow_f = NULL;
-
-   out->ncoeff_i = NULL;
-   out->power_i = NULL;
-   out->coeff_i = NULL;
-   out->mxpow_i = NULL;
-
-   out->jacobian = NULL;
-   out->lintrunc = NULL;
-
-/* Get the number of inputs and outputs of the uninverted Mapping. */
-   nin = ( (AstMapping *) in )->nin;
-   nout = ( (AstMapping *) in )->nout;
-
-/* Copy the number of coefficients associated with each output of the
-   forward transformation of the uninverted Mapping. */
-   if( in->ncoeff_f ) {
-      out->ncoeff_f = (int *) astStore( NULL, (void *) in->ncoeff_f,
-                                        sizeof( int )*(size_t) nout );
-
-/* Copy the maximum power of each input axis value used by the forward
-   transformation. */
-      out->mxpow_f = (int *) astStore( NULL, (void *) in->mxpow_f,
-                                       sizeof( int )*(size_t) nin );
-
-/* Copy the coefficient values used by the forward transformation. */
-      if( in->coeff_f ) {
-         out->coeff_f = astMalloc( sizeof( double * )*(size_t) nout );
-         if( astOK ) {
-            for( i = 0; i < nout; i++ ) {
-               out->coeff_f[ i ] = (double *) astStore( NULL, (void *) in->coeff_f[ i ],
-                                                     sizeof( double )*(size_t) in->ncoeff_f[ i ] );
-            }
-         }
-      }
-
-/* Copy the input axis powers associated with each coefficient of the forward
-   transformation. */
-      if( in->power_f ) {
-         out->power_f = astMalloc( sizeof( int ** )*(size_t) nout );
-         if( astOK ) {
-            for( i = 0; i < nout; i++ ) {
-               out->power_f[ i ] = astMalloc( sizeof( int * )*(size_t) in->ncoeff_f[ i ] );
-               if( astOK ) {
-                  for( j = 0; j < in->ncoeff_f[ i ]; j++ ) {
-                     out->power_f[ i ][ j ] = (int *) astStore( NULL, (void *) in->power_f[ i ][ j ],
-                                                                sizeof( int )*(size_t) nin );
-                  }
-               }
-            }
-         }
-      }
-   }
-
-/* Do the same for the inverse transformation. */
-   if( in->ncoeff_i ) {
-      out->ncoeff_i = (int *) astStore( NULL, (void *) in->ncoeff_i,
-                                        sizeof( int )*(size_t) nin );
-
-      out->mxpow_i = (int *) astStore( NULL, (void *) in->mxpow_i,
-                                       sizeof( int )*(size_t) nout );
-
-      if( in->coeff_i ) {
-         out->coeff_i = astMalloc( sizeof( double * )*(size_t) nin );
-         if( astOK ) {
-            for( i = 0; i < nin; i++ ) {
-               out->coeff_i[ i ] = (double *) astStore( NULL, (void *) in->coeff_i[ i ],
-                                                     sizeof( double )*(size_t) in->ncoeff_i[ i ] );
-            }
-         }
-      }
-
-      if( in->power_i ) {
-         out->power_i = astMalloc( sizeof( int ** )*(size_t) nin );
-         if( astOK ) {
-            for( i = 0; i < nin; i++ ) {
-               out->power_i[ i ] = astMalloc( sizeof( int * )*(size_t) in->ncoeff_i[ i ] );
-               if( astOK ) {
-                  for( j = 0; j < in->ncoeff_i[ i ]; j++ ) {
-                     out->power_i[ i ][ j ] = (int *) astStore( NULL, (void *) in->power_i[ i ][ j ],
-                                                                sizeof( int )*(size_t) nout );
-                  }
-               }
-            }
-         }
-      }
-   }
-
-/* Copy the linear truncation of the PolyMap - if it has been found. */
-   if( in->lintrunc ) out->lintrunc = astCopy( in->lintrunc );
-
-/* If an error has occurred, free all the resources allocated above. */
-   if( !astOK ) {
-      FreeArrays( out, 1, status );
-      FreeArrays( out, 0, status );
-   }
-
-   return;
-
-}
-
-/* Destructor. */
-/* ----------- */
-static void Delete( AstObject *obj, int *status ) {
-/*
-*  Name:
-*     Delete
-
-*  Purpose:
-*     Destructor for PolyMap objects.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     void Delete( AstObject *obj, int *status )
-
-*  Description:
-*     This function implements the destructor for PolyMap objects.
-
-*  Parameters:
-*     obj
-*        Pointer to the object to be deleted.
-*     status
-*        Pointer to the inherited status variable.
-
-*  Returned Value:
-*     void
-
-*  Notes:
-*     This function attempts to execute even if the global error status is
-*     set.
-*/
-
-/* Local Variables: */
-   AstPolyMap *this;
-   int nc;
-   int ic;
-   int lstat;
-   int error;
-
-/* Obtain a pointer to the PolyMap structure. */
-   this = (AstPolyMap *) obj;
-
-/* Free the arrays. */
-   FreeArrays( this, 1, status );
-   FreeArrays( this, 0, status );
-
-/* Free the resources used to store the Jacobian of the forward
-   transformation. */
-   if( this->jacobian ) {
-
-/* Get the number of PolyMap inputs. We need to clear any error status
-   first since astGetNin returns zero if an error has occurred. The
-   Jacobian will only be non-NULL if the number of inputs and outputs
-   are equal. */
-      error = !astOK;
-      if( error ) {
-         lstat = astStatus;
-         astClearStatus;
-      }
-      nc = astGetNin( this );
-      if( error ) astSetStatus( lstat );
-
-      for( ic = 0; ic < nc; ic++ ) {
-         (this->jacobian)[ ic ] = astAnnul( (this->jacobian)[ ic ] );
-      }
-      this->jacobian = astFree( this->jacobian );
-   }
-
-/* Free the linear truncation. */
-   if( this->lintrunc ) this->lintrunc = astAnnul( this->lintrunc );
-}
-
-/* Dump function. */
-/* -------------- */
-static void Dump( AstObject *this_object, AstChannel *channel, int *status ) {
-/*
-*  Name:
-*     Dump
-
-*  Purpose:
-*     Dump function for PolyMap objects.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     void Dump( AstObject *this, AstChannel *channel, int *status )
-
-*  Description:
-*     This function implements the Dump function which writes out data
-*     for the PolyMap class to an output Channel.
-
-*  Parameters:
-*     this
-*        Pointer to the PolyMap whose data are being written.
-*     channel
-*        Pointer to the Channel to which the data are being written.
-*     status
-*        Pointer to the inherited status variable.
-*/
-
-#define KEY_LEN 50               /* Maximum length of a keyword */
-
-/* Local Variables: */
-   AstPolyMap *this;             /* Pointer to the PolyMap structure */
-   char buff[ KEY_LEN + 1 ];     /* Buffer for keyword string */
-   char comm[ 100 ];             /* Buffer for comment string */
-   double dval;                  /* Floating point attribute value */
-   int i;                        /* Loop index */
-   int iv;                       /* Vectorised keyword index */
-   int ival;                     /* Integer value */
-   int j;                        /* Loop index */
-   int k;                        /* Loop index */
-   int nin;                      /* No. of input coords */
-   int nout;                     /* No. of output coords */
-   int set;                      /* Attribute value set? */
-
-/* Check the global error status. */
-   if ( !astOK ) return;
-
-/* Obtain a pointer to the PolyMap structure. */
-   this = (AstPolyMap *) this_object;
-
-/* Find the number of inputs and outputs of the uninverted Mapping. */
-   nin = ( (AstMapping *) this )->nin;
-   nout = ( (AstMapping *) this )->nout;
-
-/* Write out values representing the instance variables for the
-   PolyMap class.  */
-
-/* First do the forward transformation arrays. Check they are used. */
-   if( this->ncoeff_f ) {
-
-/* Store the maximum power of each input axis value used by the forward
-   transformation. */
-      for( i = 0; i < nin; i++ ){
-         (void) sprintf( buff, "MPF%d", i + 1 );
-         (void) sprintf( comm, "Max. power of input %d in any forward polynomial", i + 1 );
-          astWriteInt( channel, buff, 1, 1, (this->mxpow_f)[ i ], comm );
-      }
-
-/* Store the number of coefficients associated with each output of the forward
-   transformation. */
-      for( i = 0; i < nout; i++ ){
-         (void) sprintf( buff, "NCF%d", i + 1 );
-         (void) sprintf( comm, "No. of coeff.s for forward polynomial %d", i + 1 );
-         astWriteInt( channel, buff, 1, 1, (this->ncoeff_f)[ i ], comm );
-      }
-
-/* Store the coefficient values used by the forward transformation. */
-      iv = 1;
-      for( i = 0; i < nout; i++ ){
-         for( j = 0; j < this->ncoeff_f[ i ]; j++, iv++ ){
-            if( (this->coeff_f)[ i ][ j ] != AST__BAD ) {
-               (void) sprintf( buff, "CF%d", iv );
-               (void) sprintf( comm, "Coeff %d of forward polynomial %d", j + 1, i + 1 );
-               astWriteDouble( channel, buff, 1, 1, (this->coeff_f)[ i ][ j ], comm );
-            }
-         }
-      }
-
-/* Store the input axis powers associated with each coefficient of the forward
-   transformation. */
-      iv = 1;
-      for( i = 0; i < nout; i++ ){
-         for( j = 0; j < this->ncoeff_f[ i ]; j++ ){
-            for( k = 0; k < nin; k++, iv++ ){
-               if( (this->power_f)[ i ][ j ][ k ] > 0 ) {
-                  (void) sprintf( buff, "PF%d", iv );
-                  (void) sprintf( comm, "Power of i/p %d for coeff %d of fwd poly %d", k + 1, j + 1, i + 1 );
-                  astWriteDouble( channel, buff, 1, 1, (this->power_f)[ i ][ j ][ k ], comm );
-               }
-            }
-         }
-      }
-   }
-
-/* Now do the inverse transformation arrays. Check they are used. */
-   if( this->ncoeff_i ) {
-
-/* Store the maximum power of each output axis value used by the inverse
-   transformation. */
-      for( i = 0; i < nout; i++ ){
-         (void) sprintf( buff, "MPI%d", i + 1 );
-         (void) sprintf( comm, "Max. power of output %d in any inverse polynomial", i + 1 );
-          astWriteInt( channel, buff, 1, 1, (this->mxpow_i)[ i ], comm );
-      }
-
-/* Store the number of coefficients associated with each input of the inverse
-   transformation. */
-      for( i = 0; i < nin; i++ ){
-         (void) sprintf( buff, "NCI%d", i + 1 );
-         (void) sprintf( comm, "No. of coeff.s for inverse polynomial %d", i + 1 );
-         astWriteInt( channel, buff, 1, 1, (this->ncoeff_i)[ i ], comm );
-      }
-
-/* Store the coefficient values used by the inverse transformation. */
-      iv = 1;
-      for( i = 0; i < nin; i++ ){
-         for( j = 0; j < this->ncoeff_i[ i ]; j++, iv++ ){
-            if( (this->coeff_i)[ i ][ j ] != AST__BAD ) {
-               (void) sprintf( buff, "CI%d", iv );
-               (void) sprintf( comm, "Coeff %d of inverse polynomial %d", j + 1, i + 1 );
-               astWriteDouble( channel, buff, 1, 1, (this->coeff_i)[ i ][ j ], comm );
-            }
-         }
-      }
-
-/* Store the output axis powers associated with each coefficient of the inverse
-   transformation. */
-      iv = 1;
-      for( i = 0; i < nin; i++ ){
-         for( j = 0; j < this->ncoeff_i[ i ]; j++ ){
-            for( k = 0; k < nout; k++, iv++ ){
-               if( (this->power_i)[ i ][ j ][ k ] > 0 ) {
-                  (void) sprintf( buff, "PI%d", iv );
-                  (void) sprintf( comm, "Power of o/p %d for coeff %d of inv poly %d", k + 1, j + 1, i + 1 );
-                  astWriteDouble( channel, buff, 1, 1, (this->power_i)[ i ][ j ][ k ], comm );
-               }
-            }
-         }
-      }
-   }
-
-/* Use an iterative inverse? */
-   set = TestIterInverse( this, status );
-   ival = set ? GetIterInverse( this, status ) : astGetIterInverse( this );
-   astWriteInt( channel, "IterInv", set, 0, ival, ival ? "Use an iterative inverse" : "Do not use an iterative inverse" );
-
-/* Max number of iterations for iterative inverse. */
-   set = TestNiterInverse( this, status );
-   ival = set ? GetNiterInverse( this, status ) : astGetNiterInverse( this );
-   astWriteInt( channel, "NiterInv", set, 0, ival, "Max number of iterations for iterative inverse" );
-
-/* Target relative error for iterative inverse. */
-   set = TestTolInverse( this, status );
-   dval = set ? GetTolInverse( this, status ) : astGetTolInverse( this );
-   astWriteDouble( channel, "TolInv", set, 0, dval, "Target relative error for iterative inverse" );
-
-/* Undefine macros local to this function. */
-#undef KEY_LEN
-}
-
-/* Standard class functions. */
-/* ========================= */
-/* Implement the astIsAPolyMap and astCheckPolyMap functions using the macros
-   defined for this purpose in the "object.h" header file. */
-astMAKE_ISA(PolyMap,Mapping)
-astMAKE_CHECK(PolyMap)
-
-AstPolyMap *astPolyMap_( int nin, int nout, int ncoeff_f, const double coeff_f[],
-                         int ncoeff_i, const double coeff_i[], const char *options, int *status, ...){
-/*
-*++
-*  Name:
-c     astPolyMap
-f     AST_POLYMAP
-
-*  Purpose:
-*     Create a PolyMap.
-
-*  Type:
-*     Public function.
-
-*  Synopsis:
-c     #include "polymap.h"
-c     AstPolyMap *astPolyMap( int nin, int nout, int ncoeff_f, const double coeff_f[],
-c                             int ncoeff_i, const double coeff_i[],
-c                             const char *options, ... )
-f     RESULT = AST_POLYMAP( NIN, NOUT, NCOEFF_F, COEFF_F, NCOEFF_I, COEFF_I,
-f                           OPTIONS, STATUS )
-
-*  Class Membership:
-*     PolyMap constructor.
-
-*  Description:
-*     This function creates a new PolyMap and optionally initialises
-*     its attributes.
-*
-*     A PolyMap is a form of Mapping which performs a general polynomial
-*     transformation.  Each output coordinate is a polynomial function of
-*     all the input coordinates. The coefficients are specified separately
-*     for each output coordinate. The forward and inverse transformations
-*     are defined independantly by separate sets of coefficients. If no
-*     inverse transformation is supplied, the default behaviour is to use
-*     an iterative method to evaluate the inverse based only on the forward
-*     transformation (see attribute IterInverse).
-
-*  Parameters:
-c     nin
-f     NIN = INTEGER (Given)
-*        The number of input coordinates.
-c     nout
-f     NOUT = INTEGER (Given)
-*        The number of output coordinates.
-c     ncoeff_f
-f     NCOEFF_F = INTEGER (Given)
-*        The number of non-zero coefficients necessary to define the
-*        forward transformation of the PolyMap. If zero is supplied, the
-*        forward transformation will be undefined.
-c     coeff_f
-f     COEFF_F( * ) = DOUBLE PRECISION (Given)
-*        An array containing
-c        "ncoeff_f*( 2 + nin )" elements. Each group of "2 + nin"
-f        "NCOEFF_F*( 2 + NIN )" elements. Each group of "2 + NIN"
-*        adjacent elements describe a single coefficient of the forward
-*        transformation. Within each such group, the first element is the
-*        coefficient value; the next element is the integer index of the
-*        PolyMap output which uses the coefficient within its defining
-*        polynomial (the first output has index 1); the remaining elements
-*        of the group give the integer powers to use with each input
-*        coordinate value (powers must not be negative, and floating
-*        point values are rounded to the nearest integer).
-c        If "ncoeff_f" is zero, a NULL pointer may be supplied for "coeff_f".
-*
-*        For instance, if the PolyMap has 3 inputs and 2 outputs, each group
-*        consisting of 5 elements, A groups such as "(1.2, 2.0, 1.0, 3.0, 0.0)"
-*        describes a coefficient with value 1.2 which is used within the
-*        definition of output 2. The output value is incremented by the
-*        product of the coefficient value, the value of input coordinate
-*        1 raised to the power 1, and the value of input coordinate 2 raised
-*        to the power 3. Input coordinate 3 is not used since its power is
-*        specified as zero. As another example, the group "(-1.0, 1.0,
-*        0.0, 0.0, 0.0 )" describes adds a constant value -1.0 onto
-*        output 1 (it is a constant value since the power for every input
-*        axis is given as zero).
-*
-c        Each final output coordinate value is the sum of the "ncoeff_f" terms
-c        described by the "ncoeff_f" groups within the supplied array.
-f        Each final output coordinate value is the sum of the "NCOEFF_F" terms
-f        described by the "NCOEFF_F" groups within the supplied array.
-c     ncoeff_i
-f     NCOEFF_I = INTEGER (Given)
-*        The number of non-zero coefficients necessary to define the
-*        inverse transformation of the PolyMap. If zero is supplied, the
-*        default behaviour is to use an iterative method to evaluate the
-*        inverse based only on the forward transformation (see attribute
-*        IterInverse).
-c     coeff_i
-f     COEFF_I( * ) = DOUBLE PRECISION (Given)
-*        An array containing
-c        "ncoeff_i*( 2 + nout )" elements. Each group of "2 + nout"
-f        "NCOEFF_I*( 2 + NOUT )" elements. Each group of "2 + NOUT"
-*        adjacent elements describe a single coefficient of the inverse
-c        transformation, using the same schame as "coeff_f",
-f        transformation, using the same schame as "COEFF_F",
-*        except that "inputs" and "outputs" are transposed.
-c        If "ncoeff_i" is zero, a NULL pointer may be supplied for "coeff_i".
-c     options
-f     OPTIONS = CHARACTER * ( * ) (Given)
-c        Pointer to a null-terminated string containing an optional
-c        comma-separated list of attribute assignments to be used for
-c        initialising the new PolyMap. The syntax used is identical to
-c        that for the astSet function and may include "printf" format
-c        specifiers identified by "%" symbols in the normal way.
-f        A character string containing an optional comma-separated
-f        list of attribute assignments to be used for initialising the
-f        new PolyMap. The syntax used is identical to that for the
-f        AST_SET routine.
-c     ...
-c        If the "options" string contains "%" format specifiers, then
-c        an optional list of additional arguments may follow it in
-c        order to supply values to be substituted for these
-c        specifiers. The rules for supplying these are identical to
-c        those for the astSet function (and for the C "printf"
-c        function).
-f     STATUS = INTEGER (Given and Returned)
-f        The global status.
-
-*  Returned Value:
-c     astPolyMap()
-f     AST_POLYMAP = INTEGER
-*        A pointer to the new PolyMap.
-
-*  Notes:
-*     - A null Object pointer (AST__NULL) will be returned if this
-c     function is invoked with the AST error status set, or if it
-f     function is invoked with STATUS set to an error value, or if it
-*     should fail for any reason.
-*--
-*/
-
-/* Local Variables: */
-   astDECLARE_GLOBALS          /* Pointer to thread-specific global data */
-   AstPolyMap *new;            /* Pointer to new PolyMap */
-   va_list args;               /* Variable argument list */
-
-/* Check the global status. */
-   if ( !astOK ) return NULL;
-
-/* Get a pointer to the thread specific global data structure. */
-   astGET_GLOBALS(NULL);
-
-/* Initialise the PolyMap, allocating memory and initialising the
-   virtual function table as well if necessary. */
-   new = astInitPolyMap( NULL, sizeof( AstPolyMap ), !class_init,
-                         &class_vtab, "PolyMap", nin, nout,
-                         ncoeff_f, coeff_f, ncoeff_i, coeff_i );
-
-/* If successful, note that the virtual function table has been
-   initialised. */
-   if ( astOK ) {
-      class_init = 1;
-
-/* Obtain the variable argument list and pass it along with the options string
-   to the astVSet method to initialise the new PolyMap's attributes. */
-      va_start( args, status );
-      astVSet( new, options, NULL, args );
-      va_end( args );
-
-/* If an error occurred, clean up by deleting the new object. */
-      if ( !astOK ) new = astDelete( new );
-   }
-
-/* Return a pointer to the new PolyMap. */
-   return new;
-}
-
-AstPolyMap *astPolyMapId_( int nin, int nout, int ncoeff_f, const double coeff_f[],
-                           int ncoeff_i, const double coeff_i[], const char *options, ... ){
-/*
-*  Name:
-*     astPolyMapId_
-
-*  Purpose:
-*     Create a PolyMap.
-
-*  Type:
-*     Private function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     AstPolyMap *astPolyMap( int nin, int nout, int ncoeff_f, const double coeff_f[],
-*                             int ncoeff_i, const double coeff_i[],
-*                             const char *options, ... )
-
-*  Class Membership:
-*     PolyMap constructor.
-
-*  Description:
-*     This function implements the external (public) interface to the
-*     astPolyMap constructor function. It returns an ID value (instead
-*     of a true C pointer) to external users, and must be provided
-*     because astPolyMap_ has a variable argument list which cannot be
-*     encapsulated in a macro (where this conversion would otherwise
-*     occur).
-*
-*     The variable argument list also prevents this function from
-*     invoking astPolyMap_ directly, so it must be a re-implementation
-*     of it in all respects, except for the final conversion of the
-*     result to an ID value.
-
-*  Parameters:
-*     As for astPolyMap_.
-
-*  Returned Value:
-*     The ID value associated with the new PolyMap.
-*/
-
-/* Local Variables: */
-   astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
-   AstPolyMap *new;              /* Pointer to new PolyMap */
-   va_list args;                 /* Variable argument list */
-   int *status;                  /* Pointer to inherited status value */
-
-/* Get a pointer to the inherited status value. */
-   status = astGetStatusPtr;
-
-/* Get a pointer to the thread specific global data structure. */
-   astGET_GLOBALS(NULL);
-
-/* Check the global status. */
-   if ( !astOK ) return NULL;
-
-/* Initialise the PolyMap, allocating memory and initialising the
-   virtual function table as well if necessary. */
-   new = astInitPolyMap( NULL, sizeof( AstPolyMap ), !class_init,
-                         &class_vtab, "PolyMap", nin, nout,
-                         ncoeff_f, coeff_f, ncoeff_i, coeff_i );
-
-/* If successful, note that the virtual function table has been
-   initialised. */
-   if ( astOK ) {
-      class_init = 1;
-
-/* Obtain the variable argument list and pass it along with the options string
-   to the astVSet method to initialise the new PolyMap's attributes. */
-      va_start( args, options );
-      astVSet( new, options, NULL, args );
-      va_end( args );
-
-/* If an error occurred, clean up by deleting the new object. */
-      if ( !astOK ) new = astDelete( new );
-   }
-
-/* Return an ID value for the new PolyMap. */
-   return astMakeId( new );
-}
-
-AstPolyMap *astInitPolyMap_( void *mem, size_t size, int init,
-                             AstPolyMapVtab *vtab, const char *name,
-                             int nin, int nout, int ncoeff_f, const double coeff_f[],
-                             int ncoeff_i, const double coeff_i[], int *status ){
-/*
-*+
-*  Name:
-*     astInitPolyMap
-
-*  Purpose:
-*     Initialise a PolyMap.
-
-*  Type:
-*     Protected function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     AstPolyMap *astInitPolyMap( void *mem, size_t size, int init,
-*                                 AstPolyMapVtab *vtab, const char *name,
-*                                 int nin, int nout, int ncoeff_f,
-*                                 const double coeff_f[], int ncoeff_i,
-*                                 const double coeff_i[] )
-
-*  Class Membership:
-*     PolyMap initialiser.
-
-*  Description:
-*     This function is provided for use by class implementations to initialise
-*     a new PolyMap object. It allocates memory (if necessary) to accommodate
-*     the PolyMap plus any additional data associated with the derived class.
-*     It then initialises a PolyMap structure at the start of this memory. If
-*     the "init" flag is set, it also initialises the contents of a virtual
-*     function table for a PolyMap at the start of the memory passed via the
-*     "vtab" parameter.
-
-*  Parameters:
-*     mem
-*        A pointer to the memory in which the PolyMap is to be initialised.
-*        This must be of sufficient size to accommodate the PolyMap data
-*        (sizeof(PolyMap)) plus any data used by the derived class. If a value
-*        of NULL is given, this function will allocate the memory itself using
-*        the "size" parameter to determine its size.
-*     size
-*        The amount of memory used by the PolyMap (plus derived class data).
-*        This will be used to allocate memory if a value of NULL is given for
-*        the "mem" parameter. This value is also stored in the PolyMap
-*        structure, so a valid value must be supplied even if not required for
-*        allocating memory.
-*     init
-*        A logical flag indicating if the PolyMap's virtual function table is
-*        to be initialised. If this value is non-zero, the virtual function
-*        table will be initialised by this function.
-*     vtab
-*        Pointer to the start of the virtual function table to be associated
-*        with the new PolyMap.
-*     name
-*        Pointer to a constant null-terminated character string which contains
-*        the name of the class to which the new object belongs (it is this
-*        pointer value that will subsequently be returned by the astGetClass
-*        method).
-*     nin
-*        The number of input coordinate values per point. This is the
-*        same as the number of columns in the matrix.
-*     nout
-*        The number of output coordinate values per point. This is the
-*        same as the number of rows in the matrix.
-*     ncoeff_f
-*        The number of non-zero coefficients necessary to define the
-*        forward transformation of the PolyMap. If zero is supplied, the
-*        forward transformation will be undefined.
-*     coeff_f
-*        An array containing "ncoeff_f*( 2 + nin )" elements. Each group
-*	 of "2 + nin" adjacent elements describe a single coefficient of
-*	 the forward transformation. Within each such group, the first
-*	 element is the coefficient value; the next element is the
-*	 integer index of the PolyMap output which uses the coefficient
-*	 within its defining polynomial (the first output has index 1);
-*	 the remaining elements of the group give the integer powers to
-*	 use with each input coordinate value (powers must not be
-*	 negative)
-*
-*        For instance, if the PolyMap has 3 inputs and 2 outputs, each group
-*        consisting of 5 elements, A groups such as "(1.2, 2.0, 1.0, 3.0, 0.0)"
-*        describes a coefficient with value 1.2 which is used within the
-*        definition of output 2. The output value is incremented by the
-*        product of the coefficient value, the value of input coordinate
-*        1 raised to the power 1, and the value of input coordinate 2 raised
-*        to the power 3. Input coordinate 3 is not used since its power is
-*        specified as zero. As another example, the group "(-1.0, 1.0,
-*        0.0, 0.0, 0.0 )" describes adds a constant value -1.0 onto
-*        output 1 (it is a constant value since the power for every input
-*        axis is given as zero).
-*
-*        Each final output coordinate value is the sum of the "ncoeff_f" terms
-*        described by the "ncoeff_f" groups within the supplied array.
-*     ncoeff_i
-*        The number of non-zero coefficients necessary to define the
-*        inverse transformation of the PolyMap. If zero is supplied, the
-*        inverse transformation will be undefined.
-*     coeff_i
-*        An array containing
-*        "ncoeff_i*( 2 + nout )" elements. Each group of "2 + nout"
-*        adjacent elements describe a single coefficient of the inverse
-*        transformation, using the same schame as "coeff_f", except that
-*        "inputs" and "outputs" are transposed.
-
-*  Returned Value:
-*     A pointer to the new PolyMap.
-
-*  Notes:
-*     -  A null pointer will be returned if this function is invoked with the
-*     global error status set, or if it should fail for any reason.
-*-
-*/
-
-/* Local Variables: */
-   AstPolyMap *new;              /* Pointer to new PolyMap */
-
-/* Check the global status. */
-   if ( !astOK ) return NULL;
-
-/* If necessary, initialise the virtual function table. */
-   if ( init ) astInitPolyMapVtab( vtab, name );
-
-/* Initialise a Mapping structure (the parent class) as the first component
-   within the PolyMap structure, allocating memory if necessary. Specify that
-   the Mapping should be defined in both the forward and inverse directions. */
-   new = (AstPolyMap *) astInitMapping( mem, size, 0,
-                                        (AstMappingVtab *) vtab, name,
-                                        nin, nout, 1, 1 );
-   if ( astOK ) {
-
-/* Initialise the PolyMap data. */
-/* ---------------------------- */
-
-/* First initialise the pointers in case of errors. */
-      new->ncoeff_f = NULL;
-      new->power_f = NULL;
-      new->coeff_f = NULL;
-      new->mxpow_f = NULL;
-
-      new->ncoeff_i = NULL;
-      new->power_i = NULL;
-      new->coeff_i = NULL;
-      new->mxpow_i = NULL;
-
-/* Store the forward transformation. */
-      StoreArrays( new, 1, ncoeff_f, coeff_f, status );
-
-/* Store the inverse transformation. */
-      StoreArrays( new, 0, ncoeff_i, coeff_i, status );
-
-/* Other class attributes. */
-      new->iterinverse = -INT_MAX;
-      new->niterinverse = -INT_MAX;
-      new->tolinverse = AST__BAD;
-      new->jacobian = NULL;
-      new->lintrunc = NULL;
-
-/* If an error occurred, clean up by deleting the new PolyMap. */
-      if ( !astOK ) new = astDelete( new );
-   }
-
-/* Return a pointer to the new PolyMap. */
-   return new;
-}
-
-AstPolyMap *astLoadPolyMap_( void *mem, size_t size,
-                             AstPolyMapVtab *vtab, const char *name,
-                             AstChannel *channel, int *status ) {
-/*
-*+
-*  Name:
-*     astLoadPolyMap
-
-*  Purpose:
-*     Load a PolyMap.
-
-*  Type:
-*     Protected function.
-
-*  Synopsis:
-*     #include "polymap.h"
-*     AstPolyMap *astLoadPolyMap( void *mem, size_t size,
-*                                 AstPolyMapVtab *vtab, const char *name,
-*                                 AstChannel *channel )
-
-*  Class Membership:
-*     PolyMap loader.
-
-*  Description:
-*     This function is provided to load a new PolyMap using data read
-*     from a Channel. It first loads the data used by the parent class
-*     (which allocates memory if necessary) and then initialises a
-*     PolyMap structure in this memory, using data read from the input
-*     Channel.
-*
-*     If the "init" flag is set, it also initialises the contents of a
-*     virtual function table for a PolyMap at the start of the memory
-*     passed via the "vtab" parameter.
-
-
-*  Parameters:
-*     mem
-*        A pointer to the memory into which the PolyMap is to be
-*        loaded.  This must be of sufficient size to accommodate the
-*        PolyMap data (sizeof(PolyMap)) plus any data used by derived
-*        classes. If a value of NULL is given, this function will
-*        allocate the memory itself using the "size" parameter to
-*        determine its size.
-*     size
-*        The amount of memory used by the PolyMap (plus derived class
-*        data).  This will be used to allocate memory if a value of
-*        NULL is given for the "mem" parameter. This value is also
-*        stored in the PolyMap structure, so a valid value must be
-*        supplied even if not required for allocating memory.
-*
-*        If the "vtab" parameter is NULL, the "size" value is ignored
-*        and sizeof(AstPolyMap) is used instead.
-*     vtab
-*        Pointer to the start of the virtual function table to be
-*        associated with the new PolyMap. If this is NULL, a pointer
-*        to the (static) virtual function table for the PolyMap class
-*        is used instead.
-*     name
-*        Pointer to a constant null-terminated character string which
-*        contains the name of the class to which the new object
-*        belongs (it is this pointer value that will subsequently be
-*        returned by the astGetClass method).
-*
-*        If the "vtab" parameter is NULL, the "name" value is ignored
-*        and a pointer to the string "PolyMap" is used instead.
-
-*  Returned Value:
-*     A pointer to the new PolyMap.
-
-*  Notes:
-*     - A null pointer will be returned if this function is invoked
-*     with the global error status set, or if it should fail for any
-*     reason.
-*-
-*/
-
-#define KEY_LEN 50               /* Maximum length of a keyword */
-
-   astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
-/* Local Variables: */
-   AstPolyMap *new;              /* Pointer to the new PolyMap */
-   char buff[ KEY_LEN + 1 ];     /* Buffer for keyword string */
-   int i;                        /* Loop index */
-   int iv;                       /* Vectorised keyword index */
-   int j;                        /* Loop index */
-   int k;                        /* Loop index */
-   int nin;                      /* No. of input coords */
-   int nout;                     /* No. of output coords */
-   int undef;                    /* Is the transformation undefined? */
-
-/* Get a pointer to the thread specific global data structure. */
-   astGET_GLOBALS(channel);
-
-/* Initialise. */
-   new = NULL;
-
-/* Check the global error status. */
-   if ( !astOK ) return new;
-
-/* If a NULL virtual function table has been supplied, then this is
-   the first loader to be invoked for this PolyMap. In this case the
-   PolyMap belongs to this class, so supply appropriate values to be
-   passed to the parent class loader (and its parent, etc.). */
-   if ( !vtab ) {
-      size = sizeof( AstPolyMap );
-      vtab = &class_vtab;
-      name = "PolyMap";
-
-/* If required, initialise the virtual function table for this class. */
-      if ( !class_init ) {
-         astInitPolyMapVtab( vtab, name );
-         class_init = 1;
-      }
-   }
-
-/* Invoke the parent class loader to load data for all the ancestral
-   classes of the current one, returning a pointer to the resulting
-   partly-built PolyMap. */
-   new = astLoadMapping( mem, size, (AstMappingVtab *) vtab, name,
-                         channel );
-
-   if ( astOK ) {
-
-/* Get the number of inputs and outputs for the uninverted Mapping. */
-   nin = ( (AstMapping *) new )->nin;
-   nout = ( (AstMapping *) new )->nout;
-
-/* Read input data. */
-/* ================ */
-/* Request the input Channel to read all the input data appropriate to
-   this class into the internal "values list". */
-      astReadClassData( channel, "PolyMap" );
-
-/* Allocate memory to hold the forward arrays. */
-      new->ncoeff_f = astMalloc( sizeof( int )*(size_t) nout );
-      new->mxpow_f = astMalloc( sizeof( int )*(size_t) nin );
-      new->power_f = astMalloc( sizeof( int ** )*(size_t) nout );
-      new->coeff_f = astMalloc( sizeof( double * )*(size_t) nout );
-      if( astOK ) {
-
-/* Assume the forward transformation is defined. */
-         undef = 0;
-
-/* Get the maximum power of each input axis value used by the forward
-   transformation. Set a flag "undef" if no values relating to the
-   forward transformation are found (this indicates that the forward
-   transformation is not defined). */
-         for( i = 0; i < nin && !undef; i++ ){
-            (void) sprintf( buff, "mpf%d", i + 1 );
-            (new->mxpow_f)[ i ] = astReadInt( channel, buff, INT_MAX );
-            if( (new->mxpow_f)[ i ] == INT_MAX ) undef = 1;
-         }
-
-/* Get the number of coefficients associated with each output of the forward
-   transformation. */
-         for( i = 0; i < nout && !undef; i++ ){
-            (void) sprintf( buff, "ncf%d", i + 1 );
-            (new->ncoeff_f)[ i ] = astReadInt( channel, buff, INT_MAX );
-            if( (new->ncoeff_f)[ i ] == INT_MAX ) undef = 1;
-         }
-
-/* Get the coefficient values used by the forward transformation. This
-   uses new style vectorised key names if available. Otherwise it uses
-   old style indexed names (which were superceded by vectorised names
-   because they are shorter and so work better with FitsChans). */
-         iv = 0;
-         for( i = 0; i < nout && !undef; i++ ){
-
-            (new->coeff_f)[ i ] = astMalloc( sizeof( double )*
-                                           (size_t) new->ncoeff_f[ i ] );
-            if( astOK ) {
-               for( j = 0; j < new->ncoeff_f[ i ]; j++ ){
-                  (void) sprintf( buff, "cf%d", ++iv );
-                  (new->coeff_f)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD );
-                  if( (new->coeff_f)[ i ][ j ] == AST__BAD ) {
-                     (void) sprintf( buff, "cf%d_%d", i + 1, j + 1 );
-                     (new->coeff_f)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD );
-                  }
-               }
-            }
-         }
-
-/* Get the input axis powers associated with each coefficient of the forward
-   transformation. */
-         iv = 0;
-         for( i = 0; i < nout && !undef; i++ ){
-            (new->power_f)[ i ] = astMalloc( sizeof( int * )*
-                                             (size_t) new->ncoeff_f[ i ] );
-            if( astOK ) {
-               for( j = 0; j < new->ncoeff_f[ i ]; j++ ){
-                  (new->power_f)[ i ][ j ] = astMalloc( sizeof( int )* (size_t) nin );
-                  if( astOK ) {
-                     for( k = 0; k < nin; k++ ){
-                        (void) sprintf( buff, "pf%d", ++iv );
-                        (new->power_f)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 );
-                        if( (new->power_f)[ i ][ j ][ k ] == 0 ) {
-                           (void) sprintf( buff, "pf%d_%d_%d", i + 1, j + 1, k + 1 );
-                           (new->power_f)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 );
-                        }
-                     }
-                  }
-               }
-            }
-         }
-
-/* Free the arrays if the forward transformation is undefined. */
-         if( undef ) {
-            new->ncoeff_f = astFree( new->ncoeff_f );
-            new->mxpow_f = astFree( new->mxpow_f );
-            new->power_f = astFree( new->power_f );
-            new->coeff_f = astFree( new->coeff_f );
-         }
-      }
-
-/* Allocate memory to hold the inverse arrays. */
-      new->ncoeff_i = astMalloc( sizeof( int )*(size_t) nin );
-      new->mxpow_i = astMalloc( sizeof( int )*(size_t) nout );
-      new->power_i = astMalloc( sizeof( int ** )*(size_t) nin );
-      new->coeff_i = astMalloc( sizeof( double * )*(size_t) nin );
-      if( astOK ) {
-
-/* Assume the inverse transformation is defined. */
-         undef = 0;
-
-/* Get the maximum power of each output axis value used by the inverse
-   transformation. Set a flag "undef" if no values relating to the
-   inverse transformation are found (this indicates that the inverse
-   transformation is not defined). */
-         for( i = 0; i < nout && !undef; i++ ){
-            (void) sprintf( buff, "mpi%d", i + 1 );
-            (new->mxpow_i)[ i ] = astReadInt( channel, buff, INT_MAX );
-            if( (new->mxpow_i)[ i ] == INT_MAX ) undef = 1;
-         }
-
-/* Get the number of coefficients associated with each input of the inverse
-   transformation. */
-         for( i = 0; i < nin && !undef; i++ ){
-            (void) sprintf( buff, "nci%d", i + 1 );
-            (new->ncoeff_i)[ i ] = astReadInt( channel, buff, INT_MAX );
-            if( (new->ncoeff_i)[ i ] == INT_MAX ) undef = 1;
-         }
-
-/* Get the coefficient values used by the inverse transformation. */
-         iv = 0;
-         for( i = 0; i < nin && !undef; i++ ){
-
-            (new->coeff_i)[ i ] = astMalloc( sizeof( double )*
-                                           (size_t) new->ncoeff_i[ i ] );
-            if( astOK ) {
-               for( j = 0; j < new->ncoeff_i[ i ]; j++ ){
-                  (void) sprintf( buff, "ci%d", ++iv );
-                  (new->coeff_i)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD );
-                  if( (new->coeff_i)[ i ][ j ] == AST__BAD ) {
-                     (void) sprintf( buff, "ci%d_%d", i + 1, j + 1 );
-                     (new->coeff_i)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD );
-                  }
-               }
-            }
-         }
-
-/* Get the output axis powers associated with each coefficient of the inverse
-   transformation. */
-         iv = 0;
-         for( i = 0; i < nin && !undef; i++ ){
-            (new->power_i)[ i ] = astMalloc( sizeof( int * )*
-                                             (size_t) new->ncoeff_i[ i ] );
-            if( astOK ) {
-               for( j = 0; j < new->ncoeff_i[ i ]; j++ ){
-                  (new->power_i)[ i ][ j ] = astMalloc( sizeof( int )* (size_t) nout );
-                  if( astOK ) {
-                     for( k = 0; k < nout; k++ ){
-                        (void) sprintf( buff, "pi%d", ++iv );
-                        (new->power_i)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 );
-                        if( (new->power_i)[ i ][ j ][ k ] == 0 ) {
-                           (void) sprintf( buff, "pi%d_%d_%d", i + 1, j + 1, k + 1 );
-                           (new->power_i)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 );
-                        }
-                     }
-                  }
-               }
-            }
-         }
-
-/* Free the arrays if the inverse transformation is undefined. */
-         if( undef ) {
-            new->ncoeff_i = astFree( new->ncoeff_i );
-            new->mxpow_i = astFree( new->mxpow_i );
-            new->power_i = astFree( new->power_i );
-            new->coeff_i = astFree( new->coeff_i );
-         }
-      }
-
-/* Whether to use an iterative inverse transformation. */
-      new->iterinverse = astReadInt( channel, "iterinv", -INT_MAX );
-      if ( TestIterInverse( new, status ) ) SetIterInverse( new, new->iterinverse, status );
-
-/* Max number of iterations for iterative inverse transformation. */
-      new->niterinverse = astReadInt( channel, "niterinv", -INT_MAX );
-      if ( TestNiterInverse( new, status ) ) SetNiterInverse( new, new->niterinverse, status );
-
-/* Target relative error for iterative inverse transformation. */
-      new->tolinverse = astReadDouble( channel, "tolinv", AST__BAD );
-      if ( TestTolInverse( new, status ) ) SetTolInverse( new, new->tolinverse, status );
-
-/* The Jacobian of the PolyMap's forward transformation has not yet been
-   found. */
-      new->jacobian = NULL;
-
-/* The linear truncation of the PolyMap has not yet been found. */
-      new->lintrunc = NULL;
-
-/* If an error occurred, clean up by deleting the new PolyMap. */
-      if ( !astOK ) new = astDelete( new );
-   }
-
-/* Return the new PolyMap pointer. */
-   return new;
-
-/* Undefine macros local to this function. */
-#undef KEY_LEN
-}
-
-/* Virtual function interfaces. */
-/* ============================ */
-/* These provide the external interface to the virtual functions defined by
-   this class. Each simply checks the global error status and then locates and
-   executes the appropriate member function, using the function pointer stored
-   in the object's virtual function table (this pointer is located using the
-   astMEMBER macro defined in "object.h").
-
-   Note that the member function may not be the one defined here, as it may
-   have been over-ridden by a derived class. However, it should still have the
-   same interface. */
-
-void astPolyPowers_( AstPolyMap *this, double **work, int ncoord,
-                     const int *mxpow, double **ptr, int point, int fwd,
-                     int *status ){
-   if ( !astOK ) return;
-   (**astMEMBER(this,PolyMap,PolyPowers))( this, work, ncoord, mxpow, ptr,
-                                           point, fwd, status );
-}
-
-AstPolyMap *astPolyTran_( AstPolyMap *this, int forward, double acc,
-                          double maxacc, int maxorder, const double *lbnd,
-                          const double *ubnd, int *status ){
-   if ( !astOK ) return NULL;
-   return (**astMEMBER(this,PolyMap,PolyTran))( this, forward, acc,
-                                                maxacc, maxorder, lbnd,
-                                                ubnd, status );
-}
-
-void astPolyCoeffs_( AstPolyMap *this, int forward, int nel, double *array,
-                     int *ncoeff, int *status ){
-   if ( !astOK ) return;
-   (**astMEMBER(this,PolyMap,PolyCoeffs))( this, forward, nel,
-                                           array, ncoeff, status );
-}
-
-void astFitPoly1DInit_( AstPolyMap *this, int forward, double **table,
-                        AstMinPackData *data, double *scales,
-                        int *status ){
-   if ( !astOK ) return;
-   (**astMEMBER(this,PolyMap,FitPoly1DInit))( this, forward, table, data, scales,
-                                              status );
-}
-
-
-void astFitPoly2DInit_( AstPolyMap *this, int forward, double **table,
-                        AstMinPackData *data, double *scales,
-                        int *status ){
-   if ( !astOK ) return;
-   (**astMEMBER(this,PolyMap,FitPoly2DInit))( this, forward, table, data, scales,
-                                              status );
-}
-
-
-
-