update ast 8.6.2

1 files changed, 6107 insertions, 0 deletions

diff --git a/ast/polymap.c b/ast/polymap.c
new file mode 100644
index 0000000..498000f
--- /dev/null
+++ b/ast/polymap.c
@@ -0,0 +1,6107 @@
+/*
+*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.
+*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 "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 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 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 PolyPowers( AstPolyMap *, double **, int, const int *, double **, int, int, int * );
+static void StoreArrays( AstPolyMap *, int, int, const double *, int * );
+static void PolyCoeffs( AstPolyMap *, int, int, double *, int *, int * );
+static void FitPoly1DInit( AstPolyMap *, int, double **, AstMinPackData *, double *, int *);
+static void FitPoly2DInit( AstPolyMap *, int, double **, AstMinPackData *, 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 ) {
+
+      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 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 forward transformation of the supplied PolyMap.
+*
+*     The Jacobian matrix has "nout" rows and "nin" columns, where "nin"
+*     and "nout" are the number of inputs and outputs of the supplied 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. */
+      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 inverse transformation of a
+*     PolyMap at a set of output positions.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     void IterInverse( AstPolyMap *this, AstPointSet *out, AstPointSet *result,
+*                       int *status )
+
+*  Description:
+*     This function transforms a set of PolyMap output positions using
+*     the inverse transformation of the PolyMap, to generate the corresponding
+*     input positions. An iterative Newton-Raphson method is used which
+*     only required the forward transformation of the PolyMap to be deifned.
+
+*  Parameters:
+*     this
+*        The PolyMap.
+*     out
+*        A PointSet holding the PolyMap output positions that are to be
+*        transformed using the inverse transformation.
+*     result
+*        A PointSet into which the generated PolyMap input positions are to be
+*        stored.
+*     status
+*        Pointer to the inherited status variable.
+
+*/
+
+/* Local Variables: */
+   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. We assume initially
+   that the inverse transformation is a unit mapping, and so we just copy
+   the supplied outputs positions to the results PointSet holding the
+   corresponding input positions. */
+      for( icoord = 0; icoord < ncoord; icoord++ ) {
+         memcpy( ptr_in[ icoord ], ptr_out[ icoord ], sizeof( double )*npoint );
+      }
+
+/* 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 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 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:
+*     - 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 );
+
+/* 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 inputs */
+   int nout;                     /* Number of outputs */
+   int pow;                      /* Power extracted from coeff. description */
+
+/* Check the global status. */
+   if ( !astOK ) return;
+
+/* Get the number of inputs and outputs. */
+   nin = astGetNin( this );
+   nout = astGetNout( this );
+
+/* First Free any existing arrays. */
+   FreeArrays( this, forward, status );
+
+/* Now initialise the 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 withint he 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 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 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.
+*
+*     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;
+
+/* 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. */
+   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 );
+                  }
+               }
+            }
+         }
+      }
+   }
+
+/* If an error has occurred, free al 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 );
+   }
+}
+
+/* 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;
+
+/* 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;
+
+/* 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 );
+}
+
+
+
+