path: root/ast/polymap.c

/*
*class++
*  Name:
*     PolyMap

*  Purpose:
*     Map coordinates using polynomial functions.

*  Constructor Function:
c     astPolyMap
f     AST_POLYMAP

*  Description:
*     A PolyMap is a form of Mapping which performs a general polynomial
*     transformation.  Each output coordinate is a polynomial function of
*     all the input coordinates. The coefficients are specified separately
*     for each output coordinate. The forward and inverse transformations
*     are defined independantly by separate sets of coefficients. If no
*     inverse transformation is supplied, the default behaviour is to use
*     an iterative method to evaluate the inverse based only on the forward
*     transformation (see attribute IterInverse).

*  Inheritance:
*     The PolyMap class inherits from the Mapping class.

*  Attributes:
*     In addition to those attributes common to all Mappings, every
*     PolyMap also has the following attributes:
*
*     - IterInverse: Provide an iterative inverse transformation?
*     - NiterInverse: Maximum number of iterations for iterative inverse
*     - TolInverse: Target relative error for iterative inverse

*  Functions:
c     In addition to those functions applicable to all Objects, the
c     following functions may also be applied to all Mappings:
f     In addition to those routines applicable to all Objects, the
f     following routines may also be applied to all Mappings:
*
c     - astPolyCoeffs: Retrieve the coefficients of a PolyMap transformation
c     - astPolyTran: Fit a PolyMap inverse or forward transformation
f     - AST_POLYCOEFFS: Retrieve the coefficients of a PolyMap transformation
f     - AST_POLYTRAN: Fit a PolyMap inverse or forward transformation

*  Copyright:
*     Copyright (C) 2017 East Asian Observatory.
*     Copyright (C) 1997-2006 Council for the Central Laboratory of the
*     Research Councils
*     Copyright (C) 2011 Science & Technology Facilities Council.
*     All Rights Reserved.

*  Licence:
*     This program is free software: you can redistribute it and/or
*     modify it under the terms of the GNU Lesser General Public
*     License as published by the Free Software Foundation, either
*     version 3 of the License, or (at your option) any later
*     version.
*
*     This program is distributed in the hope that it will be useful,
*     but WITHOUT ANY WARRANTY; without even the implied warranty of
*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*     GNU Lesser General Public License for more details.
*
*     You should have received a copy of the GNU Lesser General
*     License along with this program.  If not, see
*     <http://www.gnu.org/licenses/>.

*  Authors:
*     DSB: D.S. Berry (Starlink)

*  History:
*     27-SEP-2003 (DSB):
*        Original version.
*     13-APR-2005 (DSB):
*        Changed the keys used by the Dump/astLoadPolyMap functions. They
*        used to exceed 8 characters and consequently caused problems for
*        FitsChans.
*     20-MAY-2005 (DSB):
*        Correct the indexing of keywords produced in the Dump function.
*     20-APR-2006 (DSB):
*        Guard against undefined transformations in Copy.
*     10-MAY-2006 (DSB):
*        Override astEqual.
*     4-JUL-2008 (DSB):
*        Fixed loop indexing problems in Equal function.
*     27-MAY-2011 (DSB):
*        Added public method astPolyTran.
*     18-JUL-2011 (DSB):
*        - Added attributes IterInverse, NiterInverse and TolInverse.
*        - Do not report an error if astPolyTran fails to fit an inverse.
*     15-OCT-2011 (DSB):
*        Improve argument checking and error reporting in PolyTran
*     8-MAY-2014 (DSB):
*        Move to using CMinPack for minimisations.
*     11-NOV-2016 (DSB):
*        - Fix bug in MapMerge that could cause a seg fault. It did not
*        check that the PolyMap had a defined transformation before accessing
*        the transformation's coefficient array.
*        - Fix similar bugs in Equal that could cause seg faults.
*     15-MAR-2017 (DSB):
*        - Change the GetTranForward and GetTranInverse functions so that they
*        take into account the state of the Invert attribute.
*        - Improve docs for the IterInverse attribute to explain that the
*        inverse transformation replaced is always the original inverse
*        transformation, as defined by the arguments supplied to the PolyMap
*        constructor, regardless of the state of the Invert attribute.
*     17-MAR-2017 (DSB):
*        - Add the astPolyCoeffs method.
*     30-MAR-2017 (DSB):
*        Modify the astPolyTran method so that it can be used by the
*        ChebyMap class to determine new transformations implemented as
*        Chebyshev polynomials.
*     27-JUN-2017 (DSB):
*        In SamplePoly1D/2D ensure the final sample on each axis does not
*        go above the supplied upper bound. This can happen due to rounding
*        error. This is important for ChebyMaps since points outside the
*        bounds are set bad when transformed using a ChebyMap, causing NaNs
*        to be generated in lmder1 (cminpack minimisation function).
*     3-JUL-2017 (DSB):
*        Within FitPoly1D and FitPoly2D, use an initial guess that represents
*        a unit mapping between normalised input and output values, rather
*        than unnormalised PolyMap values. This is in case the PolyMap being
*        fitted includes a change of scale (e.g. the PolyMap input is in "mm"
*        but the output is in "rads" and includes some large scaling factor
*        to do the conversion).
*     9-JAN-2018 (DSB):
*        Correct Transform to take account of AST__BAD coeffs correctly.
*        Previously bad coeffs could generate NaN output values.
*     27-APR-2018 (DSB):
*        When calculating the iterative inverse, use an initial guess based
*        on the linear truncation of the PolyMap rather than a UnitMap.
*     3-MAY-2018 (DSB):
*        Clear the IterInverse attribute in the PolyMap returned by
*        astPolyTran (why would you want to use the iterative inverse if
*        you have just created an inverse fit?).
*     4-MAY-2018 (DSB):
*        Fix various places where there was confusion over whether the
*        "forward transformation" refers to the forward transformation
*        of the original uninverted PolyMap, or the current forward
*        transformation of the PolyMap (i.e. taking the "Invert" flag into
*        account).
*class--
*/

/* Module Macros. */
/* ============== */
/* Set the name of the class we are implementing. This indicates to
   the header files that define class interfaces that they should make
   "protected" symbols available. */
#define astCLASS PolyMap

/* Include files. */
/* ============== */
/* Interface definitions. */
/* ---------------------- */

#include "globals.h"             /* Thread-safe global data access */
#include "error.h"               /* Error reporting facilities */
#include "memory.h"              /* Memory allocation facilities */
#include "object.h"              /* Base Object class */
#include "pointset.h"            /* Sets of points/coordinates */
#include "matrixmap.h"           /* Matrix multiplication mappings */
#include "shiftmap.h"            /* Shift of origin mappings */
#include "mapping.h"             /* Coordinate mappings (parent class) */
#include "cmpmap.h"              /* Compound mappings */
#include "polymap.h"             /* Interface definition for this class */
#include "unitmap.h"             /* Unit mappings */
#include "cminpack/cminpack.h"   /* Levenberg - Marquardt minimization */
#include "pal.h"                 /* SLALIB function definitions */

/* Error code definitions. */
/* ----------------------- */
#include "ast_err.h"             /* AST error codes */

/* C header files. */
/* --------------- */
#include <ctype.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <float.h>

/* Module Variables. */
/* ================= */

/* Address of this static variable is used as a unique identifier for
   member of this class. */
static int class_check;

/* Pointers to parent class methods which are extended by this class. */
static AstPointSet *(* parent_transform)( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
static const char *(* parent_getattrib)( AstObject *, const char *, int * );
static int (* parent_testattrib)( AstObject *, const char *, int * );
static void (* parent_clearattrib)( AstObject *, const char *, int * );
static void (* parent_setattrib)( AstObject *, const char *, int * );
static int (* parent_getobjsize)( AstObject *, int * );

#if defined(THREAD_SAFE)
static int (* parent_managelock)( AstObject *, int, int, AstObject **, int * );
#endif


#ifdef THREAD_SAFE
/* Define how to initialise thread-specific globals. */
#define GLOBAL_inits \
   globals->Class_Init = 0; \
   globals->GetAttrib_Buff[ 0 ] = 0;

/* Create the function that initialises global data for this module. */
astMAKE_INITGLOBALS(PolyMap)

/* Define macros for accessing each item of thread specific global data. */
#define class_init astGLOBAL(PolyMap,Class_Init)
#define class_vtab astGLOBAL(PolyMap,Class_Vtab)
#define getattrib_buff astGLOBAL(LutMap,GetAttrib_Buff)

#include <pthread.h>


#else

static char getattrib_buff[ 101 ];

/* Define the class virtual function table and its initialisation flag
   as static variables. */
static AstPolyMapVtab class_vtab;   /* Virtual function table */
static int class_init = 0;       /* Virtual function table initialised? */

#endif



/* External Interface Function Prototypes. */
/* ======================================= */
/* The following functions have public prototypes only (i.e. no
   protected prototypes), so we must provide local prototypes for use
   within this module. */
AstPolyMap *astPolyMapId_( int, int, int, const double[], int, const double[], const char *, ... );


/* Prototypes for Private Member Functions. */
/* ======================================== */
static AstMapping *LinearGuess( AstPolyMap *, int * );
static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
static AstPolyMap **GetJacobian( AstPolyMap *, int * );
static AstPolyMap *PolyTran( AstPolyMap *, int, double, double, int, const double *, const double *, int * );
static double **SamplePoly1D( AstPolyMap *, int, double **, double, double, int, int *, double[2], int * );
static double **SamplePoly2D( AstPolyMap *, int, double **, const double *, const double *, int, int *, double[4], int * );
static double *FitPoly1D( AstPolyMap *, int, int, double, int, double **, double[2], int *, double *, int * );
static double *FitPoly2D( AstPolyMap *, int, int, double, int, double **, double[4], int *, double *, int * );
static int Equal( AstObject *, AstObject *, int * );
static int GetObjSize( AstObject *, int * );
static int GetTranForward( AstMapping *, int * );
static int GetTranInverse( AstMapping *, int * );
static int MPFunc1D( void *, int, int, const double *, double *, double *, int, int );
static int MPFunc2D( void *, int, int, const double *, double *, double *, int, int );
static int MapMerge( AstMapping *, int, int, int *, AstMapping ***, int **, int * );
static int ReplaceTransformation( AstPolyMap *, int, double, double, int, const double *, const double *, int * );
static void Copy( const AstObject *, AstObject *, int * );
static void Delete( AstObject *obj, int * );
static void Dump( AstObject *, AstChannel *, int * );
static void FitPoly1DInit( AstPolyMap *, int, double **, AstMinPackData *, double *, int *);
static void FitPoly2DInit( AstPolyMap *, int, double **, AstMinPackData *, double *, int *);
static void FreeArrays( AstPolyMap *, int, int * );
static void IterInverse( AstPolyMap *, AstPointSet *, AstPointSet *, int * );
static void LMFunc1D(  const double *, double *, int, int, void * );
static void LMFunc2D(  const double *, double *, int, int, void * );
static void LMJacob1D( const double *, double *, int, int, void * );
static void LMJacob2D( const double *, double *, int, int, void * );
static void PolyCoeffs( AstPolyMap *, int, int, double *, int *, int * );
static void PolyPowers( AstPolyMap *, double **, int, const int *, double **, int, int, int * );
static void StoreArrays( AstPolyMap *, int, int, const double *, int * );

#if defined(THREAD_SAFE)
static int ManageLock( AstObject *, int, int, AstObject **, int * );
#endif

static const char *GetAttrib( AstObject *, const char *, int * );
static int TestAttrib( AstObject *, const char *, int * );
static void ClearAttrib( AstObject *, const char *, int * );
static void SetAttrib( AstObject *, const char *, int * );

static int GetIterInverse( AstPolyMap *, int * );
static int TestIterInverse( AstPolyMap *, int * );
static void ClearIterInverse( AstPolyMap *, int * );
static void SetIterInverse( AstPolyMap *, int, int * );

static int GetNiterInverse( AstPolyMap *, int * );
static int TestNiterInverse( AstPolyMap *, int * );
static void ClearNiterInverse( AstPolyMap *, int * );
static void SetNiterInverse( AstPolyMap *, int, int * );

static double GetTolInverse( AstPolyMap *, int * );
static int TestTolInverse( AstPolyMap *, int * );
static void ClearTolInverse( AstPolyMap *, int * );
static void SetTolInverse( AstPolyMap *, double, int * );


/* Member functions. */
/* ================= */

static void ClearAttrib( AstObject *this_object, const char *attrib, int *status ) {
/*
*  Name:
*     ClearAttrib

*  Purpose:
*     Clear an attribute value for a PolyMap.

*  Type:
*     Private function.

*  Synopsis:
*     #include "polymap.h"
*     void ClearAttrib( AstObject *this, const char *attrib, int *status )

*  Class Membership:
*     PolyMap member function (over-rides the astClearAttrib protected
*     method inherited from the Mapping class).

*  Description:
*     This function clears the value of a specified attribute for a
*     PolyMap, so that the default value will subsequently be used.

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     attrib
*        Pointer to a null-terminated string specifying the attribute
*        name.  This should be in lower case with no surrounding white
*        space.
*     status
*        Pointer to the inherited status variable.
*/

/* Local Variables: */
   AstPolyMap *this;             /* Pointer to the PolyMap structure */

/* Check the global error status. */
   if ( !astOK ) return;

/* Obtain a pointer to the PolyMap structure. */
   this = (AstPolyMap *) this_object;

/* Check the attribute name and clear the appropriate attribute. */

/* IterInverse. */
/* ------------ */
   if ( !strcmp( attrib, "iterinverse" ) ) {
      astClearIterInverse( this );

/* NiterInverse. */
/* ------------- */
   } else if ( !strcmp( attrib, "niterinverse" ) ) {
      astClearNiterInverse( this );

/* TolInverse. */
/* ----------- */
   } else if ( !strcmp( attrib, "tolinverse" ) ) {
      astClearTolInverse( this );

/* If the attribute is still not recognised, pass it on to the parent
   method for further interpretation. */
   } else {
      (*parent_clearattrib)( this_object, attrib, status );
   }
}

static int Equal( AstObject *this_object, AstObject *that_object, int *status ) {
/*
*  Name:
*     Equal

*  Purpose:
*     Test if two PolyMaps are equivalent.

*  Type:
*     Private function.

*  Synopsis:
*     #include "polymap.h"
*     int Equal( AstObject *this, AstObject *that, int *status )

*  Class Membership:
*     PolyMap member function (over-rides the astEqual protected
*     method inherited from the astMapping class).

*  Description:
*     This function returns a boolean result (0 or 1) to indicate whether
*     two PolyMaps are equivalent.

*  Parameters:
*     this
*        Pointer to the first Object (a PolyMap).
*     that
*        Pointer to the second Object.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     One if the PolyMaps are equivalent, zero otherwise.

*  Notes:
*     - A value of zero will be returned if this function is invoked
*     with the global status set, or if it should fail for any reason.
*/

/* Local Variables: */
   AstPolyMap *that;
   AstPolyMap *this;
   int i, j, k;
   int nin;
   int nout;
   int result;
   int tmp;

/* Initialise. */
   result = 0;

/* Check the global error status. */
   if ( !astOK ) return result;

/* Obtain pointers to the two PolyMap structures. */
   this = (AstPolyMap *) this_object;
   that = (AstPolyMap *) that_object;

/* Check the second object is a PolyMap. We know the first is a
   PolyMap since we have arrived at this implementation of the virtual
   function. */
   if( astIsAPolyMap( that ) ) {

/* Get the number of inputs and outputs and check they are the same for both. */
      nin = astGetNin( this );
      nout = astGetNout( this );
      if( astGetNin( that ) == nin && astGetNout( that ) == nout ) {

/* If the Invert flags for the two PolyMaps differ, it may still be possible
   for them to be equivalent. First compare the PolyMaps if their Invert
   flags are the same. In this case all the attributes of the two PolyMaps
   must be identical. */
         if( astGetInvert( this ) == astGetInvert( that ) ) {

/* Assume they are the same until we find a difference. */
            result = 1;

/* The "_f" and "_i" suffixes in the PolyMap structure array names refer
   to the forward and inverse transformations of the original uninverted
   PolyMap. So we need to ensure the "nout" and "nin" values also refer
   to the uninverted values. So if the PolyMaps are inverted, swap nout
   and nin. */
            if(  astGetInvert( this ) ) {
               tmp = nin;
               nin = nout;
               nout = tmp;
            }

/* Check properties of the forward transformation. */
            if( this->ncoeff_f && that->ncoeff_f ) {
               for( i = 0; i < nout && result; i++ ) {
                  if( this->ncoeff_f[ i ] != that->ncoeff_f[ i ] ){
                     result = 0;
                  }
               }
            } else if( this->ncoeff_f || that->ncoeff_f ) {
               result = 0;
            }

            if( this->mxpow_f && that->mxpow_f ) {
               for( i = 0; i < nout && result; i++ ) {
                  if( this->mxpow_f[ i ] != that->mxpow_f[ i ] ){
                     result = 0;
                  }
               }
            } else if( this->mxpow_f || that->mxpow_f ) {
               result = 0;
            }

            if( this->coeff_f && that->coeff_f ) {
               for( i = 0; i < nout && result; i++ ) {
                  for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) {
                     if( !astEQUAL( this->coeff_f[ i ][ j ],
                                    that->coeff_f[ i ][ j ] ) ) {
                        result = 0;
                     }
                  }
               }
            } else if( this->coeff_f || that->coeff_f ) {
               result = 0;
            }

            if( this->power_f && that->power_f ) {
               for( i = 0; i < nout && result; i++ ) {
                  for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) {
                     for( k = 0; k < nin && result; k++ ) {
                        if( this->power_f[ i ][ j ][ k ] !=
                            that->power_f[ i ][ j ][ k ] ) {
                           result = 0;
                        }
                     }
                  }
               }
            } else if( this->power_f || that->power_f ) {
               result = 0;
            }

/* Check properties of the inverse transformation. */
            if( this->ncoeff_i && that->ncoeff_i ) {
               for( i = 0; i < nout && result; i++ ) {
                  if( this->ncoeff_i[ i ] != that->ncoeff_i[ i ] ){
                     result = 0;
                  }
               }
            } else if( this->ncoeff_i || that->ncoeff_i ) {
               result = 0;
            }

            if( this->mxpow_i && that->mxpow_i ) {
               for( i = 0; i < nout && result; i++ ) {
                  if( this->mxpow_i[ i ] != that->mxpow_i[ i ] ){
                     result = 0;
                  }
               }
            } else if( this->mxpow_i || that->mxpow_i ) {
               result = 0;
            }

            if( this->coeff_f && that->coeff_f ) {
               for( i = 0; i < nout && result; i++ ) {
                  for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) {
                     if( !astEQUAL( this->coeff_f[ i ][ j ],
                                    that->coeff_f[ i ][ j ] ) ) {
                        result = 0;
                     }
                  }
               }
            } else if( this->coeff_f || that->coeff_f ) {
               result = 0;
            }

            if( this->power_f && that->power_f ) {
               for( i = 0; i < nout && result; i++ ) {
                  for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) {
                     for( k = 0; k < nin && result; k++ ) {
                        if( this->power_f[ i ][ j ][ k ] !=
                            that->power_f[ i ][ j ][ k ] ) {
                           result = 0;
                        }
                     }
                  }
               }
            } else if( this->power_f || that->power_f ) {
               result = 0;
            }

/* If the Invert flags for the two PolyMaps differ, the attributes of the two
   PolyMaps must be inversely related to each other. */
         } else {

/* In the specific case of a PolyMap, Invert flags must be equal. */
            result = 0;

         }
      }
   }

/* If an error occurred, clear the result value. */
   if ( !astOK ) result = 0;

/* Return the result, */
   return result;
}

static double *FitPoly1D( AstPolyMap *this, int forward, int nsamp, double acc,
                          int order, double **table, double scales[2], int *ncoeff,
                          double *racc, int *status ){
/*
*  Name:
*     FitPoly1D

*  Purpose:
*     Fit a (1-in,1-out) polynomial to a supplied set of data.

*  Type:
*     Private function.

*  Synopsis:
*     double *FitPoly1D( AstPolyMap *this, int forward, int nsamp, double acc,
*                        int order, double **table, double scales[2], int *ncoeff,
*                        double *racc, int *status )

*  Description:
*     This function fits a least squares 1D polynomial curve to the
*     positions in a supplied table, and returns the coefficients of a
*     PolyMap to describe the fit. For the purposes of this function,
*     the input to the fit is refered to as x1 and the output as y1. So
*     the returned coefficients describe a PolyMap with forward
*     transformation:
*
*     y1 = P1( x1 )

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     forward
*        Non-zero if the forward transformation of "this" is being
*        replaced. Zero if the inverse transformation of "this" is being
*        replaced.
*     nsamp
*        The number of (x1,y1) positions in the supplied table.
*     acc
*        The required accuracy, expressed as a geodesic distance within
*        the polynomials output space (not the normalised tabular space).
*     order
*        The maximum power (plus one) of x1 within P1. So for instance, a
*        value of 3 refers to a quadratic polynomial.
*     table
*        Pointer to an array of 2 pointers. Each of these pointers points
*        to an array of "nsamp" doubles, being the normalised and sampled
*        values for x1 and y1 in that order.
*     scales
*        Array holding the scaling factors used to produced the normalised
*        values in the two columns of the table. Multiplying the normalised
*        table values by the scale factor produces input or output axis
*        values for the returned PolyMap
*     ncoeff
*        Pointer to an int in which to return the number of coefficients
*        described by the returned array.
*     racc
*        Pointer to a double in which to return the achieved accuracy
*        (which may be greater than "acc").
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     A pointer to an array of doubles defining the polynomial in the
*     form required by the PolyMap contructor. The number of coefficients
*     is returned via "ncoeff". If the polynomial could not be found,
*     then NULL is returned. The returned pointer should be freed using
*     astFree when no longer needed.

*/

/* Local Variables: */
   AstMinPackData data;
   double *coeffs;
   double *pc;
   double *pr;
   double *pxp1;
   double *result;
   double *work1;
   double *work2;
   double *work4;
   double f1;
   double f2;
   double maxterm;
   double term;
   double tv;
   int *work3;
   int info;
   int k;
   int ncof;
   int w1;

/* Initialise returned value */
   result = NULL;
   *ncoeff = 0;
   *racc = 10*acc;

/* Check inherited status */
   if( !astOK ) return result;

/* Number of coefficients per poly. */
   ncof = order;

/* Initialise the elements of the structure. */
   data.order = order;
   data.nsamp = nsamp;
   data.init_jac = 1;
   data.xp1 = astMalloc( nsamp*order*sizeof( double ) );
   data.xp2 = NULL;
   data.y[ 0 ] = table[ 1 ];
   data.y[ 1 ] = NULL;

/* Work space to hold coefficients. */
   coeffs = astMalloc( ncof*sizeof( double ) );

/* Other work space. */
   work1 = astMalloc( nsamp*sizeof( double ) );
   work2 = astMalloc( ncof*nsamp*sizeof( double ) );
   work3 = astMalloc( ncof*sizeof( int ) );
   work4 = astMalloc( (5*ncof+nsamp)*sizeof( double ) );
   if( astOK ) {

/* Find all the required powers of x1 and store them in the "xp1"
   component of the data structure. The required initialisation is done
   differently for different subclasses of PolyMap, so we need to wrap
   it up in a virtual function. */
      astFitPoly1DInit( this, forward, table, &data, scales );

/* The initial guess at the coefficient values represents a unit
   transformation in (normalised) tabulated (x,y) values. Using normalised
   values means that we are, effectively, including a guess at the linear
   scaling factor between input and output of the PolyMap (e.g. the PolyMap
   may have inputs in mm and outputs in radians). */
      for( k = 0; k < ncof; k++ ) coeffs[ k ] = 0.0;
      coeffs[ 1 ] = 1.0;

/* Find the best coefficients */
      info = lmder1( MPFunc1D, &data, nsamp, ncof, coeffs, work1, work2, nsamp,
                     sqrt(DBL_EPSILON), work3, work4, (5*ncof+nsamp) );
      if( info == 0 ) astError( AST__MNPCK, "astPolyMap(PolyTran): Minpack error "
                                "detected (possible programming error).", status );

/* Return the achieved accuracy. The "work1" array holds the normalised Y
   residuals at each tabulated point. */
      pr = work1;
      tv = 0.0;
      for( k = 0; k < nsamp; k++,pr++ ) tv += (*pr)*(*pr);
      *racc = scales[ 1 ]*sqrt( tv/nsamp );

/* The best fitting polynomial coefficient found above relate to the
   polynomial between the scaled positions stored in "table". These
   scaled positions are related to PolyMap input/output axis values via
   the scale factors supplied in "scales". Find the initial factor for the
   current output. */
      f1 = scales[ 1 ];
      f2 = 1.0;

/* Look at each coefficient. */
      pc = coeffs;
      for( w1 = 0; w1 < order; w1++,pc++ ) {

/* Get a pointer to the powers of X1 appropriate for the current coefficient,
    at the first sample. */
         pxp1 = data.xp1 + w1;

/* We find the contribution which this coefficient makes to the total
   polynomial value. Find the maximum contribution made at any sample
   points. */
         maxterm = 0.0;
         for( k = 0; k < nsamp; k++ ) {

/* Get the absolute value of the polynomial term that uses the current
   coefficient. */
            term = fabs( ( *pc )*( *pxp1 ) );

/* Update the maximum term found at any sample. */
            if( term > maxterm ) maxterm = term;

/* Increment the pointers to refer to the next sample. */
            pxp1 += order;
         }

/* If the maximum contribution made by this term is less than the
   required accuracy, set the coefficient value to zero. */
         if( maxterm*f1 < acc ) {
            *pc = 0.0;

/* Scale the best fitting polynomial coefficient found above to take
   account of the fact that the tabulated input and output positions in
   "table" were are not actual PolyMap input and output axis values, but
   are scaled by the factors stored in "scales". */
         } else {
            *pc *= f1/f2;
         }

         f2 *= scales[ 0 ];
      }

/* Convert the array of coefficients into PolyMap form. */
      result = astMalloc( ncof*3*sizeof( double ) );

      *ncoeff = 0;
      pr = result;
      pc = coeffs;
      for( w1 = 0; w1 < order; w1++,pc++ ) {
         if( *pc != 0.0 ) {
            *(pr++) = *pc;
            *(pr++) = 1;
            *(pr++) = w1;
            (*ncoeff)++;
         }
      }

/* Truncate the returned array. */
      result = astRealloc( result, (*ncoeff)*3*sizeof( double ) );
   }

/* Free resources. */
   coeffs = astFree( coeffs );
   data.xp1 = astFree( data.xp1 );
   work1 = astFree( work1 );
   work2 = astFree( work2 );
   work3 = astFree( work3 );
   work4 = astFree( work4 );

/* Return the coefficient array. */
   return result;

}

static void FitPoly1DInit( AstPolyMap *this, int forward, double **table,
                           AstMinPackData *data, double *scales, int *status ){
/*
*+
*  Name:
*     astFitPoly1DInit

*  Purpose:
*     Perform initialisation needed for FitPoly1D

*  Type:
*     Protected function.

*  Synopsis:
*     #include "polymap.h"
*     void astFitPoly1DInit( AstPolyMap *this, int forward, double **table,
*                            AstMinPackData *data, double *scales,
*                            int *status )

*  Class Membership:
*     PolyMap virtual function.

*  Description:
*     This function performs initialisation needed for FitPoly1D.

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     forward
*        Non-zero if the forward transformation of "this" is being
*        replaced. Zero if the inverse transformation of "this" is being
*        replaced.
*     table
*        Pointer to an array of 2 pointers. Each of these pointers points
*        to an array of "nsamp" doubles, being the scaled and sampled values
*        for x1 and y1 in that order.
*     data
*        Pointer to a structure holding information to pass the the
*        service function invoked by the minimisation function.
*     scales
*        Array holding the scaling factors for the two columns of the table.
*        Multiplying the table values by the scale factor produces PolyMap
*        input or output axis values.
*-
*/

/* Local Variables; */
   double *px1;
   double *pxp1;
   double tv;
   double x1;
   int k;
   int w1;

/* Check the local error status. */
   if ( !astOK ) return;

/* Get pointers to the supplied x1 values. */
   px1 = table[ 0 ];

/* Get pointers to the location for the next power of x1. */
   pxp1 = data->xp1;

/* Loop round all samples. */
   for( k = 0; k < data->nsamp; k++ ) {

/* Get the current x1 value. */
      x1 = *(px1++);

/* Find all the required powers of x1 and store them in the "xp1"
   component of the data structure. */
      tv = 1.0;
      for( w1 = 0; w1 < data->order; w1++ ) {
         *(pxp1++) = tv;
         tv *= x1;
      }
   }
}

static double *FitPoly2D( AstPolyMap *this, int forward, int nsamp, double acc,
                          int order, double **table, double scales[4],
                          int *ncoeff, double *racc, int *status ){
/*
*  Name:
*     FitPoly2D

*  Purpose:
*     Fit a (2-in,2-out) polynomial to a supplied set of data.

*  Type:
*     Private function.

*  Synopsis:
*     double *FitPoly2D( AstPolyMap *this, int forward, int nsamp, double acc,
*                        int order, double **table, double scales[4],
*                        int *ncoeff, double *racc, int *status )

*  Description:
*     This function fits a pair of least squares 2D polynomial surfaces
*     to the positions in a supplied table, and returns the coefficients
*     of a PolyMap to describe the fit. For the purposes of this function,
*     the inputs to the fit is refered to as (x1,x2) and the output as
*     (y1,y2). So the returned coefficients describe a PolyMap with forward
*     transformations:
*
*     y1 = P1( x1, x2 )
*     y2 = P2( x1, x2 )
*
*     P1 and P2 have the same maximum power on each input (specified by
*     the "order" parameter).

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     forward
*        Non-zero if the forward transformation of "this" is being
*        replaced. Zero if the inverse transformation of "this" is being
*        replaced.
*     nsamp
*        The number of (x1,x2,y1,y2) positions in the supplied table.
*     acc
*        The required accuracy, expressed as a geodesic distance within
*        the polynomials output space (not the normalised tabular space).
*     order
*        The maximum power (plus one) of x1 or x2 within P1 and P2. So for
*        instance, a value of 3 refers to a quadratic polynomial. Note, cross
*        terms with total powers greater than or equal to "order" are not
*        included in the fit. So the maximum number of terms in the fitted
*        polynomial is order*(order+1)/2.
*     table
*        Pointer to an array of 4 pointers. Each of these pointers points
*        to an array of "nsamp" doubles, being the normalised and sampled
*        values for x1, x2, y1 or y2 in that order.
*     scales
*        Array holding the scaling factors used to produced the normalised
*        values in the four columns of the table. Multiplying the normalised
*        table values by the scale factor produces input or output axis
*        values for the returned PolyMap
*     ncoeff
*        Pointer to an ant in which to return the number of coefficients
*        described by the returned array.
*     racc
*        Pointer to a double in which to return the achieved accuracy
*        (which may be greater than "acc").
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     A pointer to an array of doubles defining the polynomial in the
*     form required by the PolyMap contructor. The number of coefficients
*     is returned via "ncoeff". If the polynomial could not be found with
*     sufficient accuracy , then NULL is returned. The returned pointer
*     should be freed using astFree when no longer needed.

*/

/* Local Variables: */
   AstMinPackData data;
   double *coeffs;
   double *pc;
   double *pr;
   double *pxp1;
   double *pxp2;
   double *result;
   double *work1;
   double *work2;
   double *work4;
   double f1;
   double f20;
   double f2;
   double f3;
   double facc;
   double maxterm;
   double term;
   double tv;
   int *work3;
   int info;
   int iout;
   int k;
   int ncof;
   int w12;
   int w1;
   int w2;

/* Initialise returned value */
   result = NULL;
   *ncoeff = 0;
   *racc = 10*acc;

/* Check inherited status */
   if( !astOK ) return result;

/* Number of coefficients per poly. */
   ncof = order*( order + 1 )/2;

/* Initialise the elements of the structure. */
   data.order = order;
   data.nsamp = nsamp;
   data.init_jac = 1;
   data.xp1 = astMalloc( nsamp*order*sizeof( double ) );
   data.xp2 = astMalloc( nsamp*order*sizeof( double ) );
   data.y[ 0 ] = table[ 2 ];
   data.y[ 1 ] = table[ 3 ];

/* Work space to hold coefficients. */
   coeffs = astMalloc( 2*ncof*sizeof( double ) );

/* Other work space. */
   work1 = astMalloc( 2*nsamp*sizeof( double ) );
   work2 = astMalloc( 4*ncof*nsamp*sizeof( double ) );
   work3 = astMalloc( 2*ncof*sizeof( int ) );
   work4 = astMalloc( 2*(5*ncof+nsamp)*sizeof( double ) );
   if( astOK ) {

/* Find all the required powers of x1 and x2 and store them in the "xp1"
   and "xp2" components of the data structure. The required initialisation
   is done differently for different subclasses of PolyMap, so we need to
   wrap it up in a virtual function. */
      astFitPoly2DInit( this, forward, table, &data, scales );

/* The initial guess at the coefficient values represents a unit
   transformation in (normalised) tabulated (x,y) values. Using normalised
   values means that we are, effectively, including a guess at the linear
   scaling factor between input and output of the PolyMap (e.g. the PolyMap
   may have inputs in mm and outputs in radians). */
      for( k = 0; k < 2*ncof; k++ ) coeffs[ k ] = 0.0;
      coeffs[ 1 ] = 1.0;
      coeffs[ 5 ] = 1.0;

/* Find the best coefficients */
      info = lmder1( MPFunc2D, &data, 2*nsamp, 2*ncof, coeffs, work1, work2,
                     2*nsamp, sqrt(DBL_EPSILON), work3, work4, 2*(5*ncof+nsamp) );
      if( info == 0 ) astError( AST__MNPCK, "astPolyMap(PolyTran): Minpack error "
                                "detected (possible programming error).", status );

/* Return the achieved accuracy. */
      pr = work1;
      tv = 0.0;
      for( k = 0; k < 2*nsamp; k++,pr++ ) tv += (*pr)*(*pr);
      facc = 1.0/(scales[2]*scales[2]) + 1.0/(scales[3]*scales[3]);
      *racc = sqrt( tv/(2*nsamp*facc) );

/* Pointer to the first coefficient. */
      pc = coeffs;

/* Look at coefficients for each output in turn. */
      for( iout = 0; iout < 2 && astOK; iout++ ) {

/* The best fitting polynomial coefficient found above relate to the
   polynomial between the scaled positions stored in "table". These
   scaled positions are related to PolyMap input/output axis values via
   the scale factors supplied in "scales". Find the initial factor for the
   current output. */
         f1 = scales[ 2 + iout ];

/* Look at each coefficient for the current output. */
         f20 = 1.0;
         for( w12 = 0; w12 < order; w12++ ) {
            f3 = 1.0;
            f2 = f20;
            for( w2 = 0; w2 <= w12; w2++,pc++ ) {
               w1 = w12 - w2;

/* Get pointers to the powers of X1 and X2 appropriate for the current
   coefficient, at the first sample. */
               pxp1 = data.xp1 + w1;
               pxp2 = data.xp2 + w2;

/* We find the contribution which this coefficient makes to the total
   polynomial value. Find the maximum contribution made at any sample
   points. */
               maxterm = 0.0;
               for( k = 0; k < nsamp; k++ ) {

/* Get the absolute value of the polynomial term that uses the current
   coefficient. */
                  term = fabs( ( *pc )*( *pxp1 )*( *pxp2 ) );

/* Update the maximum term found at any sample. */
                  if( term > maxterm ) maxterm = term;

/* Increment the pointers to refer to the next sample. */
                  pxp1 += order;
                  pxp2 += order;
               }

/* If the maximum contribution made by this term is less than the
   required accuracy, set the coefficient value to zero. */
               if( maxterm*f1 < acc ) {
                  *pc = 0.0;

/* Scale the best fitting polynomial coefficient found above to take
   account of the fact that the tabulated input and output positions in
   "table" were are not actual PolyMap input and output axis values, but
   are scaled by the factors stored in "scales". */
               } else {
                  *pc *= f1/( f2*f3 );
               }

               f2 /= scales[ 0 ];
               f3 *= scales[ 1 ];
            }

            f20 *= scales[ 0 ];
         }
      }

/* Convert the array of coefficients into PolyMap form. */
      result = astMalloc( 2*ncof*4*sizeof( double ) );

      *ncoeff = 0;
      pr = result;
      pc = coeffs;
      for( iout = 0; iout < 2 && astOK; iout++ ) {
         for( w12 = 0; w12 < order; w12++ ) {
            for( w2 = 0; w2 <= w12; w2++,pc++ ) {
               w1 = w12 - w2;
               if( *pc != 0.0 ) {
                  *(pr++) = *pc;
                  *(pr++) = iout + 1;
                  *(pr++) = w1;
                  *(pr++) = w2;
                  (*ncoeff)++;
               }
            }
         }
      }

/* Truncate the returned array. */
      result = astRealloc( result, (*ncoeff)*4*sizeof( double ) );
   }

/* Free resources. */
   coeffs = astFree( coeffs );
   data.xp1 = astFree( data.xp1 );
   data.xp2 = astFree( data.xp2 );
   work1 = astFree( work1 );
   work2 = astFree( work2 );
   work3 = astFree( work3 );
   work4 = astFree( work4 );

/* Return the coefficient array. */
   return result;

}

static void FitPoly2DInit( AstPolyMap *this, int forward, double **table,
                           AstMinPackData *data, double *scales, int *status ){
/*
*+
*  Name:
*     astFitPoly2DInit

*  Purpose:
*     Perform initialisation needed for FitPoly2D

*  Type:
*     Protected function.

*  Synopsis:
*     #include "polymap.h"
*     void astFitPoly2DInit( AstPolyMap *this, int forward, double **table,
*                            AstMinPackData *data, double *scales,
*                            int *status )

*  Class Membership:
*     PolyMap virtual function.

*  Description:
*     This function performs initialisation needed for FitPoly2D.

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     forward
*        Non-zero if the forward transformation of "this" is being
*        replaced. Zero if the inverse transformation of "this" is being
*        replaced.
*     table
*        Pointer to an array of 4 pointers. Each of these pointers points
*        to an array of "nsamp" doubles, being the scaled and sampled values
*        for x1, x2, y1 or y2 in that order.
*     data
*        Pointer to a structure holding information to pass the the
*        service function invoked by the minimisation function.
*     scales
*        Array holding the scaling factors for the four columns of the table.
*        Multiplying the table values by the scale factor produces PolyMap
*        input or output axis values.
*-
*/

/* Local Variables; */
   double *px1;
   double *px2;
   double *pxp1;
   double *pxp2;
   double tv;
   double x1;
   double x2;
   int k;
   int w1;
   int w2;

/* Check the local error status. */
   if ( !astOK ) return;

/* Get pointers to the supplied x1 and x2 values. */
   px1 = table[ 0 ];
   px2 = table[ 1 ];

/* Get pointers to the location for the next power of x1 and x2. */
   pxp1 = data->xp1;
   pxp2 = data->xp2;

/* Loop round all samples. */
   for( k = 0; k < data->nsamp; k++ ) {

/* Get the current x1 and x2 values. */
      x1 = *(px1++);
      x2 = *(px2++);

/* Find all the required powers of x1 and store them in the "xp1"
   component of the data structure. */
      tv = 1.0;
      for( w1 = 0; w1 < data->order; w1++ ) {
         *(pxp1++) = tv;
         tv *= x1;
      }

/* Find all the required powers of x2 and store them in the "xp2"
   comonent of the data structure. */
      tv = 1.0;
      for( w2 = 0; w2 < data->order; w2++ ) {
         *(pxp2++) = tv;
         tv *= x2;
      }
   }
}

static void FreeArrays( AstPolyMap *this, int forward, int *status ) {
/*
*  Name:
*     FreeArrays

*  Purpose:
*     Free the dynamic arrays contained within a PolyMap

*  Type:
*     Private function.

*  Synopsis:
*     void FreeArrays( AstPolyMap *this, int forward, int *status )

*  Description:
*     This function frees all the dynamic arrays allocated as part of a
*     PolyMap.

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     forward
*        If non-zero, the arrays for the forward transformation are freed.
*        Otherwise, the arrays for the inverse transformation are freed.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     void

*  Notes:
*     This function attempts to execute even if the global error status is
*     set.
*/

/* Local Variables: */
   int nin;                     /* Number of inputs */
   int nout;                    /* Number of outputs */
   int i;                       /* Loop count */
   int j;                       /* Loop count */

/* Get the number of inputs and outputs of the uninverted Mapping. */
   nin = ( (AstMapping *) this )->nin;
   nout = ( (AstMapping *) this )->nout;

/* Free the dynamic arrays for the forward transformation. */
   if( forward != astGetInvert( this ) ) {

      if( this->coeff_f ) {
         for( i = 0; i < nout; i++ ) {
            this->coeff_f[ i ] = astFree( this->coeff_f[ i ] );
         }
         this->coeff_f = astFree( this->coeff_f );
      }

      if( this->power_f ) {
         for( i = 0; i < nout; i++ ) {
            if( this->ncoeff_f && this->power_f[ i ] ) {
               for( j = 0; j < this->ncoeff_f[ i ]; j++ ) {
                  this->power_f[ i ][ j ] = astFree( this->power_f[ i ][ j ] );
               }
            }
            this->power_f[ i ] = astFree( this->power_f[ i ] );
         }
         this->power_f = astFree( this->power_f );
      }

      this->ncoeff_f = astFree( this->ncoeff_f );
      this->mxpow_f = astFree( this->mxpow_f );

/* Free the dynamic arrays for the inverse transformation. */
   } else {

      if( this->coeff_i ) {
         for( i = 0; i < nin; i++ ) {
            this->coeff_i[ i ] = astFree( this->coeff_i[ i ] );
         }
         this->coeff_i = astFree( this->coeff_i );
      }

      if(this->power_i ) {
         for( i = 0; i < nin; i++ ) {
            if( this->ncoeff_i && this->power_i[ i ] ) {
               for( j = 0; j < this->ncoeff_i[ i ]; j++ ) {
                  this->power_i[ i ][ j ] = astFree( this->power_i[ i ][ j ] );
               }
            }
            this->power_i[ i ] = astFree( this->power_i[ i ] );
         }
         this->power_i = astFree( this->power_i );
      }

      this->ncoeff_i = astFree( this->ncoeff_i );
      this->mxpow_i = astFree( this->mxpow_i );
   }
}

static const char *GetAttrib( AstObject *this_object, const char *attrib, int *status ) {
/*
*  Name:
*     GetAttrib

*  Purpose:
*     Get the value of a specified attribute for a PolyMap.

*  Type:
*     Private function.

*  Synopsis:
*     #include "polymap.h"
*     const char *GetAttrib( AstObject *this, const char *attrib, int *status )

*  Class Membership:
*     PolyMap member function (over-rides the protected astGetAttrib
*     method inherited from the Mapping class).

*  Description:
*     This function returns a pointer to the value of a specified
*     attribute for a PolyMap, formatted as a character string.

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     attrib
*        Pointer to a null-terminated string containing the name of
*        the attribute whose value is required. This name should be in
*        lower case, with all white space removed.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     - Pointer to a null-terminated string containing the attribute
*     value.

*  Notes:
*     - The returned string pointer may point at memory allocated
*     within the PolyMap, or at static memory. The contents of the
*     string may be over-written or the pointer may become invalid
*     following a further invocation of the same function or any
*     modification of the PolyMap. A copy of the string should
*     therefore be made if necessary.
*     - A NULL pointer will be returned if this function is invoked
*     with the global error status set, or if it should fail for any
*     reason.
*/

/* Local Variables: */
   astDECLARE_GLOBALS           /* Pointer to thread-specific global data */
   AstPolyMap *this;            /* Pointer to the PolyMap structure */
   const char *result;          /* Pointer value to return */
   double dval;                 /* Floating point attribute value */
   int ival;                    /* Integer attribute value */

/* Initialise. */
   result = NULL;

/* Check the global error status. */
   if ( !astOK ) return result;

/* Get a pointer to the thread specific global data structure. */
   astGET_GLOBALS(this_object);

/* Obtain a pointer to the PolyMap structure. */
   this = (AstPolyMap *) this_object;

/* Compare "attrib" with each recognised attribute name in turn,
   obtaining the value of the required attribute. If necessary, write
   the value into "getattrib_buff" as a null-terminated string in an appropriate
   format.  Set "result" to point at the result string. */

/* IterInverse. */
/* ------------ */
   if ( !strcmp( attrib, "iterinverse" ) ) {
      ival = astGetIterInverse( this );
      if ( astOK ) {
         (void) sprintf( getattrib_buff, "%d", ival );
         result = getattrib_buff;
      }

/* NiterInverse. */
/* ------------- */
   } else if ( !strcmp( attrib, "niterinverse" ) ) {
      ival = astGetNiterInverse( this );
      if ( astOK ) {
         (void) sprintf( getattrib_buff, "%d", ival );
         result = getattrib_buff;
      }

/* TolInverse. */
/* ----------- */
   } else if ( !strcmp( attrib, "tolinverse" ) ) {
      dval = astGetTolInverse( this );
      if ( astOK ) {
         (void) sprintf( getattrib_buff, "%.*g", AST__DBL_DIG, dval );
         result = getattrib_buff;
      }

/* If the attribute name was not recognised, pass it on to the parent
   method for further interpretation. */
   } else {
      result = (*parent_getattrib)( this_object, attrib, status );
   }

/* Return the result. */
   return result;

}

static AstPolyMap **GetJacobian( AstPolyMap *this, int *status ){
/*
*  Name:
*     GetJacobian

*  Purpose:
*     Get a description of a Jacobian matrix for the original
*     fwd transformation of a PolyMap.

*  Type:
*     Private function.

*  Synopsis:
*     AstPolyMap **GetJacobian( AstPolyMap *this, int *status )

*  Description:
*     This function returns a set of PolyMaps which define the Jacobian
*     matrix of the original forward transformation of the supplied PolyMap
*     (i.e. the Negated attribute is assumed to be zero).
*
*     The Jacobian matrix has "nout" rows and "nin" columns, where "nin"
*     and "nout" are the number of inputs and outputs of the original PolyMap.
*     Row "i", column "j" of the matrix holds the rate of change of the
*     i'th PolyMap output with respect to the j'th PolyMap input.
*
*     Since the values in the Jacobian matrix vary across the input space
*     of the PolyMap, the matrix is returned in the form of a set of new
*     PolyMaps which generate the elements of the Jacobian for any given
*     position in the input space. The "nout" values in a single column of
*     the Jacobian matrix are generated by the "nout" outputs from a single
*     new PolyMap. The whole matrix is described by "nin" PolyMaps.
*
*     The returned PolyMaps are cached in the supplied PolyMap object in
*     order to speed up subsequent calls to this function.

*  Parameters:
*     this
*        The PolyMap for which the Jacbian is required.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     A pointer to an array of "nin" PolyMap pointers, where "nin" is the number
*     of inputs for the sipplied PolyMap. The returned array should not be changed
*     in any way, and the PolyMaps should not be freed (they will be freed when
*     the supplied PolyMap is deleted).

*/

/* Local Variables: */
   double *coeffs;
   double *pc;
   int icof;
   int icol;
   int iin;
   int irow;
   int ncof;
   int ncof_row;
   int ncof_total;
   int nin;
   int nout;
   int power;

/* Check inherited status */
   if( !astOK ) return NULL;

/* Ensure there is a Jacobian stored in the PolyMap. */
   if( !this->jacobian ) {

/* Get the number of inputs and outputs of the original forward
   transformation. */
      if( astGetInvert( this ) ){
         nout = astGetNin( this );
         nin = astGetNout( this );
      } else {
         nin = astGetNin( this );
         nout = astGetNout( this );
      }

/* Allocate memory to hold pointers to the PolyMaps used to describe the
   Jacobian matrix. */
      this->jacobian = astCalloc( nin, sizeof(AstPolyMap *) );

/* Find the total number of coefficients used to describe the supplied
   PolyMap (forward transformation) and allocate work space to hold the
   coefficients for a single new PolyMap forward transformation. */
      ncof = 0;
      for( irow = 0; irow <  nout; irow++ ) {
         ncof += this->ncoeff_f[ irow ];
      }
      coeffs = astMalloc( ncof*( 2 + nin )*sizeof( double ) );

/* Check pointers can be used safely. */
      if( astOK ) {

/* The Jacobian matrix has "nout" rows and "nin" columns. The "nout" values
   in a single column of the Jacobian matrix corresponds to the "nout" outputs
   from a single PolyMap. The whole matrix is described by "nin" PolyMaps.
   Loop over each column of the Matrix, creating the corresponding PolyMap
   for each. */
         for( icol = 0; icol <  nin; icol++ ) {

/* Initialise the total number of coefficients used to describe the
   element of the PolyMap. */
            ncof_total = 0;

/* Loop over each row of the Jacobian matrix (i.e. each PolyMap output). */
            pc = coeffs;
            for( irow = 0; irow <  nout; irow++ ) {

/* Loop over each coefficient used in the polynomial that generates the
   current PolyMap output. */
               ncof_row = this->ncoeff_f[ irow ];
               for( icof = 0; icof <  ncof_row; icof++ ) {

/* Get the power of input "icol" associated with the current coefficient. */
                  power = (int)( this->power_f[ irow ][ icof ][ icol ] + 0.5 );

/* We can skip the coefficient if the power is zero. */
                  if( power > 0 ) {
                     ncof_total++;

/* Store the coefficient value, modified so that it describes a
   polynomial that has been differentiated with respect to input "icol". */
                     *(pc++) = this->coeff_f[ irow ][ icof ]*power;

/* Store the output PolyMap to which the coeff relates. */
                     *(pc++) = irow + 1;

/* Store the powers of the inputs associated with the coeff. These are
   the same as the original powers, except that the power of "icol"
   (the input with respect to which the output has been differentiated)
   is reduced by one. */
                     for( iin = 0; iin <  nin; iin++ ) {
                        if( iin != icol ) {
                           *(pc++) = this->power_f[ irow ][ icof ][ iin ];
                        } else {
                           *(pc++) = this->power_f[ irow ][ icof ][ iin ] - 1;
                        }
                     }
                  }
               }
            }

/* Create the PolyMap and store a pointer to it in the jacobian array in
   the supplied PolyMap. */
            (this->jacobian)[ icol ] = astPolyMap( nin, nout, ncof_total, coeffs,
                                                   0, NULL, " ", status );
         }
      }

/* Free resources */
      coeffs = astFree( coeffs );
   }

/* Return the Jacobian. */
   return this->jacobian;
}

static int GetObjSize( AstObject *this_object, int *status ) {
/*
*  Name:
*     GetObjSize

*  Purpose:
*     Return the in-memory size of an Object.

*  Type:
*     Private function.

*  Synopsis:
*     #include "polymap.h"
*     int GetObjSize( AstObject *this, int *status )

*  Class Membership:
*     PolyMap member function (over-rides the astGetObjSize protected
*     method inherited from the parent class).

*  Description:
*     This function returns the in-memory size of the supplied PolyMap,
*     in bytes.

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     The Object size, in bytes.

*  Notes:
*     - A value of zero will be returned if this function is invoked
*     with the global status set, or if it should fail for any reason.
*/

/* Local Variables: */
   AstPolyMap *this;
   int ic;
   int nc;
   int result;

/* Initialise. */
   result = 0;

/* Check the global error status. */
   if ( !astOK ) return result;

/* Obtain a pointers to the PolyMap structure. */
   this = (AstPolyMap *) this_object;

/* Invoke the GetObjSize method inherited from the parent class, and then
   add on any components of the class structure defined by this class
   which are stored in dynamically allocated memory. */
   result = (*parent_getobjsize)( this_object, status );

   if( this->jacobian ) {
      nc = astGetNin( this );
      for( ic = 0; ic < nc; ic++ ) {
         result +=  astGetObjSize( (this->jacobian)[ ic ] );
      }
      result += sizeof( AstPolyMap * )*nc;
   }

/* If an error occurred, clear the result value. */
   if ( !astOK ) result = 0;

/* Return the result, */
   return result;
}

static int GetTranForward( AstMapping *this, int *status ) {
/*
*
*  Name:
*     GetTranForward

*  Purpose:
*     Determine if a PolyMap defines a forward coordinate transformation.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mapping.h"
*     int GetTranForward( AstMapping *this, int *status )

*  Class Membership:
*     PolyMap member function (over-rides the astGetTranForward method
*     inherited from the Mapping class).

*  Description:
*     This function returns a value indicating if the PolyMap is able
*     to perform a forward coordinate transformation.

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     Zero if the forward coordinate transformation is not defined, or 1 if it
*     is.

*  Notes:
*     -  A value of zero will be returned if this function is invoked with the
*     global error status set, or if it should fail for any reason.
*/

/* Local Variables: */
   AstPolyMap *map;            /* Pointer to PolyMap to be queried */
   int result;                 /* The returned value */

/* Check the global error status. */
   if ( !astOK ) return 0;

/* Obtain a pointer to the PolyMap. */
   map = (AstPolyMap *) this;

/* First deal with cases where the PolyMap has not been inverted. */
   if( !astGetInvert( this ) ) {

/* The PolyMap has a defined forward transformation if one or more
   coefficients values were supplied for the original forward
   transformation. It is not possible to replace the original forward
   transformation with an iterative algorithm. */
      result = map->ncoeff_f ? 1 : 0;

/* Now deal with cases where the PolyMap has been inverted. */
   } else {

/* The PolyMap has a defined forward transformation if one or more
   coefficients values were supplied for the original inverse
   transformation, or if the original inverse transformation is being
   approximated using an iterative algorithm. */
      result = ( map->ncoeff_i || astGetIterInverse( map ) ) ? 1 : 0;
   }

/* Return the result. */
   return result;
}

static int GetTranInverse( AstMapping *this, int *status ) {
/*
*
*  Name:
*     GetTranInverse

*  Purpose:
*     Determine if a PolyMap defines an inverse coordinate transformation.

*  Type:
*     Private function.

*  Synopsis:
*     #include "matrixmap.h"
*     int GetTranInverse( AstMapping *this, int *status )

*  Class Membership:
*     PolyMap member function (over-rides the astGetTranInverse method
*     inherited from the Mapping class).

*  Description:
*     This function returns a value indicating if the PolyMap is able
*     to perform an inverse coordinate transformation.

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     Zero if the inverse coordinate transformation is not defined, or 1 if it
*     is.

*  Notes:
*     -  A value of zero will be returned if this function is invoked with the
*     global error status set, or if it should fail for any reason.
*/

/* Local Variables: */
   AstPolyMap *map;            /* Pointer to PolyMap to be queried */
   int result;                 /* The returned value */

/* Check the global error status. */
   if ( !astOK ) return 0;

/* Obtain a pointer to the PolyMap. */
   map = (AstPolyMap *) this;

/* First deal with cases where the PolyMap has not been inverted. */
   if( !astGetInvert( this ) ) {

/* The PolyMap has a defined inverse transformation if one or more
   coefficients values were supplied for the original inverse
   transformation, or if the original inverse transformation is being
   approximated using an iterative algorithm. */
      result = ( map->ncoeff_i || astGetIterInverse( map ) ) ? 1 : 0;

/* Now deal with cases where the PolyMap has been inverted. */
   } else {

/* The PolyMap has a defined inverse transformation if one or more
   coefficients values were supplied for the original forward
   transformation. It is not possible to replace the original forward
   transformation with an iterative algorithm. */
      result = map->ncoeff_f ? 1 : 0;
   }

/* Return the result. */
   return result;
}

void astInitPolyMapVtab_(  AstPolyMapVtab *vtab, const char *name, int *status ) {
/*
*+
*  Name:
*     astInitPolyMapVtab

*  Purpose:
*     Initialise a virtual function table for a PolyMap.

*  Type:
*     Protected function.

*  Synopsis:
*     #include "polymap.h"
*     void astInitPolyMapVtab( AstPolyMapVtab *vtab, const char *name )

*  Class Membership:
*     PolyMap vtab initialiser.

*  Description:
*     This function initialises the component of a virtual function
*     table which is used by the PolyMap class.

*  Parameters:
*     vtab
*        Pointer to the virtual function table. The components used by
*        all ancestral classes will be initialised if they have not already
*        been initialised.
*     name
*        Pointer to a constant null-terminated character string which contains
*        the name of the class to which the virtual function table belongs (it
*        is this pointer value that will subsequently be returned by the Object
*        astClass function).
*-
*/

/* Local Variables: */
   astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
   AstObjectVtab *object;        /* Pointer to Object component of Vtab */
   AstMappingVtab *mapping;      /* Pointer to Mapping component of Vtab */

/* Check the local error status. */
   if ( !astOK ) return;


/* Get a pointer to the thread specific global data structure. */
   astGET_GLOBALS(NULL);

/* Initialize the component of the virtual function table used by the
   parent class. */
   astInitMappingVtab( (AstMappingVtab *) vtab, name );

/* Store a unique "magic" value in the virtual function table. This
   will be used (by astIsAPolyMap) to determine if an object belongs
   to this class.  We can conveniently use the address of the (static)
   class_check variable to generate this unique value. */
   vtab->id.check = &class_check;
   vtab->id.parent = &(((AstMappingVtab *) vtab)->id);

/* Initialise member function pointers. */
/* ------------------------------------ */
/* Store pointers to the member functions (implemented here) that provide
   virtual methods for this class. */
   vtab->PolyPowers = PolyPowers;
   vtab->FitPoly1DInit = FitPoly1DInit;
   vtab->FitPoly2DInit = FitPoly2DInit;
   vtab->PolyTran = PolyTran;
   vtab->PolyCoeffs = PolyCoeffs;

   vtab->ClearIterInverse = ClearIterInverse;
   vtab->GetIterInverse = GetIterInverse;
   vtab->SetIterInverse = SetIterInverse;
   vtab->TestIterInverse = TestIterInverse;

   vtab->ClearNiterInverse = ClearNiterInverse;
   vtab->GetNiterInverse = GetNiterInverse;
   vtab->SetNiterInverse = SetNiterInverse;
   vtab->TestNiterInverse = TestNiterInverse;

   vtab->ClearTolInverse = ClearTolInverse;
   vtab->GetTolInverse = GetTolInverse;
   vtab->SetTolInverse = SetTolInverse;
   vtab->TestTolInverse = TestTolInverse;

/* Save the inherited pointers to methods that will be extended, and
   replace them with pointers to the new member functions. */
   object = (AstObjectVtab *) vtab;
   mapping = (AstMappingVtab *) vtab;

   parent_getobjsize = object->GetObjSize;
   object->GetObjSize = GetObjSize;

#if defined(THREAD_SAFE)
   parent_managelock = object->ManageLock;
   object->ManageLock = ManageLock;
#endif

   parent_clearattrib = object->ClearAttrib;
   object->ClearAttrib = ClearAttrib;
   parent_getattrib = object->GetAttrib;
   object->GetAttrib = GetAttrib;
   parent_setattrib = object->SetAttrib;
   object->SetAttrib = SetAttrib;
   parent_testattrib = object->TestAttrib;
   object->TestAttrib = TestAttrib;

   parent_transform = mapping->Transform;
   mapping->Transform = Transform;
   mapping->GetTranForward = GetTranForward;
   mapping->GetTranInverse = GetTranInverse;

/* Store replacement pointers for methods which will be over-ridden by
   new member functions implemented here. */
   object->Equal = Equal;
   mapping->MapMerge = MapMerge;

/* Declare the destructor and copy constructor. */
   astSetDelete( (AstObjectVtab *) vtab, Delete );
   astSetCopy( (AstObjectVtab *) vtab, Copy );

/* Declare the class dump function. */
   astSetDump( vtab, Dump, "PolyMap", "Polynomial transformation" );

/* If we have just initialised the vtab for the current class, indicate
   that the vtab is now initialised, and store a pointer to the class
   identifier in the base "object" level of the vtab. */
   if( vtab == &class_vtab ) {
      class_init = 1;
      astSetVtabClassIdentifier( vtab, &(vtab->id) );
   }
}

static void IterInverse( AstPolyMap *this, AstPointSet *out, AstPointSet *result,
                         int *status ){
/*
*  Name:
*     IterInverse

*  Purpose:
*     Use an iterative method to evaluate the original inverse transformation
*     of a PolyMap at a set of (original) output positions.

*  Type:
*     Private function.

*  Synopsis:
*     void IterInverse( AstPolyMap *this, AstPointSet *out, AstPointSet *result,
*                       int *status )

*  Description:
*     This function transforms a set of PolyMap positions using the original
*     inverse transformation of the PolyMap (i.e. the Negated attribute
*     is assumed to be zero). An iterative Newton-Raphson method is used
*     which only required the original forward transformation of the PolyMap
*     to be defined.

*  Parameters:
*     this
*        The PolyMap.
*     out
*        A PointSet holding the positions that are to be transformed using
*        the original inverse transformation. These corresponds to
*        outputs of the original (i.e. uninverted) PolyMap
*     result
*        A PointSet into which the transformed positions are to be stored.
*        These corresponds to inputs of the original (i.e. uninverted) PolyMap

*     status
*        Pointer to the inherited status variable.

*/

/* Local Variables: */
   AstMapping *lintrunc;
   AstPointSet *work;
   AstPointSet **ps_jac;
   AstPolyMap **jacob;
   double *vec;
   double *pb;
   double **ptr_work;
   double ***ptr_jac;
   double *mat;
   double **ptr_out;
   double **ptr_in;
   double *pa;
   double det;
   double maxerr;
   double vlensq;
   double xlensq;
   double xx;
   int *flags;
   int *iw;
   int fwd;
   int icol;
   int icoord;
   int ipoint;
   int irow;
   int iter;
   int maxiter;
   int nconv;
   int ncoord;
   int npoint;
   int sing;

/* Check inherited status */
   if( !astOK ) return;

/* Check the PolyMap has equal numbers of inputs and outputs. */
   ncoord = astGetNin( this );
   if( ncoord != astGetNout( this ) ) {
      astError( AST__INTER, "astTransform(%s): Supplied %s has unequal numbers"
                " of inputs and outputs and therefore an iterative inverse "
                "cannot be used (internal AST Programming errpr).", status,
                astGetClass(this), astGetClass(this) );
   }

/* Get information about the Jacobian matrix for the forward polynomial
   transformation. This matrix is a ncoord X ncoord matrix, in which
   element (row=I,col=J) is the rate of change of output coord I with
   respect to input coord J, of the supplied PolyMap's forward transformation.
   The numerical values of the matrix vary depending on where it is
   evaluated within the input space of the PolyMap. For this reason, the
   "jacob" variable holds a vector of "ncoord" PolyMaps. The outputs of
   each of these PolyMaps corresponds to a single column in the Jacobian
   matrix. */
   jacob = GetJacobian( this, status );

/* Get the number of points to be transformed. */
   npoint = astGetNpoint( out );

/* Get another PointSet to hold intermediate results. */
   work = astPointSet( npoint, ncoord, " ", status );

/* See if the PolyMap has been inverted.*/
   fwd = !astGetInvert( this );

/* Get pointers to the data arrays for all PointSets. Note, here "in" and
   "out" refer to inputs and outputs of the PolyMap (i.e. the forward
   transformation). These are respectively *outputs* and *inputs* of the
   inverse transformation. */
   ptr_in = astGetPoints( result );  /* Returned input positions */
   ptr_out = astGetPoints( out );    /* Supplied output positions */
   ptr_work = astGetPoints( work );  /* Work space */

/* Allocate an array of PointSets to hold the elements of the Jacobian
   matrix. */
   ptr_jac = astMalloc( sizeof( double ** )*ncoord );
   ps_jac = astCalloc( ncoord, sizeof( AstPointSet * ) );
   if( astOK ) {
      for( icoord = 0; icoord < ncoord; icoord++ ) {
         ps_jac[ icoord ] = astPointSet( npoint, ncoord, " ", status );
         ptr_jac[ icoord ] = astGetPoints( ps_jac[ icoord ] );
      }
   }

/* Allocate an array to hold flags indicating if each position has
   converged. Initialise it to hold zero at every element. */
   flags = astCalloc( npoint, sizeof( int ) );

/* Allocate memory to hold the Jacobian matrix at a single point. */
   mat = astMalloc( sizeof( double )*ncoord*ncoord );

/* Allocate memory to hold the offset vector. */
   vec = astMalloc( sizeof( double )*ncoord );

/* Allocate memory to hold work space for palDmat. */
   iw = astMalloc( sizeof( int )*ncoord );

/* Check pointers can be used safely. */
   if( astOK ) {

/* Store the initial guess at the required input positions. These are
   determined by transforming the supplied output positions using the
   inverse of a linear truncation of the PolyMap's forward
   transformation. */
      lintrunc = LinearGuess( this, status );
      (void) astTransform( lintrunc, out, 0, result );
      lintrunc = astAnnul( lintrunc );

/* Get the maximum number of iterations to perform. */
      maxiter = astGetNiterInverse( this );

/* Get the target relative error for the returned input axis values, and
   square it. */
      maxerr = astGetTolInverse( this );
      maxerr *= maxerr;

/* Initialise the number of positions which have reached the required
   accuracy. */
      nconv = 0;

/* Loop round doing iterations of a Newton-Raphson algorithm, until
   all points have achieved the required relative error, or the
   maximum number of iterations have been performed. */
      for( iter = 0; iter < maxiter && nconv < npoint && astOK; iter++ ) {

/* Use the original forward transformation of the supplied PolyMap to
   transform the current guesses at the required input positions into
   the corresponding output positions. Store the results in the "work"
   PointSet. */
         (void) astTransform( this, result, fwd, work );

/* Modify the work PointSet so that it holds the offsets from the output
   positions produced by the current input position guesses, and the
   required output positions. */
         for( icoord = 0; icoord < ncoord; icoord++ ) {
            pa = ptr_out[ icoord ];
            pb = ptr_work[ icoord ];
            for( ipoint = 0; ipoint< npoint; ipoint++,pa++,pb++ ) {
               if( *pa != AST__BAD && *pb != AST__BAD ){
                  *pb = *pa - *pb;
               } else {
                  *pb = AST__BAD;
               }
            }
         }

/* Evaluate the elements of the Jacobian matrix at the current input
   position guesses. */
         for( icoord = 0; icoord < ncoord; icoord++ ) {
            (void) astTransform( jacob[ icoord ], result, 1, ps_jac[ icoord ] );
         }

/* For each position, we now invert the matrix equation

    Dy = Jacobian.Dx

   to find a guess at the vector (dx) holding the offsets from the
   current input positions guesses to their required values. Loop over all
   points. */
         for( ipoint = 0; ipoint < npoint; ipoint++ ) {

/* Do not change positions that have already converged. */
            if( !flags[ ipoint ] ) {

/* Get the numerical values for the elements of the Jacobian matrix at
   the current point. */
               pa = mat;
               for( irow = 0; irow < ncoord; irow++ ) {
                  for( icol = 0; icol < ncoord; icol++ ) {
                     *(pa++) = ptr_jac[ icol ][ irow ][ ipoint ];
                  }

/* Store the offset from the current output position to the required
   output position. */
                  vec[ irow ] = ptr_work[ irow ][ ipoint ];
               }

/* Find the corresponding offset from the current input position to the required
   input position. */
               palDmat( ncoord, mat, vec, &det, &sing, iw );

/* If the matrix was singular, the input position cannot be evaluated so
   store a bad value for it and indicate it has converged. */
               if( sing ) {
                  for( icoord = 0; icoord < ncoord; icoord++ ) {
                     ptr_in[ icoord ][ ipoint ] = AST__BAD;
                  }
                  flags[ ipoint ] = 1;
                  nconv++;

/* Otherwise, update the input position guess. */
               } else {
                  vlensq = 0.0;
                  xlensq = 0.0;
                  pa = vec;
                  for( icoord = 0; icoord < ncoord; icoord++,pa++ ) {
                     xx = ptr_in[ icoord ][ ipoint ] + (*pa);
                     ptr_in[ icoord ][ ipoint ] = xx;
                     xlensq += xx*xx;
                     vlensq += (*pa)*(*pa);
                  }

/* Check for convergence. */
                  if( vlensq <= maxerr*xlensq ) {
                     flags[ ipoint ] = 1;
                     nconv++;
                  }
               }
            }
         }
      }
   }

/* Free resources. */
   vec = astFree( vec );
   iw = astFree( iw );
   mat = astFree( mat );
   flags = astFree( flags );
   work = astAnnul( work );

   if( ps_jac ) {
      for( icoord = 0; icoord < ncoord; icoord++ ) {
         ps_jac[ icoord ] = astAnnul( ps_jac[ icoord ] );
      }
      ps_jac = astFree( ps_jac );
   }

   ptr_jac = astFree( ptr_jac );

}

static AstMapping *LinearGuess( AstPolyMap *this, int *status ){
/*
*  Name:
*     LinearTruncation

*  Purpose:
*     Get a Mapping representing a linear approximation of a PolyMap

*  Type:
*     Private function.

*  Synopsis:
*     AstMapping *LinearGuess( AstPolyMap *this, int *status )

*  Description:
*     This function returns a linear Mapping approximating the original
*     forward transformation of the supplied PolyMap (i.e. the Invert
*     flag is assume to be zero). The linear and constant terms in the
*     supplied PolyMap are used for all outputs that have such terms. Any
*     other outputs are assumed to be equal to the corresponding inputs.
*     If the forward transformation of the PolyMap is defined, then it is
*     used to determine the returned Mapping. Otherwise, the inverse
*     transformation of the PolyMap is used. The forward transformation of
*     the returned Mapping always represents the forward transformation of
*     the original (i.e. uninverted) PolyMap.

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     The returned linear Mapping, or NULL if an error occurs.

*/

/* Local Variables: */
   AstMapping *result;           /* Pointer to returned Mapping */
   AstMapping *mm;               /* MatrixMap representing linear terms */
   AstMapping *sm;               /* ShiftMap representing offset terms */
   double **coeff;               /* Pointer to coefficient value arrays */
   double *matrix;               /* Linear scaling matrix */
   double *vector;               /* Offset vector */
   int ***power;                 /* Pointer to coefficient power arrays */
   int *ncoeff;                  /* Pointer to no. of coefficients */
   int flag;                     /* Indicates nature of current coeff */
   int ico;                      /* Coefficient index */
   int iin;                      /* Index of transformation input */
   int inverse;                  /* PolyMap inverse transformation selected? */
   int iout;                     /* Index of transformation output */
   int jin;                      /* Transformation input for current coeff */
   int nin;                      /* No. of i/ps for selected transformation */
   int nout;                     /* No. of o/ps for selected transformation */
   int ok;                       /* Current output is defined? */
   int pow;                      /* Power for current coeff and input */

/* Initialise returned value */
   result = NULL;

/* Check inherited status */
   if( !astOK ) return result;

/* If the required Mapping has already been determined, return a pointer
   to it. */
   if( this->lintrunc ) return astClone( this->lintrunc );

/* Get pointers to the arrays holding the required coefficient
   values and powers. Use the forward transformation of the original
   (uninverted) PolyMap if it is defined and the inverse transformation
   otherwise. */
   if ( this->ncoeff_f ){
      ncoeff = this->ncoeff_f;
      coeff = this->coeff_f;
      power = this->power_f;
      nin = astGetNin( this );
      nout = astGetNout( this );
      inverse = 0;

   } else {
      ncoeff = this->ncoeff_i;
      coeff = this->coeff_i;
      power = this->power_i;
      nout = astGetNin( this );
      nin = astGetNout( this );
      inverse = 1;
   }

/* Allocate memory to hold a matrix holding the linear terms of the
   PolyMap transformation selected above. Initialise it to hold zeros. */
   matrix = astCalloc( nin*nout, sizeof( double * ) );

/* Allocate memory to hold a vector holding the offset terms of the
   PolyMap transformation selected above. Initialise it to hold zeros. */
   vector = astCalloc( nout, sizeof( double * ) );
   if( astOK ){

/* Loop round all the outputs of the transformation selected above. */

      for( iout = 0; iout < nout; iout++ ){

/* Loop round all the coefficients for the current transformation output. */
         for( ico = 0; ico < ncoeff[iout]; ico++ ){

/* Initialise a tristate flag that indicates if the current coefficient is
   constant (0), linear (1) or something else (2). */
            flag = 0;

/* Look for the first input axis for which the current coeff uses a
   non-zero power. If the power is not one, the current coeff is neither
   linear nor constant, and so we can pass on to the next coeff. If any
   subsequent non-zero power is found, the coeff is not linear. */
            jin = 0;
            for( iin = 0; iin < nin; iin++ ){
               pow = power[iout][ico][iin];
               if( pow ) {
                  if( pow != 1 ) {
                     flag = 2;
                     break;
                  } else if( flag == 1 ) {
                     flag = 2;
                     break;
                  } else {
                     flag = 1;
                     jin = iin;
                  }
               }
            }

/* If the coeff is a constant term, insert it into the offset vector. */
            if( flag == 0 ) {
               vector[ iout ] = coeff[ iout ][ ico ];

/* If the coeff is a linear term, insert it into the matrix. */
            } else if( flag == 1 ) {
               matrix[ jin + iout*nin ] = coeff[ iout ][ ico ];
            }
         }
      }

/* If any PolyMap output has no linear terms, the matrix will not be
   invertable. So insert a 1.0 value on the diagonal of such rows. */
      for( iout = 0; iout < nout; iout++ ){
         ok = 0;
         for( iin = 0; iin < nin; iin++ ) {
            if( matrix[ iin + iout*nin ] != 0.0 ) {
               ok = 1;
               break;
            }
         }
         if( ! ok ) matrix[ iout + iout*nin ] = 1.0;
      }

/* Create a MatrixMap to represent the matrix. */
      mm = (AstMapping *) astMatrixMap( nin, nout, 0, matrix, " ", status );

/* If the MatrixMap cannot be inverted, use a unit map instead. */
      if( !astGetTranInverse( mm ) ) {
         mm = astAnnul( mm );
         mm = (AstMapping *) astUnitMap( nin, " ", status );
      }

/* Create a ShiftMap to represent the offset. */
      sm = (AstMapping *) astShiftMap( nout, vector, " ", status );

/* Combine them in series into a single compound Mapping. */
      result = (AstMapping *) astCmpMap( mm, sm, 1, " ", status );

/* If this CmpMap represents the inverse transformation of the original
   (i.e. uninverted) PolyMap, invert it so that it represents the forward
   transformation. */
      if( inverse ) astInvert( result );

/* Free resources */
      mm = astAnnul( mm );
      sm = astAnnul( sm );

/* Store the Mapping in the PolyMap structure so we do not need to
   recalculate it every time it is needed. */
      if( astOK ) this->lintrunc = astClone( result );
   }

/* Free resources */
   matrix = astFree( matrix );
   vector = astFree( vector );

/* Return the required Mapping. */
   return result;

}

static void LMFunc1D( const double *p, double *hx, int m, int n, void *adata ){
/*
*  Name:
*     LMFunc1D

*  Purpose:
*     Evaluate a test 1D polynomal.

*  Type:
*     Private function.

*  Synopsis:
*     void LMFunc1D( const double *p, double *hx, int m, int n, void *adata )

*  Description:
*     This function finds the residuals implied by a supplied set of
*     candidate polynomial coefficients. Each residual is a candidate
*     polynomial evaluated at one of the sample points, minus the
*     supplied target value for the polynomial at that test point.
*
*     The minimisation process minimises the sum of the squared residuals.

*  Parameters:
*     p
*        An array of "m" coefficients for the candidate polynomial. The
*        coefficients are ordered C0, C1, C2, etc.
*     hx
*        An array in which to return the "n" residuals. The residual at
*        sample "k" is returned in element (k).
*     m
*        The length of the "p" array. This should be equal to order.
*     n
*        The length of the "hx" array. This should be equal to nsamp.
*     adata
*        Pointer to a structure holding the sample positions and values,
*        and other information.

*/

/* Local Variables: */
   AstMinPackData *data;
   double *px1;
   double *py;
   const double *vp;
   double *vr;
   double res;
   int k;
   int w1;

/* Get a pointer to the data structure. */
   data = (AstMinPackData *) adata;

/* Initialise a pointer to the current returned residual value. */
   vr = hx;

/* Initialise a pointer to the sampled Y values for the polynomial output. */
   py = data->y[ 0 ];

/* Initialise a pointer to the powers of the input X values at the curent
   (i.e. first) sample. */
   px1 = data->xp1;

/* Loop over the index of the sample to which this residual refers. */
   for( k = 0; k < data->nsamp; k++ ) {

/* Initialise a pointer to the first coefficient (the  constant term) for the
   polynomial output coordinate. */
      vp = p;

/* Initialise this residual to hold the sampled Y value. Increment the
   pointer to the next sampled value for the current polynomial output. */
      res = -( *(py++) );

/* Loop over the coefficients. */
      for( w1 = 0; w1 < data->order; w1++ ) {

/* Increment the residual by the value of the current term Cw1*(x1^w1).
   Increment the pointer to the next coefficient (C). Also increment the
   pointer to the next higher power of X1. */
         res += ( *(vp++) )*( *(px1++) );
      }

/* Store the complete residual in the returned array, and increment the
   pointer to the next residual. */
      *(vr++) = res;
   }
}

static void LMFunc2D( const double *p, double *hx, int m, int n, void *adata ){
/*
*  Name:
*     LMFunc2D

*  Purpose:
*     Evaluate a test 2D polynomal.

*  Type:
*     Private function.

*  Synopsis:
*     void LMFunc2D( const double *p, double *hx, int m, int n, void *adata )

*  Description:
*     This function finds the residuals implied by a supplied set of
*     candidate polynomial coefficients. Each residual is a candidate
*     polynomial (either P1 or P2) evaluated at one of the sample points
*     (x1,x2), minus the supplied target value for the polynomial at that
*     test point.
*
*     The minimisation process minimises the sum of the squared residuals.

*  Parameters:
*     p
*        An array of "m" coefficients for the candidate polynomial. All the
*        coefficients for polynomial P1 come first, followed by those for P2.
*        Within each each polynomial, the coefficients are order C00, C10,
*        C01, C20, C11, C02, C30, C21, C12, C03, etc. So the coefficient
*        of (x1^j*x2^k) (=Cjk) for polynomial Pi is stored in element
*        [k + (j + k)*(j + k + 1)/2 + i*order*(order+1)/2] of the "p" array.
*     hx
*        An array in which to return the "n" residuals. The residual at
*        sample "k" for polynomial "i" is returned in element (k + nsamp*i).
*     m
*        The length of the "p" array. This should be equal to order*(order+1).
*     n
*        The length of the "hx" array. This should be equal to 2*nsamp.
*     adata
*        Pointer to a structure holding the sample positions and values,
*        and other information.

*/

/* Local Variables: */
   AstMinPackData *data;
   const double *vp0;
   const double *vp;
   double *py;
   double *px1;
   double *px10;
   double *px20;
   double *vr;
   double res;
   double *px2;
   int iout;
   int k;
   int w12;
   int w2;

/* Get a pointer to the data structure. */
   data = (AstMinPackData *) adata;

/* Initialise a pointer to the current returned residual value. */
   vr = hx;

/* Initialise a pointer to the first coefficient (the  constant term) for the
   current (i.e. first) polynomial output  coordinate. */
   vp0 = p;

/* Loop over each polynomial output coordinate. */
   for( iout = 0; iout < 2; iout++ ) {

/* Initialise a pointer to the sampled Y values for the first polynomial
   output. */
      py = data->y[ iout ];

/* Initialise pointers to the powers of the input X values at the curent
   (i.e. first) sample. */
      px10 = data->xp1;
      px20 = data->xp2;

/* Loop over the index of the sample to which this residual refers. */
      for( k = 0; k < data->nsamp; k++ ) {

/* Reset the pointer to the first coefficient (the  constant term)
   for the current polynomial output  coordinate. */
         vp = vp0;

/* Initialise this residual to hold the sampled Y value. Increment the
   pointer to the next sampled value for the current polynomial output. */
         res = -( *(py++) );

/* The w12 value is the sum of the powers of X1 and X2. So w12=0
   corresponds to the constant term in the polynomial, and (e.g.) w12=6
   corresponds to all the terms for which the sum of the powerss (w1+w2)
   is 6. Loop over all possible w12 values. */
         for( w12 = 0; w12 < data->order; w12++ ) {

/* The next coeff refers to (x1^w12)*(x2^0). Get pointers to the values
   holding x1^w12 and x2^0. */
            px1 = px10++;
            px2 = px20;

/* Loop over powers of x2. The corresponding power of x1 is "w12-w2", but
   is not explicitly needed here. So x1 moves down from w12 to zero, as
   w2 moves up from zero to w12. */
            for( w2 = 0; w2 <= w12; w2++ ) {

/* Increment the residual by the value of the current term Cw1w2*(x1^w1)*(x2^w2).
   Increment the pointer tio the next coefficient (C). Also decrement the
   pointer to the next lower power of X1, and increment the pointer to the next
   higher power of X2. */
               res += ( *(vp++) )*( *(px1--) )*( *(px2++) );
            }
         }

/* Move on to the x2 powers for the next sample. Don't need to do this
   for x1 since px10 is incremented within the above loop. */
         px20 += data->order;

/* Store the complete residual in the returned array, and increment the
   pointer to the next residual. */
         *(vr++) = res;
      }

/* Get a pointer to the first coefficient (the  constant term) for the
   next polynomial output  coordinate. */
      vp0 += data->order*( 1 + data->order )/2;
   }
}

static void LMJacob1D( const double *p, double *jac, int m, int n, void *adata ){
/*
*  Name:
*     LMJacob1D

*  Purpose:
*     Evaluate the Jacobian matrix of a test 1D polynomal.

*  Type:
*     Private function.

*  Synopsis:
*     void LMJacob1D( const double *p, double *jac, int m, int n, void *adata )

*  Description:
*     This function finds the Jacobian matrix that describes the rate of
*     change of every residual with respect to every polynomial coefficient.
*     Each residual is a candidate polynomial evaluated at one of the sample
*     points minus the supplied target value for the polynomial at that test
*     point.
*
*     For a polynomial the Jacobian matrix is constant (i.e. does not
*     depend on the values of the polynomial coefficients). So we only
*     evaluate it on the first call.

*  Parameters:
*     p
*        An array of "m" coefficients for the candidate polynomial.
*     jac
*        An array in which to return the "m*n" elements of the Jacobian
*        matrix. The rate of change of residual "r" with respect to
*        coefficient "c" is returned in element "r + c*n". The residual
*        at sample "k" of polynomial Pi has an "r" index of (k + nsamp*i).
*        The coefficient of (x1^j) for polynomial Pi has a "c" index
*        of j.
*     m
*        The length of the "p" array. This should be equal to order.
*     n
*        The number of residuals. This should be equal to nsamp.
*     adata
*        Pointer to a structure holding the sample positions and values,
*        and other information.

*/

/* Local Variables: */
   AstMinPackData *data;
   int k;
   int w1;

/* Get a pointer to the data structure. */
   data = (AstMinPackData *) adata;

/* The Jacobian of the residuals with respect to the polynomial
   coefficients is constant (i.e. does not depend on the values of the
   polynomial coefficients). So we only need to calculate it once. If
   this is the first call, calculate the Jacobian and return it in "jac".
   otherwise, just return immediately retaining the supplied "jac" values
   (which will be the values returned by the previous call to this
   function). */
   if( data->init_jac ) {
      data->init_jac = 0;

/* Loop over all residuals. */
      for( k = 0; k < n; k++ ) {

/* Loop over all parameters (i.e. polynomial coefficients). */
         for( w1 = 0; w1 < m; w1++ ) {

/* Store the Jacobian. */
            jac[ k + w1*n ] = (data->xp1)[ w1 + k*data->order ];
         }
      }
   }
}

static void LMJacob2D( const double *p, double *jac, int m, int n, void *adata ){
/*
*  Name:
*     LMJacob2D

*  Purpose:
*     Evaluate the Jacobian matrix of a test 2D polynomal.

*  Type:
*     Private function.

*  Synopsis:
*     void LMJacob2D( const double *p, double *jac, int m, int n, void *adata )

*  Description:
*     This function finds the Jacobian matrix that describes the rate of
*     change of every residual with respect to every polynomial coefficient.
*     Each residual is a candidate polynomial (either P1 or P2) evaluated
*     at one of the sample points (x1,x2), minus the supplied target value
*     for the polynomial at that test point.
*
*     For a polynomial the Jacobian matrix is constant (i.e. does not
*     depend on the values of the polynomial coefficients). So we only
*     evaluate it on the first call.

*  Parameters:
*     p
*        An array of "m" coefficients for the candidate polynomial. All the
*        coefficients for polynomial P1 come first, followed by those for P2.
*        Within each each polynomial, the coefficients are order C00, C10,
*        C01, C20, C11, C02, C30, C21, C12, C03, etc. So the coefficient
*        of (x1^j*x2^k) (=Cjk) for polynomial Pi is stored in element
*        [k + (j + k)*(j + k + 1)/2 + i*order*(order+1)/2] of the "p" array.
*     jac
*        An array in which to return the "m*n" elements of the Jacobian
*        matrix. The rate of change of residual "r" with respect to
*        coefficient "c" is returned in element "r + c*n". The residual
*        at sample "k" of polynomial Pi has an "r" index of (k + nsamp*i).
*        The coefficient of (x1^j*x2^k) for polynomial Pi has a "c" index
*        of [k + (j + k)*(j + k + 1)/2 + i*order*(order+1)/2].
*     m
*        The length of the "p" array. This should be equal to order*(order+1).
*     n
*        The number of residuals. This should be equal to 2*nsamp.
*     adata
*        Pointer to a structure holdin gthe sample positions and values,
*        and other information.

*/

/* Local Variables: */
   AstMinPackData *data;
   int iout;
   int k;
   int ncof;
   int vp;
   int vr;
   int w1;
   int w12;
   int w2;

/* Get a pointer to the data structure. */
   data = (AstMinPackData *) adata;

/* The Jacobian of the residuals with respect to the polynomial
   coefficients is constant (i.e. does not depend on the values of the
   polynomial coefficients). So we only need to calculate it once. If
   this is the first call, calculate the Jacobian and return it in "jac".
   otherwise, just return immediately retaining the supplied "jac" values
   (which will be the values returned by the previous call to this
   function). */
   if( data->init_jac ) {
      data->init_jac = 0;

/* Store the number of coefficients in one polynomial. */
      ncof = data->order*( 1 + data->order )/2;

/* Loop over all residuals. */
      for( vr = 0; vr < n; vr++ ) {

/* Calculate the polynomial output index, and sample index, that creates
   the current residual. */
         iout = vr/data->nsamp;
         k = vr - iout*data->nsamp;

/* Loop over all parameters (i.e. polynomial coefficients). */
         for( vp = 0; vp < m; vp++ ) {

/* If this coefficient is not used in the creation of the current
   polynomial output value, then the Jacobian value is zero. */
            if( vp/ncof != iout ) {
               jac[ vr + vp*n ] = 0.0;

/* Otherwise, get the powers of the two polynomial inputs, to which
   the current coefficient relates. */
            } else {
               w12 = ( -1.0 + sqrt( 1.0 + 8.0*(vp - iout*ncof) ) )/2.0;
               w2 = vp - iout*ncof - w12*( w12 + 1 )/2;
               w1 = w12 - w2;

/* Store the Jacobian. */
               jac[ vr + vp*n ] = (data->xp1)[ w1 + k*data->order ]*
                                  (data->xp2)[ w2 + k*data->order ];
            }
         }
      }
   }
}

#if defined(THREAD_SAFE)
static int ManageLock( AstObject *this_object, int mode, int extra,
                       AstObject **fail, int *status ) {
/*
*  Name:
*     ManageLock

*  Purpose:
*     Manage the thread lock on an Object.

*  Type:
*     Private function.

*  Synopsis:
*     #include "object.h"
*     AstObject *ManageLock( AstObject *this, int mode, int extra,
*                            AstObject **fail, int *status )

*  Class Membership:
*     PolyMap member function (over-rides the astManageLock protected
*     method inherited from the parent class).

*  Description:
*     This function manages the thread lock on the supplied Object. The
*     lock can be locked, unlocked or checked by this function as
*     deteremined by parameter "mode". See astLock for details of the way
*     these locks are used.

*  Parameters:
*     this
*        Pointer to the Object.
*     mode
*        An integer flag indicating what the function should do:
*
*        AST__LOCK: Lock the Object for exclusive use by the calling
*        thread. The "extra" value indicates what should be done if the
*        Object is already locked (wait or report an error - see astLock).
*
*        AST__UNLOCK: Unlock the Object for use by other threads.
*
*        AST__CHECKLOCK: Check that the object is locked for use by the
*        calling thread (report an error if not).
*     extra
*        Extra mode-specific information.
*     fail
*        If a non-zero function value is returned, a pointer to the
*        Object that caused the failure is returned at "*fail". This may
*        be "this" or it may be an Object contained within "this". Note,
*        the Object's reference count is not incremented, and so the
*        returned pointer should not be annulled. A NULL pointer is
*        returned if this function returns a value of zero.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*    A local status value:
*        0 - Success
*        1 - Could not lock or unlock the object because it was already
*            locked by another thread.
*        2 - Failed to lock a POSIX mutex
*        3 - Failed to unlock a POSIX mutex
*        4 - Bad "mode" value supplied.

*  Notes:
*     - This function attempts to execute even if an error has already
*     occurred.
*/

/* Local Variables: */
   AstPolyMap *this;
   int ic;
   int result;
   int nc;

/* Initialise */
   result = 0;

/* Check the supplied pointer is not NULL. */
   if( !this_object ) return result;

/* Obtain a pointer to the PolyMap structure. */
   this = (AstPolyMap *) this_object;

/* Invoke the ManageLock method inherited from the parent class. */
   if( !result ) result = (*parent_managelock)( this_object, mode, extra,
                                                fail, status );

/* Invoke the astManageLock method on any Objects contained within
   the supplied Object. */
   if( this->jacobian ) {
      nc = astGetNin( this );
      for( ic = 0; ic < nc && result; ic++ ) {
         result = astManageLock( (this->jacobian)[ ic ], mode,
                                 extra, fail );
      }
   }

   return result;

}
#endif

static int MapMerge( AstMapping *this, int where, int series, int *nmap,
                     AstMapping ***map_list, int **invert_list, int *status ) {
/*
*  Name:
*     MapMerge

*  Purpose:
*     Simplify a sequence of Mappings containing a PolyMap.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mapping.h"
*     int MapMerge( AstMapping *this, int where, int series, int *nmap,
*                   AstMapping ***map_list, int **invert_list, int *status )

*  Class Membership:
*     PolyMap method (over-rides the protected astMapMerge method
*     inherited from the Mapping class).

*  Description:
*     This function attempts to simplify a sequence of Mappings by
*     merging a nominated PolyMap in the sequence with its neighbours,
*     so as to shorten the sequence if possible.
*
*     In many cases, simplification will not be possible and the
*     function will return -1 to indicate this, without further
*     action.
*
*     In most cases of interest, however, this function will either
*     attempt to replace the nominated PolyMap with a Mapping which it
*     considers simpler, or to merge it with the Mappings which
*     immediately precede it or follow it in the sequence (both will
*     normally be considered). This is sufficient to ensure the
*     eventual simplification of most Mapping sequences by repeated
*     application of this function.
*
*     In some cases, the function may attempt more elaborate
*     simplification, involving any number of other Mappings in the
*     sequence. It is not restricted in the type or scope of
*     simplification it may perform, but will normally only attempt
*     elaborate simplification in cases where a more straightforward
*     approach is not adequate.

*  Parameters:
*     this
*        Pointer to the nominated PolyMap which is to be merged with
*        its neighbours. This should be a cloned copy of the PolyMap
*        pointer contained in the array element "(*map_list)[where]"
*        (see below). This pointer will not be annulled, and the
*        PolyMap it identifies will not be modified by this function.
*     where
*        Index in the "*map_list" array (below) at which the pointer
*        to the nominated PolyMap resides.
*     series
*        A non-zero value indicates that the sequence of Mappings to
*        be simplified will be applied in series (i.e. one after the
*        other), whereas a zero value indicates that they will be
*        applied in parallel (i.e. on successive sub-sets of the
*        input/output coordinates).
*     nmap
*        Address of an int which counts the number of Mappings in the
*        sequence. On entry this should be set to the initial number
*        of Mappings. On exit it will be updated to record the number
*        of Mappings remaining after simplification.
*     map_list
*        Address of a pointer to a dynamically allocated array of
*        Mapping pointers (produced, for example, by the astMapList
*        method) which identifies the sequence of Mappings. On entry,
*        the initial sequence of Mappings to be simplified should be
*        supplied.
*
*        On exit, the contents of this array will be modified to
*        reflect any simplification carried out. Any form of
*        simplification may be performed. This may involve any of: (a)
*        removing Mappings by annulling any of the pointers supplied,
*        (b) replacing them with pointers to new Mappings, (c)
*        inserting additional Mappings and (d) changing their order.
*
*        The intention is to reduce the number of Mappings in the
*        sequence, if possible, and any reduction will be reflected in
*        the value of "*nmap" returned. However, simplifications which
*        do not reduce the length of the sequence (but improve its
*        execution time, for example) may also be performed, and the
*        sequence might conceivably increase in length (but normally
*        only in order to split up a Mapping into pieces that can be
*        more easily merged with their neighbours on subsequent
*        invocations of this function).
*
*        If Mappings are removed from the sequence, any gaps that
*        remain will be closed up, by moving subsequent Mapping
*        pointers along in the array, so that vacated elements occur
*        at the end. If the sequence increases in length, the array
*        will be extended (and its pointer updated) if necessary to
*        accommodate any new elements.
*
*        Note that any (or all) of the Mapping pointers supplied in
*        this array may be annulled by this function, but the Mappings
*        to which they refer are not modified in any way (although
*        they may, of course, be deleted if the annulled pointer is
*        the final one).
*     invert_list
*        Address of a pointer to a dynamically allocated array which,
*        on entry, should contain values to be assigned to the Invert
*        attributes of the Mappings identified in the "*map_list"
*        array before they are applied (this array might have been
*        produced, for example, by the astMapList method). These
*        values will be used by this function instead of the actual
*        Invert attributes of the Mappings supplied, which are
*        ignored.
*
*        On exit, the contents of this array will be updated to
*        correspond with the possibly modified contents of the
*        "*map_list" array.  If the Mapping sequence increases in
*        length, the "*invert_list" array will be extended (and its
*        pointer updated) if necessary to accommodate any new
*        elements.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     If simplification was possible, the function returns the index
*     in the "map_list" array of the first element which was
*     modified. Otherwise, it returns -1 (and makes no changes to the
*     arrays supplied).

*  Notes:
*     - A value of -1 will be returned if this function is invoked
*     with the global error status set, or if it should fail for any
*     reason.
*/

/* Local Variables: */
   AstPolyMap *pmap0;    /* Nominated PolyMap */
   AstPolyMap *pmap1;    /* Neighbouring PolyMap */
   int i;                /* Index of neighbour */
   int nin;              /* Number of input coordinates for nominated PolyMap */
   int nout;             /* Number of output coordinates for nominated PolyMap */
   int ok;               /* Are PolyMaps equivalent? */
   int oldinv0;          /* Original Invert value in pmap0 */
   int oldinv1;          /* Original Invert value in pmap1 */
   int result;           /* Result value to return */

/* Initialise. */
   result = -1;

/* Check the global error status. */
   if ( !astOK ) return result;

/* Save a pointer to the nominated PolyMap. */
   pmap0 = (AstPolyMap *) ( *map_list )[ where ];

/* The only simplification which can currently be performed is to merge a PolyMap
   with its own inverse. This can only be done in series. Obviously,
   there are potentially other simplications which could be performed, but
   time does not currently allow these to be coded. */
   if( series ) {

/* Temporarily set the Invert flag of the nominated PolyMap to the
   required value, first saving the original value so that it can be
   re-instated later. */
      oldinv0 = astGetInvert( pmap0 );
      astSetInvert( pmap0, ( *invert_list )[ where ] );

/* Get the number of inputs and outputs to the nominated PolyMap. */
      nin = astGetNin( pmap0 );
      nout = astGetNout( pmap0 );

/* Check each neighbour. */
      for( i = where - 1; i <= where + 1; i += 2 ) {

/* Continue with the next pass if the neighbour does not exist. */
         if( i < 0 || i >= *nmap ) continue;

/* Continue with the next pass if this neighbour is not a PermMap. */
         if( ! astIsAPolyMap( ( *map_list )[ i ] ) ) continue;

/* Get a pointer to it. */
         pmap1 = (AstPolyMap *) ( *map_list )[ i ];

/* Check it is used in the opposite direction to the nominated PolyMap. */
         if( ( *invert_list )[ i ] == ( *invert_list )[ where ] ) continue;

/* We use the astEqual method to check that the two PolyMaps are equal.
   But at the moment they may not be equal because they may have
   different Invert flags. Therefore, temporarily set the invert flag
   of the neighbour so that it is the same as the nominated PolyMap,
   first saving the original value so that it can be re-instated later.
   Note, we have already checked that the two PolyMaps are used in opposite
   directions within the CmpMap. */
         oldinv1 = astGetInvert( pmap1 );
         astSetInvert( pmap1, ( *invert_list )[ where ] );

/* Use astEqual to check that the coefficients etc are equal in the two
   PolyMaps. */
         ok = astEqual( pmap0, pmap1 );

/* Re-instate the original value of the Invert flag in the neighbour. */
         astSetInvert( pmap1, oldinv1 );

/* Pass on to the next neighbour if the current neighbour differs from
   the nominated PolyMap. */
         if( !ok ) continue;

/* If we get this far, then the nominated PolyMap and the current
   neighbour cancel each other out, so replace each by a UnitMap. */
         pmap0 = astAnnul( pmap0 );
         pmap1 = astAnnul( pmap1 );
         if( i < where ) {
            ( *map_list )[ where ] = (AstMapping *) astUnitMap( nout, "", status );
            ( *map_list )[ i ] = (AstMapping *) astUnitMap( nout, "", status );
            ( *invert_list )[ where ] = 0;
            ( *invert_list )[ i ] = 0;
            result = i;
         } else {
            ( *map_list )[ where ] = (AstMapping *) astUnitMap( nin, "", status );
            ( *map_list )[ i ] = (AstMapping *) astUnitMap( nin, "", status );
            ( *invert_list )[ where ] = 0;
            ( *invert_list )[ i ] = 0;
            result = where;
         }

/* Leave the loop. */
         break;
      }

/* If the nominated PolyMap was not replaced by a UnitMap, then
   re-instate its original value for the Invert flag. */
      if( pmap0 ) astSetInvert( pmap0, oldinv0 );
   }

/* Return the result. */
   return result;
}

static int MPFunc1D( void *p, int m, int n, const double *x, double *fvec,
                     double *fjac, int ldfjac, int iflag ) {
/*
*  Name:
*     MPFunc1D

*  Purpose:
*     Evaluate a test 1D polynomal or its Jacobian.

*  Type:
*     Private function.

*  Synopsis:
*     int MPFunc1D( void *p, int m, int n, const double *x, double *fvec,
*                   double *fjac, int ldfjac, int iflag )

*  Description:
*     This function finds the residuals or Jacobian implied by a supplied
*     set of candidate polynomial coefficients. Each residual is a candidate
*     polynomial evaluated at one of the sample points, minus the
*     supplied target value for the polynomial at that test point.
*
*     The minimisation process minimises the sum of the squared residuals.

*  Parameters:
*     p
*        A pointer to data needed to evaluate the required residuals or
*        Jacobians.
*     m
*        The number of residuals.
*     n
*        The number of coefficients.
*     x
*        An array of "n" coefficients for the candidate polynomial.
*     fvec
*        An array in which to return the "m" residuals.
*     fjac
*        An array in which to return the "mxn" Jacobian values.
*     ldflac
*        Unused (should always be equal to "m").
*     iflag
*        If 1 calculate the residuals, otherwise calculate the Jacobians.

*/

/* If required, calculate the function values at X, and return this
   vector in fvec. Do not alter fjac. */
   if( iflag == 1 ) {
      LMFunc1D( x, fvec, n, m, p );

/* Otherwise, calculate the Jacobian values at X, and return this
   matrix in fjac. Do not alter fvec. */
   } else {
      LMJacob1D( x, fjac, n, m, p );
   }

   return 0;
}

static int MPFunc2D( void *p, int m, int n, const double *x, double *fvec,
                     double *fjac, int ldfjac, int iflag ) {
/*
*  Name:
*     MPFunc1D

*  Purpose:
*     Evaluate a test 2D polynomal or its Jacobian.

*  Type:
*     Private function.

*  Synopsis:
*     int MPFunc2D( void *p, int m, int n, const double *x, double *fvec,
*                   double *fjac, int ldfjac, int iflag )

*  Description:
*     This function finds the residuals or Jacobian implied by a supplied
*     set of candidate polynomial coefficients. Each residual is a candidate
*     polynomial evaluated at one of the sample points, minus the
*     supplied target value for the polynomial at that test point.
*
*     The minimisation process minimises the sum of the squared residuals.

*  Parameters:
*     p
*        A pointer to data needed to evaluate the required residuals or
*        Jacobians.
*     m
*        The number of residuals.
*     n
*        The number of coefficients.
*     x
*        An array of "n" coefficients for the candidate polynomial.
*     fvec
*        An array in which to return the "m" residuals.
*     fjac
*        An array in which to return the "mxn" Jacobian values.
*     ldflac
*        Unused (should always be equal to "m").
*     iflag
*        If 1 calculate the residuals, otherwise calculate the Jacobians.

*/

/* If required, calculate the function values at X, and return this
   vector in fvec. Do not alter fjac. */
   if( iflag == 1 ) {
      LMFunc2D( x, fvec, n, m, p );

/* Otherwise, calculate the Jacobian values at X, and return this
   matrix in fjac. Do not alter fvec. */
   } else {
      LMJacob2D( x, fjac, n, m, p );

   }

   return 0;
}

static void PolyCoeffs( AstPolyMap *this, int forward, int nel, double *coeffs,
                        int *ncoeff, int *status ){
/*
*++
*  Name:
c     astPolyCoeffs
f     AST_POLYCOEFFS

*  Purpose:
*     Retrieve the coefficient values used by a PolyMap.

*  Type:
*     Public function.

*  Synopsis:
c     #include "polymap.h"
c     void astPolyCoeffs( AstPolyMap *this, int forward, int nel, double *coeffs,
c                         int *ncoeff )
f     CALL AST_POLYCOEFFS( THIS, FORWARD, NEL, COEFFS, NCOEFF, STATUS )

*  Class Membership:
*     PolyMap method.

*  Description:
*     This function returns the coefficient values used by either the
*     forward or inverse transformation of a PolyMap, in the same form
*     that they are supplied to the PolyMap constructor.
*
*     Usually, you should call this method first with
c     "nel"
f     NEL
*     set to zero to determine the number of coefficients used by the
*     PolyMap. This allows you to allocate an array of the correct size to
*     hold all coefficient data. You should then call this method a
*     second time to get the coefficient data.

*  Parameters:
c     this
f     THIS = INTEGER (Given)
*        Pointer to the original Mapping.
c     forward
f     FORWARD = LOGICAL (Given)
c        If non-zero,
f        If .TRUE.,
*        the coefficients of the forward PolyMap transformation are
*        returned. Otherwise the inverse transformation coefficients are
*        returned.
c     nel
f     NEL = INTEGER (Given)
*        The length of the supplied
c        "coeffs"
f        COEFFS
*        array. It should be at least "ncoeff*( nin + 2 )" if
c        "foward" is non-zero,
f        FORWARD is .TRUE.,
*        and "ncoeff*( nout + 2 )" otherwise, where "ncoeff" is the
*        number of coefficients to be returned. If a value of zero
*        is supplied, no coefficient values are returned, but the
*        number of coefficients used by the transformation is still
*        returned in
c        "ncoeff".
f        NCOEFF.
c     coeffs
f     COEFFS( NEL ) = DOUBLE PRECISION (Returned)
*        An array in which to return the coefficients used by the
*        requested transformation of the PolyMap. Ignored if
c        "nel" is zero.
f        NEL is zero.
*        The coefficient data is returned in the form in which it is
*        supplied to the PolyMap constructor. That is, each group of
*        "2 + nin" or "2 + nout" adjacent elements describe a single
*        coefficient of the forward or inverse transformation. See the
*        PolyMap constructor documentation for further details.
*
*        If the supplied array is too short to hold all the coefficients,
*        trailing coefficients are excluded. If the supplied array is
*        longer than needed to hold all the coefficients, trailing
*        elements are filled with zeros.
c     ncoeff
f     NCOEFF = INTEGER (Returned)
*        The number of coefficients used by the requested transformation.
*        A value of zero is returned if the transformation does not
*        have any defining polynomials. A value is returned for this argument
*        even if
c        "nel" is zero.
f        NEL is zero.
f     STATUS = INTEGER (Given and Returned)
f        The global status.

*--
*/

/* Local Variables: */
   double **coeff;
   int ***power;
   int *nco;
   int *pp;
   int iax;
   int icoeff;
   int iel;
   int ipoly;
   int nax;
   int npoly;

/* Initialise */
   *ncoeff = 0;

/* Check the inherited status. */
   if ( !astOK ) return;

/* Fill any supplied array with zeros. */
   if( nel ) memset( coeffs, 0, nel*sizeof( *coeffs ) );

/* Get the values to use, taking account of whether the PolyMap has been
   inverted or not. */
   if( forward != astGetInvert( this ) ){
      nco = this->ncoeff_f;
      power = this->power_f;
      coeff = this->coeff_f;
      npoly = astGetNout( this );
      nax = astGetNin( this );
   } else {
      nco = this->ncoeff_i;
      power = this->power_i;
      coeff = this->coeff_i;
      npoly = astGetNin( this );
      nax = astGetNout( this );
   }

/* Notheg to do if there are no coeffs. */
   if( nco && power && coeff ) {

/* Initialise index of next value to store in returned array. */
      iel = 0;

/* Loop round each 1D polynomial. */
      for( ipoly = 0; ipoly < npoly; ipoly++ ){

/* Loop round coefficients. */
         for( icoeff = 0; icoeff < nco[ ipoly ]; icoeff++ ) {

/* Store the coefficient value in the next element of the returned array,
   if the array is not already full. Increment the pointer to the next
   element. */
            if( iel < nel ) coeffs[ iel++ ] = coeff[ipoly][icoeff];

/* Next store the integer index of the PolyMap input or output which uses
   the coefficient within its defining polynomial (the first axis has
   index 1). */
            if( iel < nel ) coeffs[ iel++ ] = ipoly + 1;

/* The remaining elements of the group give the integer powers to use
   with each input or output coordinate value. */
            pp = power[ipoly][icoeff];
            for( iax = 0; iax < nax; iax++,pp++ ) {
               if( iel < nel ) coeffs[ iel++ ] = *pp;
            }
         }

/* Increment the total number of coefficients used by the transformation. */
         *ncoeff += nco[ ipoly ];
      }
   }
}

static void PolyPowers( AstPolyMap *this, double **work, int ncoord,
                        const int *mxpow, double **ptr, int point,
                        int fwd, int *status ){
/*
*+
*  Name:
*     astPolyPowers

*  Purpose:
*     Find the required powers of the input axis values.

*  Type:
*     Protected function.

*  Synopsis:
*     #include "polymap.h"
*     void astPolyPowers( AstPolyMap *this, double **work, int ncoord,
*                         const int *mxpow, double **ptr, int point,
*                         int fwd )

*  Class Membership:
*     PolyMap virtual function.

*  Description:
*     This function is used by astTransform to calculate the powers of
*     the axis values for a single input position. In the case of
*     sub-classes, the powers may not be simply powers of the supplied
*     axis values but may be more complex quantities such as a Chebyshev
*     polynomial of the required degree evaluated at the input axis values.

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     work
*        An array of "ncoord" pointers, each pointing to an array of
*        length "max(2,mxpow)". The required values are placed in this
*        array on exit.
*     ncoord
*        The number of axes.
*     mxpow
*        Pointer to an array holding the maximum power required of each
*        axis value. Should have "ncoord" elements.
*     ptr
*        An array of "ncoord" pointers, each pointing to an array holding
*        the axis values. Each of these arrays of axis values must have
*        at least "point+1" elements.
*     point
*        The zero based index of the point within "ptr" that holds the
*        axis values to be exponentiated.
*     fwd
*        Do the supplied coefficients define the foward transformation of
*        the PolyMap?
*-
*/

/* Local Variables; */
   double *pwork;
   double x;
   int coord;
   int ip;

/* Check the local error status. */
   if ( !astOK ) return;

/* For the base PolyMap class, this method simply raises each input axis
   value to the required power. Loop over all input axes. */
   for( coord = 0; coord < ncoord; coord++ ) {

/* Get a pointer to the array in which the powers of the current axis
   value are to be returned. */
      pwork = work[ coord ];

/* Anything to the power zero is 1.0. */
      pwork[ 0 ] = 1.0;

/* Get the input axis value. If it is bad, store bad values for all
   remaining powers. */
      x = ptr[ coord ][ point ];
      if( x == AST__BAD ) {
         for( ip = 1; ip <= mxpow[ coord ]; ip++ ) pwork[ ip ] = AST__BAD;

/* Otherwise, form and store the required powers of the input axis value. */
      } else {
         for( ip = 1; ip <= mxpow[ coord ]; ip++ ) {
            pwork[ ip ] = pwork[ ip - 1 ]*x;
         }
      }
   }
}

static AstPolyMap *PolyTran( AstPolyMap *this, int forward, double acc,
                             double maxacc, int maxorder, const double *lbnd,
                             const double *ubnd, int *status ){
/*
*++
*  Name:
c     astPolyTran
f     AST_POLYTRAN

*  Purpose:
*     Fit a PolyMap inverse or forward transformation.

*  Type:
*     Public function.

*  Synopsis:
c     #include "polymap.h"
c     AstPolyMap *astPolyTran( AstPolyMap *this, int forward, double acc,
c                              double maxacc, int maxorder, const double *lbnd,
c                              const double *ubnd )
f     RESULT = AST_POLYTRAN( THIS, FORWARD, ACC, MAXACC, MAXORDER, LBND,
f                            UBND, STATUS )

*  Class Membership:
*     PolyMap method.

*  Description:
*     This function creates a new PolyMap which is a copy of the supplied
*     PolyMap, in which a specified transformation (forward or inverse)
*     has been replaced by a new polynomial transformation. The
*     coefficients of the new transformation are estimated by sampling
*     the other transformation and performing a least squares polynomial
*     fit in the opposite direction to the sampled positions and values.
*
*     This method can only be used on (1-input,1-output) or (2-input,2-output)
*     PolyMaps.
*
*     The transformation to create is specified by the
c     "forward" parameter.
f     FORWARD parameter.
*     In what follows "X" refers to the inputs of the PolyMap, and "Y" to
*     the outputs of the PolyMap. The forward transformation transforms
*     input values (X) into output values (Y), and the inverse transformation
*     transforms output values (Y) into input values (X). Within a PolyMap,
*     each transformation is represented by an independent set of
*     polynomials, P_f or P_i: Y=P_f(X) for the forward transformation and
*     X=P_i(Y) for the inverse transformation.
*
c     The "forward"
f     The FORWARD
*     parameter specifies the transformation to be replaced. If it is
c     non-zero,
f     is .TRUE.,
*     a new forward transformation is created
*     by first finding the input values (X) using the inverse transformation
*     (which must be available) at a regular grid of points (Y) covering a
*     rectangular region of the PolyMap's output space. The coefficients of
*     the required forward polynomial, Y=P_f(X), are chosen in order to
*     minimise the sum of the squared residuals between the sampled values
*     of Y and P_f(X).
*
c     If "forward" is zero (probably the most likely case),
f     If FORWARD is .FALSE. (probably the most likely case),
*     a new inverse transformation is created by
*     first finding the output values (Y) using the forward transformation
*     (which must be available) at a regular grid of points (X) covering a
*     rectangular region of the PolyMap's input space. The coefficients of
*     the required inverse polynomial, X=P_i(Y), are chosen in order to
*     minimise the sum of the squared residuals between the sampled values
*     of X and P_i(Y).
*
*     This fitting process is performed repeatedly with increasing
*     polynomial orders (starting with linear) until the target
*     accuracy is achieved, or a specified maximum order is reached. If
*     the target accuracy cannot be achieved even with this maximum-order
*     polynomial, the best fitting maximum-order polynomial is returned so
*     long as its accuracy is better than
c     "maxacc".
f     MAXACC.
*     If it is not, a NULL pointer is returned but no error is reported.

*  Parameters:
c     this
f     THIS = INTEGER (Given)
*        Pointer to the original Mapping.
c     forward
f     FORWARD = LOGICAL (Given)
c        If non-zero,
f        If .TRUE.,
*        the forward PolyMap transformation is replaced. Otherwise the
*        inverse transformation is replaced.
c     acc
f     ACC = DOUBLE (Given)
*        The target accuracy, expressed as a geodesic distance within
*        the PolyMap's input space (if
c        "forward" is zero)  or output space (if "forward" is non-zero).
f        FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.).
c     maxacc
f     MAXACC = DOUBLE (Given)
*        The maximum allowed accuracy for an acceptable polynomial,
*        expressed as a geodesic distance within the PolyMap's input
*        space (if
c        "forward" is zero)  or output space (if "forward" is non-zero).
f        FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.).
c     maxorder
f     MAXORDER = INTEGER (Given)
*        The maximum allowed polynomial order. This is one more than the
*        maximum power of either input axis. So for instance, a value of
*        3 refers to a quadratic polynomial. Note, cross terms with total
*        powers greater than or equal to
c        maxorder
f        MAXORDER
*        are not inlcuded in the fit. So the maximum number of terms in
*        each of the fitted polynomials is
c        maxorder*(maxorder+1)/2.
f        MAXORDER*(MAXORDER+1)/2.
c     lbnd
f     LBND( * ) = DOUBLE PRECISION (Given)
c        Pointer to an
f        An
*        array holding the lower bounds of a rectangular region within
*        the PolyMap's input space (if
c        "forward" is zero)  or output space (if "forward" is non-zero).
f        FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.).
*        The new polynomial will be evaluated over this rectangle. The
*        length of this array should equal the value of the PolyMap's Nin
*        or Nout attribute, depending on
c        "forward".
f        FORWARD.
c     ubnd
f     UBND( * ) = DOUBLE PRECISION (Given)
c        Pointer to an
f        An
*        array holding the upper bounds of a rectangular region within
*        the PolyMap's input space (if
c        "forward" is zero)  or output space (if "forward" is non-zero).
f        FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.).
*        The new polynomial will be evaluated over this rectangle.  The
*        length of this array should equal the value of the PolyMap's Nin
*        or Nout attribute, depending on
c        "forward".
f        FORWARD.
f     STATUS = INTEGER (Given and Returned)
f        The global status.

*  Returned Value:
c     astPolyTran()
f     AST_POLYTRAN = INTEGER
*        A pointer to the new PolyMap.
c        A NULL pointer
f        AST__NULL
*        will be returned if the fit fails to achieve the accuracy
*        specified by
c        "maxacc",
f        MAXACC,
*        but no error will be reported.

*  Applicability:
*     PolyMap
*        All PolyMaps have this method.
*     ChebyMap
c        The ChebyMap implementation of this method allows
c        NULL pointers to be supplied for "lbnd" and/or "ubnd",
c        in which case the corresponding bounds supplied when the ChebyMap
c        was created are used.
*        The returned PolyMap will be a ChebyMap, and the new transformation
*        will be defined as a weighted sum of Chebyshev functions of the
*        first kind.

*  Notes:
*     - The IterInverse attribute is always cleared in the returned PolyMap.
*     This means that the returned PolyMap will always use the new fit by
*     default, rather than the iterative inverse, regardless of the setting
*     of IterInverse in the supplied PolyMap.
*     - This function can only be used on 1D or 2D PolyMaps which have
*     the same number of inputs and outputs.
*     - A null Object pointer (AST__NULL) will be returned if this
c     function is invoked with the AST error status set, or if it
f     function is invoked with STATUS set to an error value, or if it
*     should fail for any reason.
*--
*/

/* Local Variables: */
   AstPolyMap *result;
   int ok;

/* Initialise. */
   result = NULL;

/* Check the inherited status. */
   if ( !astOK ) return result;

/* Take a copy of the supplied PolyMap. */
   result = astCopy( this );

/* Replace the required transformation. */
   ok = ReplaceTransformation( result, forward, acc, maxacc, maxorder, lbnd,
                               ubnd, status );

/* If an error occurred, or the fit was not good enough, annul the returned
   PolyMap. */
   if ( !ok || !astOK ) {
      result = astAnnul( result );

/* Otherwise, ensure the new fit is used in preference to the iterative
   inverse by clearing the IterInverse attribute. */
   } else {
      astClearIterInverse( result );
   }

/* Return the result. */
   return result;
}

static int ReplaceTransformation( AstPolyMap *this, int forward, double acc,
                                  double maxacc, int maxorder, const double *lbnd,
                                  const double *ubnd, int *status ){
/*
*  Name:
*     ReplaceTransformation

*  Purpose:
*     Create a new inverse or forward transformation for a PolyMap.

*  Type:
*     Private function.

*  Synopsis:
*     int ReplaceTransformation( AstPolyMap *this, int forward, double acc,
*                                double maxacc, int maxorder, const double *lbnd,
*                                const double *ubnd, int *status )

*  Description:
*     This function creates a new forward or inverse transformation for
*     the supplied PolyMap (replacing any existing transformation), by
*     sampling the other transformation and performing a least squares
*     polynomial fit to the sample positions and values.
*
*     The transformation to create is specified by the "forward" parameter.
*     In what follows "X" refers to the inputs of the PolyMap, and "Y" to
*     the outputs of the PolyMap. The forward transformation transforms
*     input values (X) into output values (Y), and the inverse transformation
*     transforms output values (Y) into input values (X). Within a PolyMap,
*     each transformation is represented by an independent set of
*     polynomials: Y=P_f(X) for the forward transformation and X=P_i(Y)
*     for the inverse transformation.
*
*     If "forward" is zero then a new inverse transformation is created by
*     first finding the output values (Y) using the forward transformation
*     (which must be available) at a regular grid of points (X) covering a
*     rectangular region of the PolyMap's input space. The coefficients of
*     the required inverse polynomial, X=P_i(Y), are chosen in order to
*     minimise the sum of the squared residuals between the sampled values
*     of X and P_i(Y).
*
*     If "forward" is non-zero then a new forward transformation is created
*     by first finding the input values (X) using the inverse transformation
*     (which must be available) at a regular grid of points (Y) covering a
*     rectangular region of the PolyMap's output space. The coefficients of
*     the required forward polynomial, Y=P_f(X), are chosen in order to
*     minimise the sum of the squared residuals between the sampled values
*     of Y and P_f(X).
*
*     This fitting process is performed repeatedly with increasing
*     polynomial orders (starting with linear) until the target
*     accuracy is achieved, or a specified maximum order is reached. If
*     the target accuracy cannot be achieved even with this maximum-order
*     polynomial, the best fitting maximum-order polynomial is returned so
*     long as its accuracy is better than "maxacc".

*  Parameters:
*     this
*        The PolyMap.
*     forward
*        If non-zero, then the forward PolyMap transformation is
*        replaced. Otherwise the inverse transformation is replaced.
*     acc
*        The target accuracy, expressed as a geodesic distance within
*        the PolyMap's input space (if "forward" is zero) or output
*        space (if "forward" is non-zero).
*     maxacc
*        The maximum allowed accuracy for an acceptable polynomial,
*        expressed as a geodesic distance within the PolyMap's input
*        space (if "forward" is zero) or output space (if "forward" is
*        non-zero).
*     maxorder
*        The maximum allowed polynomial order. This is one more than the
*        maximum power of either input axis. So for instance, a value of
*        3 refers to a quadratic polynomial. Note, cross terms with total
*        powers greater than or equal to maxorder are not inlcuded in the
*        fit. So the maximum number of terms in each of the fitted polynomials
*        is maxorder*(maxorder+1)/2.
*     lbnd
*        An array holding the lower bounds of a rectangular region within
*        the PolyMap's input space (if "forward" is zero) or output
*        space (if "forward" is non-zero). The new polynomial will be
*        evaluated over this rectangle.
*     ubnd
*        An array holding the upper bounds of a rectangular region within
*        the PolyMap's input space (if "forward" is zero) or output
*        space (if "forward" is non-zero). The new polynomial will be
*        evaluated over this rectangle.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     - Non-zero if a fit was performed succesfully, to at least the
*     maximum allowed accuracy. Zero if the fit failed, or an error
*     occurred.

*  Notes:
*     - No error is reported if the fit fails to achieve the required
*     maximum accuracy.
*     - An error is reported if the transformation that is not being
*     replaced is not defined.
*     - An error is reported if the PolyMap does not have equal numbers
*     of inputs and outputs.
*     - An error is reported if the PolyMap has more than 2 inputs or outputs.

*/

/* Local Variables: */
   double **table;
   double *cofs;
   double racc;
   double scales[ 4 ];
   int ndim;
   int ncof;
   int nsamp;
   int order;
   int result;

/* Check inherited status */
   if( !astOK ) return 0;

/* Check the PolyMap can be used. */
   ndim = astGetNin( this );
   if( astGetNout( this ) != ndim && astOK ) {
      astError( AST__BADNI, "astPolyTran(%s): Supplied %s has "
                "different number of inputs (%d) and outputs (%d).",
                status, astGetClass( this ), astGetClass( this ), ndim,
                astGetNout( this ) );

   } else if( ndim > 2 && astOK ) {
      astError( AST__BADNI, "astPolyTran(%s): Supplied %s has "
                "too many inputs and outputs (%d) - must be 1 or 2.",
                status, astGetClass( this ), astGetClass( this ), ndim );
   }

   if( forward != astGetInvert( this ) ){
      if( ! this->ncoeff_i && astOK ) {
         astError( AST__NODEF, "astPolyTran(%s): Supplied %s has "
                   "no inverse transformation.", status, astGetClass( this ),
                   astGetClass( this ) );
      }
   } else {
      if( ! this->ncoeff_f && astOK ) {
         astError( AST__NODEF, "astPolyTran(%s): Supplied %s has "
                   "no forward transformation.", status, astGetClass( this ),
                   astGetClass( this ) );
      }
   }

/* Check the bounds can be used. */
   if( !lbnd || !ubnd ) {
      astError( AST__NODEF, "astPolyTran(%s): No upper and/or lower bounds "
                "supplied.", status, astGetClass( this ) );
   } else if( lbnd[ 0 ] >= ubnd[ 0 ] && astOK ) {
      astError( AST__NODEF, "astPolyTran(%s): Supplied upper "
                "bound for the first axis (%g) is less than or equal to the "
                "supplied lower bound (%g).", status, astGetClass( this ),
                lbnd[ 0 ], ubnd[ 0 ] );
   } else if( ndim == 2 && lbnd[ 1 ] >= ubnd[ 1 ] && astOK ) {
      astError( AST__NODEF, "astPolyTran(%s): Supplied upper "
                "bound for the second axis (%g) is less than or equal to the "
                "supplied lower bound (%g).", status, astGetClass( this ),
                lbnd[ 1 ], ubnd[ 1 ] );
   }

/* Initialise pointer to work space. */
   table = NULL;

/* Loop over increasing polynomial orders until the required accuracy is
   achieved, up to the specified maximum order. The "order" value is one more
   than the maximum power in the polynomial (so a quadratic has "order" 3). */
   if( maxorder < 2 ) maxorder = 2;
   for( order = 2; order <= maxorder; order++ ) {

/* First do 2D PolyMaps. */
      if( ndim == 2 ) {

/* Sample the requested polynomial transformation at a grid of points. This
   grid covers the user-supplied region, using 2*order points on each
   axis. */
         table = SamplePoly2D( this, !forward, table, lbnd, ubnd, 2*order,
                               &nsamp, scales, status );

/* Fit the polynomial. Always fit a linear polynomial ("order" 2) to any
   dummy second axis. If successfull, replace the PolyMap transformation
   and break out of the order loop. */
         cofs = FitPoly2D( this, forward,  nsamp, acc, order, table, scales,
                           &ncof, &racc, status );

/* Now do 1D PolyMaps. */
      } else {
         table = SamplePoly1D( this, !forward, table, lbnd[ 0 ], ubnd[ 0 ],
                               2*order, &nsamp, scales, status );
         cofs = FitPoly1D( this, forward, nsamp, acc, order, table, scales,
                           &ncof, &racc, status );
      }

/* If the fit was succesful, replace the PolyMap transformation and break
   out of the order loop. */
      if( cofs && ( racc < acc || ( racc < maxacc && order == maxorder ) ) ) {
         StoreArrays( this, forward, ncof, cofs, status );
         break;
      } else {
        cofs = astFree( cofs );
      }
   }

/* If no fit was produced, return zero. */
   result = cofs ? 1 : 0;

/* Free resources. */
   cofs = astFree( cofs );
   table = astFreeDouble( table );

/* Return the result. */
   return result;
}

static double **SamplePoly1D( AstPolyMap *this, int forward, double **table,
                              double lbnd, double ubnd, int npoint, int *nsamp,
                              double scales[2], int *status ){
/*
*  Name:
*     SamplePoly1D

*  Purpose:
*     Create a table of input and output positions for a 1D PolyMap.

*  Type:
*     Private function.

*  Synopsis:
*     double **SamplePoly1D( AstPolyMap *this, int forward, double **table,
*                            double lbnd, double ubnd, int npoint, int *nsamp,
*                            double scales[2], int *status )

*  Description:
*     This function creates a table containing samples of the requested
*     polynomial transformation at a grid of input points. This grid covers
*     the user-supplied region, using "npoint" points.

*  Parameters:
*     this
*        The PolyMap.
*     forward
*        If non-zero, then the forward PolyMap transformation is sampled.
*        Otherwise the inverse transformation is sampled.
*     table
*        Pointer to a previous table created by this function, which is
*        to be re-used, or NULL.
*     lbnd
*        The lower bounds of the region within the PolyMap's 1D input space
*        (if "forward" is non-zero) or output space (if "forward" is zero).
*        The new polynomial will be evaluated over this region.
*     ubnd
*        The upper bounds of the region within the PolyMap's 1D input space
*        (if "forward" is non-zero) or output space (if "forward" is zero).
*        The new polynomial will be evaluated over this region.
*     npoint
*        The number of points to use.
*     nsamp
*        Address of an int in which to return the total number of samples
*        in the returned table.
*     scales
*        Array in which to return the scaling factors for the two columns
*        of the returned table. Multiplying the returned table values by
*        the scale factor produces PolyMap input or output axis values.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*        Pointer to an array of 2 pointers. Each of these pointers points
*        to an array of "nsamp" doubles, being the sampled values for y1
*        and x1 in that order. Here x1 is the input value for the sampled
*        transformation (these are spaced on the regular grid specified
*        by lbnd, ubnd and npoint), and y1 is the output position produced
*        by the sampled transformation. The returned values are scaled so
*        that each column has an RMS value of 1.0. The scaling factors that
*        convert scaled values into original values are returned in "scales".
*        The returned pointer should be freed using astFreeDouble when no
*        longer needed.

*/

/* Local Variables: */
   AstPointSet *ps1;
   AstPointSet *ps2;
   double **result;
   double *p0;
   double *p1;
   double *ptr1[ 1 ];
   double *ptr2[ 1 ];
   double delta0;
   double rms;
   double sum;
   double val0;
   int i;
   int icol;

/* Initialise returned value */
   result = table;
   *nsamp = 0;

/* Check inherited status */
   if( !astOK ) return result;

/* Ensure we have a table of the correct size. */
   *nsamp = npoint;
   if( !result ) result = astCalloc( 2, sizeof( double * ) );
   if( result ) {
      for( i = 0; i < 2; i++ ) {
         result[ i ] = astRealloc( result[ i ] , (*nsamp)*sizeof( double ) );
      }
   }

/* Work out the step sizes for the grid. */
   delta0 = ( ubnd - lbnd )/( npoint - 1 );

/* Create a PointSet to hold the grid of input positions. Use column 1
   of the table to hold the PointSet values. */
   ps1 = astPointSet( *nsamp, 1, " ", status );
   ptr1[ 0 ] = result[ 1 ];
   astSetPoints( ps1, ptr1 );

/* Create a PointSet to hold the grid of output positions. Use column 0
   of the table to hold the PointSet values. */
   ps2 = astPointSet( *nsamp, 1, " ", status );
   ptr2[ 0 ] = result[ 0 ];
   astSetPoints( ps2, ptr2 );
   if( astOK ) {

/* Calculate the grid of input positions and store in the PointSet and
   therefore also in the returned table. Rounding error may cause the
   final point to be just over the upper bound, so we leave the final
   point out of the loop and set it explicitly to the upper bound
   afterwards. */
      val0 = lbnd;
      p0 = ptr1[ 0 ];
      for( i = 0; i < npoint-1; i++ ) {
         *(p0++) = val0;
         val0 += delta0;
      }
      *p0 = ubnd;

/* Transform the input grid to get the output grid. */
      (void) astTransform( this, ps1, forward, ps2 );

/* Scale each column in turn. */
      for( icol = 0; icol < 2; icol++ ) {

/* Find the RMS of the values in the column. */
         sum = 0.0;
         p0 = result[ icol ];
         p1 = p0 + (*nsamp);
         for( ; p0 < p1; p0++ ) sum += ( *p0 )*( *p0 );
         rms = sqrt( sum/(*nsamp) );

/* Divide the table values by the RMS. */
         p0 = result[ icol ];
         p1 = p0 + (*nsamp);
         for( ; p0 < p1; p0++ ) *p0 /= rms;

/* Return the RMS as the scale factor. */
         scales[ icol ] = rms;
      }
   }

/* Free resources */
   ps1 = astAnnul( ps1 );
   ps2 = astAnnul( ps2 );

/* If an error occurred, free the returned array. */
   if( !astOK ) result = astFreeDouble( result );

/* Return a pointer to the table. */
   return result;
}

static double **SamplePoly2D( AstPolyMap *this, int forward, double **table,
                              const double *lbnd, const double *ubnd, int npoint,
                              int *nsamp, double scales[4], int *status ){
/*
*  Name:
*     SamplePoly2D

*  Purpose:
*     Create a table of input and output positions for a 2D PolyMap.

*  Type:
*     Private function.

*  Synopsis:
*     double **SamplePoly2D( AstPolyMap *this, int forward, double **table,
*                            const double *lbnd, const double *ubnd, int npoint,
*                            int *nsamp, double scales[4], int *status )

*  Description:
*     This function creates a table containing samples of the requested
*     polynomial transformation at a grid of input points. This grid covers
*     the user-supplied region, using "npoint" points on each axis.

*  Parameters:
*     this
*        The PolyMap.
*     forward
*        If non-zero, then the forward PolyMap transformation is sampled.
*        Otherwise the inverse transformation is sampled.
*     table
*        Pointer to a previous table created by this function, which is
*        to be re-used, or NULL.
*     lbnd
*        An array holding the lower bounds of a rectangular region within
*        the PolyMap's input space (if "forward" is non-zero) or output
*        space (if "forward" is zero). The new polynomial will be
*        evaluated over this rectangle.
*     ubnd
*        An array holding the upper bounds of a rectangular region within
*        the PolyMap's input space (if "forward" is non-zero) or output
*        space (if "forward" is zero). The new polynomial will be
*        evaluated over this rectangle.
*     npoint
*        The number of points along each edge of the grid.
*     nsamp
*        Address of an int in which to return the total number of samples
*        in the returned table.
*     scales
*        Array in which to return the scaling factors for the four
*        columns of the returned table. Multiplying the returned table
*        values by the scale factor produces PolyMap input or output axis
*        values.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*        Pointer to an array of 4 pointers. Each of these pointers points
*        to an array of "nsamp" doubles, being the sampled values for y1,
*        y2, x1 or x2 in that order. Here (x1,x2) are the input values
*        for the sampled transformation (these are spaced on the regular
*        grid specified by lbnd, ubnd and npoint), and (y1,y2) are the
*        output positions produced by the sampled transformation. The
*        returned values are scaled so that each column has an RMS value
*        of 1.0. The scaling factors that convert scaled values into
*        original values are returned in "scales". The returned pointer
*        should be freed using astFreeDouble when no longer needed.

*/

/* Local Variables: */
   AstPointSet *ps1;
   AstPointSet *ps2;
   double **result;
   double *p0;
   double *p1;
   double *ptr1[ 2 ];
   double *ptr2[ 2 ];
   double delta1;
   double delta0;
   double rms;
   double sum;
   double val0;
   double val1;
   int i;
   int icol;
   int j;

/* Initialise returned value */
   result = table;
   *nsamp = 0;

/* Check inherited status */
   if( !astOK ) return result;

/* Ensure we have a table of the correct size. */
   *nsamp = npoint*npoint;
   if( !result ) result = astCalloc( 4, sizeof( double * ) );
   if( result ) {
      for( i = 0; i < 4; i++ ) {
         result[ i ] = astRealloc( result[ i ] , (*nsamp)*sizeof( double ) );
      }
   }

/* Work out the step sizes for the grid. */
   delta0 = ( ubnd[ 0 ] - lbnd[ 0 ] )/( npoint - 1 );
   delta1 = ( ubnd[ 1 ] - lbnd[ 1 ] )/( npoint - 1 );

/* Create a PointSet to hold the grid of input positions. Use columns 2
   and 3 of the table to hold the PointSet values. */
   ps1 = astPointSet( *nsamp, 2, " ", status );
   ptr1[ 0 ] = result[ 2 ];
   ptr1[ 1 ] = result[ 3 ];
   astSetPoints( ps1, ptr1 );

/* Create a PointSet to hold the grid of output positions. Use columns 0
   and 1 of the table to hold the PointSet values. */
   ps2 = astPointSet( *nsamp, 2, " ", status );
   ptr2[ 0 ] = result[ 0 ];
   ptr2[ 1 ] = result[ 1 ];
   astSetPoints( ps2, ptr2 );
   if( astOK ) {

/* Calculate the grid of input positions and store in the PointSet and
   therefore also in the returned table. Rounding error may cause the
   final point on each axis to be just over the upper bound, so we leave
   the final point out of the loop and set it explicitly to the upper bound
   afterwards. */
      val0 = lbnd[ 0 ];
      p0 = ptr1[ 0 ];
      p1 = ptr1[ 1 ];
      for( i = 0; i < npoint - 1; i++ ) {
         val1 = lbnd[ 1 ];
         for( j = 0; j < npoint-1; j++ ) {
             *(p0++) = val0;
             *(p1++) = val1;
             val1 += delta1;
         }
         *(p0++) = val0;
         *(p1++) = ubnd[ 1 ];
         val0 += delta0;
      }


      val1 = lbnd[ 1 ];
      for( j = 0; j < npoint-1; j++ ) {
          *(p0++) = ubnd[ 0 ];
          *(p1++) = val1;
          val1 += delta1;
      }
      *(p0++) = ubnd[ 0 ];
      *(p1++) = ubnd[ 1 ];


/* Transform the input grid to get the output grid. */
      (void) astTransform( this, ps1, forward, ps2 );

/* Scale each pair of columns in turn. Use the same scale factor for
   each axis in order to ensure an isotropic metric. */
      for( icol = 0; icol < 4; icol += 2 ) {

/* Find the RMS of the values in the two columns. */
         sum = 0.0;
         p0 = result[ icol ];
         p1 = p0 + (*nsamp);
         for( ; p0 < p1; p0++ ) sum += ( *p0 )*( *p0 );

         p0 = result[ icol + 1 ];
         p1 = p0 + (*nsamp);
         for( ; p0 < p1; p0++ ) sum += ( *p0 )*( *p0 );

         rms = sqrt( sum/(2*(*nsamp)) );

/* Divide the table values by the RMS. */
         p0 = result[ icol ];
         p1 = p0 + (*nsamp);
         for( ; p0 < p1; p0++ ) *p0 /= rms;

         p0 = result[ icol + 1 ];
         p1 = p0 + (*nsamp);
         for( ; p0 < p1; p0++ ) *p0 /= rms;

/* Return the RMS as the scale factor. */
         scales[ icol ] = rms;
         scales[ icol + 1 ] = rms;
      }
   }

/* Free resources */
   ps1 = astAnnul( ps1 );
   ps2 = astAnnul( ps2 );

/* If an error occurred, free the returned array. */
   if( !astOK ) result = astFreeDouble( result );

/* Return a pointer to the table. */
   return result;
}

static void SetAttrib( AstObject *this_object, const char *setting, int *status ) {
/*
*  Name:
*     SetAttrib

*  Purpose:
*     Set an attribute value for a PolyMap.

*  Type:
*     Private function.

*  Synopsis:
*     #include "polymap.h"
*     void SetAttrib( AstObject *this, const char *setting )

*  Class Membership:
*     PolyMap member function (over-rides the astSetAttrib protected
*     method inherited from the Mapping class).

*  Description:
*     This function assigns an attribute value for a PolyMap, the
*     attribute and its value being specified by means of a string of
*     the form:
*
*        "attribute= value "
*
*     Here, "attribute" specifies the attribute name and should be in
*     lower case with no white space present. The value to the right
*     of the "=" should be a suitable textual representation of the
*     value to be assigned and this will be interpreted according to
*     the attribute's data type.  White space surrounding the value is
*     only significant for string attributes.

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     setting
*        Pointer to a null-terminated string specifying the new attribute
*        value.
*/

/* Local Variables: */
   AstPolyMap *this;             /* Pointer to the PolyMap structure */
   double dval;                  /* Floating point attribute value */
   int ival;                     /* Integer attribute value */
   int len;                      /* Length of setting string */
   int nc;                       /* Number of characters read by astSscanf */

/* Check the global error status. */
   if ( !astOK ) return;

/* Obtain a pointer to the PolyMap structure. */
   this = (AstPolyMap *) this_object;

/* Obtain the length of the setting string. */
   len = (int) strlen( setting );

/* Test for each recognised attribute in turn, using "astSscanf" to parse
   the setting string and extract the attribute value (or an offset to
   it in the case of string values). In each case, use the value set
   in "nc" to check that the entire string was matched. Once a value
   has been obtained, use the appropriate method to set it. */

/* IterInverse. */
/* ------------ */
   if ( nc = 0,
        ( 1 == astSscanf( setting, "iterinverse= %d %n", &ival, &nc ) )
        && ( nc >= len ) ) {
      astSetIterInverse( this, ival );

/* NiterInverse. */
/* ------------- */
   } else if ( nc = 0,
        ( 1 == astSscanf( setting, "niterinverse= %d %n", &ival, &nc ) )
        && ( nc >= len ) ) {
      astSetNiterInverse( this, ival );

/* TolInverse. */
/* ----------- */
   } else if ( nc = 0,
        ( 1 == astSscanf( setting, "tolinverse= %lg %n", &dval, &nc ) )
        && ( nc >= len ) ) {
      astSetTolInverse( this, dval );

/* If the attribute is still not recognised, pass it on to the parent
   method for further interpretation. */
   } else {
      (*parent_setattrib)( this_object, setting, status );
   }
}

static void StoreArrays( AstPolyMap *this, int forward, int ncoeff,
                         const double *coeff, int *status ){
/*
*  Name:
*     StoreArrays

*  Purpose:
*     Store the dynamic arrays for a single transformation within a PolyMap

*  Type:
*     Private function.

*  Synopsis:
*     #include "polymap.h"
*     void StoreArrays( AstPolyMap *this, int forward, int ncoeff,
*                       const double *coeff, int *status )

*  Class Membership:
*     PolyMap initialiser.

*  Description:
*     This function sets up the arrays within a PolyMap structure that
*     describes either the forward or inverse transformation.

*  Parameters:
*     this
*        The PolyMap.
*     forward
*        If non-zero, replace the forward transformation. Otherwise,
*        replace the inverse transformation.
*     ncoeff
*        The number of non-zero coefficients necessary to define the
*        specified transformation of the PolyMap. If zero is supplied, the
*        transformation will be undefined.
*     coeff
*        An array containing "ncof*( 2 + nin )" elements. Each group
*	 of "2 + nin" adjacent elements describe a single coefficient of
*	 the transformation. Within each such group, the first
*	 element is the coefficient value; the next element is the
*	 integer index of the PolyMap output which uses the coefficient
*	 within its defining polynomial (the first output has index 1);
*	 the remaining elements of the group give the integer powers to
*	 use with each input coordinate value (powers must not be
*	 negative).
*     status
*        Pointer to inherited status.
*/

/* Local Variables: */
   const double *group;          /* Pointer to start of next coeff. description */
   int *pows;                    /* Pointer to powers for current coeff. */
   int gsize;                    /* Length of each coeff. description */
   int i;                        /* Loop count */
   int ico;                      /* Index of next coeff. for current input or output */
   int iin;                      /* Input index extracted from coeff. description */
   int iout;                     /* Output index extracted from coeff. description */
   int j;                        /* Loop count */
   int nin;                      /* Number of orignal inputs */
   int nout;                     /* Number of original outputs */
   int pow;                      /* Power extracted from coeff. description */

/* Check the global status. */
   if ( !astOK ) return;

/* First Free any existing arrays. */
   FreeArrays( this, forward, status );

/* Get the number of inputs and outputs of the original uninverted PolyMap.
   Swap the transformation to be set if the PolyMap has been inverted. */
   if( astGetInvert( this ) ) {
      nout = astGetNin( this );
      nin = astGetNout( this );
      forward = !forward;
   } else {
      nin = astGetNin( this );
      nout = astGetNout( this );
   }

/* Now initialise the original forward transformation, if required. */
   if( forward && ncoeff > 0 ) {

/* Create the arrays decribing the forward transformation. */
      this->ncoeff_f = astMalloc( sizeof( int )*(size_t) nout );
      this->mxpow_f = astMalloc( sizeof( int )*(size_t) nin );
      this->power_f = astMalloc( sizeof( int ** )*(size_t) nout );
      this->coeff_f = astMalloc( sizeof( double * )*(size_t) nout );
      if( astOK ) {

/* Initialise the count of coefficients for each output coordinate to zero. */
         for( i = 0; i < nout; i++ ) this->ncoeff_f[ i ] = 0;

/* Initialise max power for each input coordinate to zero. */
         for( j = 0; j < nin; j++ ) this->mxpow_f[ j ] = 0;

/* Scan through the supplied forward coefficient array, counting the
   number of coefficients which relate to each output. Also find the
   highest power used for each input axis. Report errors if any unusable
   values are found in the supplied array. */
         group = coeff;
         gsize = 2 + nin;
         for( i = 0; i < ncoeff && astOK; i++, group += gsize ) {

            iout = floor( group[ 1 ] + 0.5 );
            if( iout < 1 || iout > nout ) {
               astError( AST__BADCI, "astInitPolyMap(%s): Forward "
                         "coefficient %d referred to an illegal output "
                         "coordinate %d.", status, astGetClass( this ), i + 1,
                         iout );
               astError( AST__BADCI, "This number should be in the "
                         "range 1 to %d.", status, nout );
               break;
            }

            this->ncoeff_f[ iout - 1 ]++;

            for( j = 0; j < nin; j++ ) {
               pow = floor( group[ 2 + j ] + 0.5 );
               if( pow < 0 ) {
                  astError( AST__BADPW, "astInitPolyMap(%s): Forward "
                            "coefficient %d has a negative power (%d) "
                            "for input coordinate %d.", status,
                            astGetClass( this ), i + 1, pow, j + 1 );
                  astError( AST__BADPW, "All powers should be zero or "
                            "positive." , status);
                  break;
               }
               if( pow > this->mxpow_f[ j ] ) this->mxpow_f[ j ] = pow;
            }
         }

/* Allocate the arrays to store the input powers associated with each
   coefficient, and the coefficient values. Reset the coefficient count
   for each axis to zero afterwards so that we can use the array as an index
   to the next vacant slot within the following loop. */
         for( i = 0; i < nout; i++ ) {
            this->power_f[ i ] = astMalloc( sizeof( int * )*
                                           (size_t) this->ncoeff_f[ i ] );
            this->coeff_f[ i ] = astMalloc( sizeof( double )*
                                           (size_t) this->ncoeff_f[ i ] );
            this->ncoeff_f[ i ] = 0;
         }

         if( astOK ) {

/* Extract the coefficient values and powers form the supplied array and
   store them in the arrays created above. */
            group = coeff;
            for( i = 0; i < ncoeff && astOK; i++, group += gsize ) {
               iout = floor( group[ 1 ] + 0.5 ) - 1;
               ico = ( this->ncoeff_f[ iout ] )++;
               this->coeff_f[ iout ][ ico ] = group[ 0 ];

               pows = astMalloc( sizeof( int )*(size_t) nin );
               this->power_f[ iout ][ ico ] = pows;
               if( astOK ) {
                  for( j = 0; j < nin; j++ ) {
                     pows[ j ] = floor( group[ 2 + j ] + 0.5 );
                  }
               }
            }
         }
      }
   }

/* Now initialise the original inverse transformation, if required. */
   if( !forward && ncoeff > 0 ) {

/* Create the arrays decribing the inverse transformation. */
      this->ncoeff_i = astMalloc( sizeof( int )*(size_t) nin );
      this->mxpow_i = astMalloc( sizeof( int )*(size_t) nout );
      this->power_i = astMalloc( sizeof( int ** )*(size_t) nin );
      this->coeff_i = astMalloc( sizeof( double * )*(size_t) nin );
      if( astOK ) {

/* Initialise the count of coefficients for each input coordinate to zero. */
         for( i = 0; i < nin; i++ ) this->ncoeff_i[ i ] = 0;

/* Initialise max power for each output coordinate to zero. */
         for( j = 0; j < nout; j++ ) this->mxpow_i[ j ] = 0;

/* Scan through the supplied inverse coefficient array, counting the
   number of coefficients which relate to each input. Also find the
   highest power used for each output axis. Report errors if any unusable
   values are found in the supplied array. */
         group = coeff;

         gsize = 2 + nout;
         for( i = 0; i < ncoeff && astOK; i++, group += gsize ) {

            iin = floor( group[ 1 ] + 0.5 );
            if( iin < 1 || iin > nin ) {
               astError( AST__BADCI, "astInitPolyMap(%s): Inverse "
                         "coefficient %d referred to an illegal input "
                         "coordinate %d.", status, astGetClass( this ),
                         i + 1, iin );
               astError( AST__BADCI, "This number should be in the "
                         "range 1 to %d.", status, nin );
               break;
            }

            this->ncoeff_i[ iin - 1 ]++;

            for( j = 0; j < nout; j++ ) {
               pow = floor( group[ 2 + j ] + 0.5 );
               if( pow < 0 ) {
                  astError( AST__BADPW, "astInitPolyMap(%s): Inverse "
                            "coefficient %d has a negative power (%d) "
                            "for output coordinate %d.", status,
                            astGetClass( this ), i + 1, pow, j + 1 );
                  astError( AST__BADPW, "All powers should be zero or "
                            "positive." , status);
                  break;
               }
               if( pow > this->mxpow_i[ j ] ) this->mxpow_i[ j ] = pow;
            }
         }

/* Allocate the arrays to store the output powers associated with each
   coefficient, and the coefficient values. Reset the coefficient count
   for each axis to zero afterwards so that we can use the array as an index
   to the next vacant slot within the following loop. */
         for( i = 0; i < nin; i++ ) {
            this->power_i[ i ] = astMalloc( sizeof( int * )*
                                           (size_t) this->ncoeff_i[ i ] );
            this->coeff_i[ i ] = astMalloc( sizeof( double )*
                                           (size_t) this->ncoeff_i[ i ] );
            this->ncoeff_i[ i ] = 0;
         }

         if( astOK ) {

/* Extract the coefficient values and powers form the supplied array and
   store them in the arrays created above. */
            group = coeff;
            for( i = 0; i < ncoeff && astOK; i++, group += gsize ) {
               iin = floor( group[ 1 ] + 0.5 ) - 1;
               ico = ( this->ncoeff_i[ iin ] )++;
               this->coeff_i[ iin ][ ico ] = group[ 0 ];

               pows = astMalloc( sizeof( int )*(size_t) nout );
               this->power_i[ iin ][ ico ] = pows;
               if( astOK ) {
                  for( j = 0; j < nout; j++ ) {
                     pows[ j ] = floor( group[ 2 + j ] + 0.5 );
                  }
               }
            }
         }
      }
   }
}

static int TestAttrib( AstObject *this_object, const char *attrib, int *status ) {
/*
*  Name:
*     TestAttrib

*  Purpose:
*     Test if a specified attribute value is set for a PolyMap.

*  Type:
*     Private function.

*  Synopsis:
*     #include "polymap.h"
*     int TestAttrib( AstObject *this, const char *attrib, int *status )

*  Class Membership:
*     PolyMap member function (over-rides the astTestAttrib protected
*     method inherited from the Mapping class).

*  Description:
*     This function returns a boolean result (0 or 1) to indicate whether
*     a value has been set for one of a PolyMap's attributes.

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     attrib
*        Pointer to a null-terminated string specifying the attribute
*        name.  This should be in lower case with no surrounding white
*        space.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     One if a value has been set, otherwise zero.

*  Notes:
*     - A value of zero will be returned if this function is invoked
*     with the global status set, or if it should fail for any reason.
*/

/* Local Variables: */
   AstPolyMap *this;             /* Pointer to the PolyMap structure */
   int result;                   /* Result value to return */

/* Initialise. */
   result = 0;

/* Check the global error status. */
   if ( !astOK ) return result;

/* Obtain a pointer to the PolyMap structure. */
   this = (AstPolyMap *) this_object;

/* Check the attribute name and test the appropriate attribute. */

/* IterInverse. */
/* ------------ */
   if ( !strcmp( attrib, "iterinverse" ) ) {
      result = astTestIterInverse( this );

/* NiterInverse. */
/* ------------- */
   } else if ( !strcmp( attrib, "niterinverse" ) ) {
      result = astTestNiterInverse( this );

/* TolInverse. */
/* ----------- */
   } else if ( !strcmp( attrib, "tolinverse" ) ) {
      result = astTestTolInverse( this );

/* If the attribute is still not recognised, pass it on to the parent
   method for further interpretation. */
   } else {
      result = (*parent_testattrib)( this_object, attrib, status );
   }

/* Return the result, */
   return result;
}

static AstPointSet *Transform( AstMapping *this, AstPointSet *in,
                               int forward, AstPointSet *out, int *status ) {
/*
*  Name:
*     Transform

*  Purpose:
*     Apply a PolyMap to transform a set of points.

*  Type:
*     Private function.

*  Synopsis:
*     #include "polymap.h"
*     AstPointSet *Transform( AstMapping *this, AstPointSet *in,
*                             int forward, AstPointSet *out, int *status )

*  Class Membership:
*     PolyMap member function (over-rides the astTransform protected
*     method inherited from the Mapping class).

*  Description:
*     This function takes a PolyMap and a set of points encapsulated in a
*     PointSet and transforms the points.

*  Parameters:
*     this
*        Pointer to the PolyMap.
*     in
*        Pointer to the PointSet holding the input coordinate data.
*     forward
*        A non-zero value indicates that the forward coordinate transformation
*        should be applied, while a zero value requests the inverse
*        transformation.
*     out
*        Pointer to a PointSet which will hold the transformed (output)
*        coordinate values. A NULL value may also be given, in which case a
*        new PointSet will be created by this function.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     Pointer to the output (possibly new) PointSet.

*  Notes:
*     -  A null pointer will be returned if this function is invoked with the
*     global error status set, or if it should fail for any reason.
*     -  The number of coordinate values per point in the input PointSet must
*     match the number of columns in the PolyMap being applied.
*     -  The number of coordinate values per point in the output PointSet will
*     equal the number of rows in the PolyMap being applied.
*     -  If an output PointSet is supplied, it must have space for sufficient
*     number of points and coordinate values per point to accommodate the
*     result. Any excess space will be ignored.
*/

/* Local Variables: */
   AstPointSet *result;          /* Pointer to output PointSet */
   AstPolyMap *map;              /* Pointer to PolyMap to be applied */
   double **coeff;               /* Pointer to coefficient value arrays */
   double **ptr_in;              /* Pointer to input coordinate data */
   double **ptr_out;             /* Pointer to output coordinate data */
   double **work;                /* Pointer to exponentiated axis values */
   double *outcof;               /* Pointer to next coefficient value */
   double outval;                /* Output axis value */
   double term;                  /* Term to be added to output value */
   double xp;                    /* Exponentiated input axis value */
   int ***power;                 /* Pointer to coefficient power arrays */
   int **outpow;                 /* Pointer to next set of axis powers */
   int *mxpow;                   /* Pointer to max used power for each input */
   int *ncoeff;                  /* Pointer to no. of coefficients */
   int in_coord;                 /* Index of output coordinate */
   int ico;                      /* Coefficient index */
   int nc;                       /* No. of coefficients in polynomial */
   int ncoord_in;                /* Number of coordinates per input point */
   int ncoord_out;               /* Number of coordinates per output point */
   int npoint;                   /* Number of points */
   int out_coord;                /* Index of output coordinate */
   int point;                    /* Loop counter for points */
   int pow;                      /* Next axis power */

/* Check the global error status. */
   if ( !astOK ) return NULL;

/* Obtain a pointer to the PolyMap. */
   map = (AstPolyMap *) this;

/* Apply the parent mapping using the stored pointer to the Transform member
   function inherited from the parent Mapping class. This function validates
   all arguments and generates an output PointSet if necessary, but does not
   actually transform any coordinate values. */
   result = (*parent_transform)( this, in, forward, out, status );

/* Determine whether to apply the original forward or inverse mapping,
   according to the direction specified and whether the mapping has been
   inverted. */
   if ( astGetInvert( map ) ) forward = !forward;

/* We will now extend the parent astTransform method by performing the
   calculations needed to generate the output coordinate values. */

/* If we are using the original inverse transformatiom, and the IterInverse
   attribute is non-zero, use an iterative inverse algorithm rather than any
   inverse transformation defined within the PolyMap. */
   if( !forward && astGetIterInverse(map) ) {
      IterInverse( map, in, result, status );

/* Otherwise, determine the numbers of points and coordinates per point from
   the input and output PointSets and obtain pointers for accessing the input
   and output coordinate values. */
   } else {
      ncoord_in = astGetNcoord( in );
      ncoord_out = astGetNcoord( result );
      npoint = astGetNpoint( in );
      ptr_in = astGetPoints( in );
      ptr_out = astGetPoints( result );

/* Get a pointer to the arrays holding the required coefficient
   values and powers, according to the direction of mapping required. */
      if ( forward ) {
         ncoeff = map->ncoeff_f;
         coeff = map->coeff_f;
         power = map->power_f;
         mxpow = map->mxpow_f;
      } else {
         ncoeff = map->ncoeff_i;
         coeff = map->coeff_i;
         power = map->power_i;
         mxpow = map->mxpow_i;
      }

/* Allocate memory to hold the required powers of the input axis values. */
      work = astMalloc( sizeof( double * )*(size_t) ncoord_in );
      for( in_coord = 0; in_coord < ncoord_in; in_coord++ ) {
         work[ in_coord ] = astMalloc( sizeof( double )*
                           (size_t) ( astMAX( 2, mxpow[in_coord]+1 ) ) );
      }

/* Perform coordinate arithmetic. */
/* ------------------------------ */
      if ( astOK ) {

/* Loop to apply the polynomial to each point in turn.*/
         for ( point = 0; point < npoint; point++ ) {

/* Find the required powers of the input axis values and store them
   in the work array. Note, using a virtual method here slows the PolyMap
   Transform function down by about 5%, compared to doing the equivalent
   calculations in-line. But we need some way to allow the ChebyMap class
   to over-ride the calculation of the powers, so we must do something
   like this. If the 5% slow-down is too much, it can be reduced down to
   about 2% by replacing the invocation of the astPolyPowers_ interface
   function with a direct call to the implementation function itself.
   This involves replacing the astPolyPowers call below with this:

   (**astMEMBER(this,PolyMap,PolyPowers))( (AstPolyMap *) this, work, ncoord_in,
                                           mxpow, ptr_in, point, forward, status );

   The above could be wrapped up in an alternative implementation of the
   astPolyPowers macro, so that it looks the same as the existing code.
   In fact, this scheme could be more widely used to speed up invocation
   of virtual functions within AST. The disadvantage is that the interface
   functions for some virtual methods includes some extra processing,
   over and above simply invoking the implementation function.

   Of course the other way to get rid of the 5% slow down, is to
   revert to using in-line code below, and then replicate this entire
   function in the ChebyMap class, making suitable changes to use
   Chebyshev functions in place of simple powers. But that is bad
   structuring... */
            astPolyPowers( this, work, ncoord_in, mxpow, ptr_in, point,
                           forward );

/* Loop round each output. */
            for( out_coord = 0; out_coord < ncoord_out; out_coord++ ) {

/* Initialise the output value. */
               outval = 0.0;

/* Get pointers to the coefficients and powers for this output. */
               outcof = coeff[ out_coord ];
               outpow = power[ out_coord ];

/* Loop round all polynomial coefficients.*/
               nc = ncoeff[ out_coord ];
               for ( ico = 0; ico < nc && outval != AST__BAD;
                     ico++, outcof++, outpow++ ) {

/* Initialise the current term to be equal to the value of the coefficient.
   If it is bad, store a bad output value. */
                  term = *outcof;
                  if( term == AST__BAD ) {
                     outval = AST__BAD;

/* Otherwise, loop round all inputs */
                  } else {
                     for( in_coord = 0; in_coord < ncoord_in; in_coord++ ) {

/* Get the power of the current input axis value used by the current
   coefficient. If it is zero, pass on. */
                        pow = (*outpow)[ in_coord ];
                        if( pow > 0 ) {

/* Get the axis value raised to the appropriate power. */
                           xp = work[ in_coord ][ pow ];

/* If bad, set the output value bad and break. */
                           if( xp == AST__BAD ) {
                              outval = AST__BAD;
                              break;

/* Otherwise multiply the current term by the exponentiated axis value. */
                           } else {
                              term *= xp;
                           }
                        }
                     }
                  }

/* Increment the output value by the current term of the polynomial. */
                  if( outval != AST__BAD ) outval += term;

               }

/* Store the output value. */
               ptr_out[ out_coord ][ point ] = outval;

            }
         }
      }

/* Free resources. */
      for( in_coord = 0; in_coord < ncoord_in; in_coord++ ) {
         work[ in_coord ] = astFree( work[ in_coord ] );
      }
      work = astFree( work );
   }

/* Return a pointer to the output PointSet. */
   return result;
}

/* Functions which access class attributes. */
/* ---------------------------------------- */
/* Implement member functions to access the attributes associated with
   this class using the macros defined for this purpose in the
   "object.h" file. For a description of each attribute, see the class
   interface (in the associated .h file). */

/*
*att++
*  Name:
*     IterInverse

*  Purpose:
*     Provide an iterative inverse transformation?

*  Type:
*     Public attribute.

*  Synopsis:
*     Integer (boolean).

*  Description:
*     This attribute indicates whether the original inverse transformation
*     of the PolyMap should be implemented via an iterative Newton-Raphson
*     approximation that uses the forward transformation to transform
*     candidate input positions until an output position is found which
*     is close to the required output position. By default, an iterative
*     inverse is provided if, and only if, no inverse polynomial was supplied
*     when the PolyMap was constructed.
*
*     Note, the term "inverse transformation" here refers to the inverse
*     transformation of the original PolyMap, ignoring any subsequent
*     inversions. Also, "input" and "output" refer to the inputs and
*     outputs of the original PolyMap.

*     The NiterInverse and TolInverse attributes provide parameters that
*     control the behaviour of the inverse approximation method.
*
*     The iterative inverse returns AST__BAD axis values at positions
*     for which the inverse transformation is undefined. For instance,
*     if the forward transformation is y = x*x, the iterative inverse
*     will return x = AST__BAD at y = -1. If the inverse transformation
*     is multiply defined, the position returned by the iterative inverse
*     will be the position of the solution that is closest to the
*     supplied position. For instance, using the above example, y = x*x,
*     the iterative inverse will return x = +2 at y = 4, because x = +2
*     is the closest solution to 4 (the other solution is x = -2).

*  Applicability:
*     PolyMap
*        All PolyMaps have this attribute.
*     ChebyMap
*        The ChebyMap class does not currently provide an option for an
*        iterative inverse, and so the IterInverse value is always zero.
*        Setting or clearing the IterInverse attribute of a ChebyMap has
*        no effect.

*  Notes:
*     - The transformation replaced by the iterative algorithm is the
*     transformation from the original PolyMap output space to the
*     original PolyMap input space (i.e. the input and output spaces as
*     defined by the arguments of the PolyMap constructor). This is still
*     the case even if the PolyMap has subsequently been inverted. In
*     other words if a PolyMap is created and then inverted, setting
*     the IterInverse to a non-zero value will replace the forward
*     transformation of the inverted PolyMap (i.e. the inverse
*     transformation of the original PolyMap). It is not possible to
*     replace the other transformation (i.e. from the original PolyMap
*     input space to the original PolyMap output space) with an iterative
*     algorithm.
*     - If a PolyMap that has an iterative inverse transformation is
*     subsequently inverted, the inverted PolyMap will have an iterative
*     forward transformation.
*     - An iterative inverse can only be used if the PolyMap has equal
*     numbers of inputs and outputs, as given by the Nin and Nout
*     attributes. An error will be reported if IterInverse is set non-zero
*     for a PolyMap that does not meet this requirement.

*att--
*/
astMAKE_CLEAR(PolyMap,IterInverse,iterinverse,-INT_MAX)
astMAKE_GET(PolyMap,IterInverse,int,0,( ( this->iterinverse == -INT_MAX ) ?
                                          (this->ncoeff_i == 0) : this->iterinverse ))
astMAKE_SET(PolyMap,IterInverse,int,iterinverse,
  (((astGetNin(this)==astGetNout(this))||!value)?((value?1:0)):(astError(AST__ATTIN,"astSetIterInverse(%s):"
  "Cannot use an iterative inverse because the %s has unequal numbers of "
  "inputs and outputs.", status, astGetClass(this),astGetClass(this)),this->iterinverse)))
astMAKE_TEST(PolyMap,IterInverse,( this->iterinverse != -INT_MAX ))

/* NiterInverse. */
/* --------- */
/*
*att++
*  Name:
*     NiterInverse

*  Purpose:
*     Maximum number of iterations for the iterative inverse transformation.

*  Type:
*     Public attribute.

*  Synopsis:
*     Integer.

*  Description:
*     This attribute controls the iterative inverse transformation
*     used if the IterInverse attribute is non-zero.
*
*     Its value gives the maximum number of iterations of the
*     Newton-Raphson algorithm to be used for each transformed position.
*     The default value is 4. See also attribute TolInverse.

*  Applicability:
*     PolyMap
*        All PolyMaps have this attribute.

*att--
*/
astMAKE_CLEAR(PolyMap,NiterInverse,niterinverse,-INT_MAX)
astMAKE_GET(PolyMap,NiterInverse,int,0,( this->niterinverse == -INT_MAX ? 4 : this->niterinverse))
astMAKE_SET(PolyMap,NiterInverse,int,niterinverse,value)
astMAKE_TEST(PolyMap,NiterInverse,( this->niterinverse != -INT_MAX ))

/* TolInverse. */
/* ----------- */
/*
*att++
*  Name:
*     TolInverse

*  Purpose:
*     Target relative error for the iterative inverse transformation.

*  Type:
*     Public attribute.

*  Synopsis:
*     Floating point.

*  Description:
*     This attribute controls the iterative inverse transformation
*     used if the IterInverse attribute is non-zero.
*
*     Its value gives the target relative error in the axis values of
*     each transformed position. Further iterations will be performed
*     until the target relative error is reached, or the maximum number
*     of iterations given by attribute NiterInverse is reached.

*     The default value is 1.0E-6.

*  Applicability:
*     PolyMap
*        All PolyMaps have this attribute.
*att--
*/
astMAKE_CLEAR(PolyMap,TolInverse,tolinverse,AST__BAD)
astMAKE_GET(PolyMap,TolInverse,double,0.0,( this->tolinverse == AST__BAD ? 1.0E-6 : this->tolinverse))
astMAKE_SET(PolyMap,TolInverse,double,tolinverse,value)
astMAKE_TEST(PolyMap,TolInverse,( this->tolinverse != AST__BAD ))

/* Copy constructor. */
/* ----------------- */
static void Copy( const AstObject *objin, AstObject *objout, int *status ) {
/*
*  Name:
*     Copy

*  Purpose:
*     Copy constructor for PolyMap objects.

*  Type:
*     Private function.

*  Synopsis:
*     void Copy( const AstObject *objin, AstObject *objout, int *status )

*  Description:
*     This function implements the copy constructor for PolyMap objects.

*  Parameters:
*     objin
*        Pointer to the object to be copied.
*     objout
*        Pointer to the object being constructed.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     void

*  Notes:
*     -  This constructor makes a deep copy, including a copy of the
*     coefficients associated with the input PolyMap.
*/


/* Local Variables: */
   AstPolyMap *in;               /* Pointer to input PolyMap */
   AstPolyMap *out;              /* Pointer to output PolyMap */
   int nin;                      /* No. of input coordinates */
   int nout;                     /* No. of output coordinates */
   int i;                        /* Loop count */
   int j;                        /* Loop count */

/* Check the global error status. */
   if ( !astOK ) return;

/* Obtain pointers to the input and output PolyMaps. */
   in = (AstPolyMap *) objin;
   out = (AstPolyMap *) objout;

/* Nullify the pointers stored in the output object since these will
   currently be pointing at the input data (since the output is a simple
   byte-for-byte copy of the input). Otherwise, the input data could be
   freed by accidient if the output object is deleted due to an error
   occuring in this function. */
   out->ncoeff_f = NULL;
   out->power_f = NULL;
   out->coeff_f = NULL;
   out->mxpow_f = NULL;

   out->ncoeff_i = NULL;
   out->power_i = NULL;
   out->coeff_i = NULL;
   out->mxpow_i = NULL;

   out->jacobian = NULL;
   out->lintrunc = NULL;

/* Get the number of inputs and outputs of the uninverted Mapping. */
   nin = ( (AstMapping *) in )->nin;
   nout = ( (AstMapping *) in )->nout;

/* Copy the number of coefficients associated with each output of the
   forward transformation of the uninverted Mapping. */
   if( in->ncoeff_f ) {
      out->ncoeff_f = (int *) astStore( NULL, (void *) in->ncoeff_f,
                                        sizeof( int )*(size_t) nout );

/* Copy the maximum power of each input axis value used by the forward
   transformation. */
      out->mxpow_f = (int *) astStore( NULL, (void *) in->mxpow_f,
                                       sizeof( int )*(size_t) nin );

/* Copy the coefficient values used by the forward transformation. */
      if( in->coeff_f ) {
         out->coeff_f = astMalloc( sizeof( double * )*(size_t) nout );
         if( astOK ) {
            for( i = 0; i < nout; i++ ) {
               out->coeff_f[ i ] = (double *) astStore( NULL, (void *) in->coeff_f[ i ],
                                                     sizeof( double )*(size_t) in->ncoeff_f[ i ] );
            }
         }
      }

/* Copy the input axis powers associated with each coefficient of the forward
   transformation. */
      if( in->power_f ) {
         out->power_f = astMalloc( sizeof( int ** )*(size_t) nout );
         if( astOK ) {
            for( i = 0; i < nout; i++ ) {
               out->power_f[ i ] = astMalloc( sizeof( int * )*(size_t) in->ncoeff_f[ i ] );
               if( astOK ) {
                  for( j = 0; j < in->ncoeff_f[ i ]; j++ ) {
                     out->power_f[ i ][ j ] = (int *) astStore( NULL, (void *) in->power_f[ i ][ j ],
                                                                sizeof( int )*(size_t) nin );
                  }
               }
            }
         }
      }
   }

/* Do the same for the inverse transformation. */
   if( in->ncoeff_i ) {
      out->ncoeff_i = (int *) astStore( NULL, (void *) in->ncoeff_i,
                                        sizeof( int )*(size_t) nin );

      out->mxpow_i = (int *) astStore( NULL, (void *) in->mxpow_i,
                                       sizeof( int )*(size_t) nout );

      if( in->coeff_i ) {
         out->coeff_i = astMalloc( sizeof( double * )*(size_t) nin );
         if( astOK ) {
            for( i = 0; i < nin; i++ ) {
               out->coeff_i[ i ] = (double *) astStore( NULL, (void *) in->coeff_i[ i ],
                                                     sizeof( double )*(size_t) in->ncoeff_i[ i ] );
            }
         }
      }

      if( in->power_i ) {
         out->power_i = astMalloc( sizeof( int ** )*(size_t) nin );
         if( astOK ) {
            for( i = 0; i < nin; i++ ) {
               out->power_i[ i ] = astMalloc( sizeof( int * )*(size_t) in->ncoeff_i[ i ] );
               if( astOK ) {
                  for( j = 0; j < in->ncoeff_i[ i ]; j++ ) {
                     out->power_i[ i ][ j ] = (int *) astStore( NULL, (void *) in->power_i[ i ][ j ],
                                                                sizeof( int )*(size_t) nout );
                  }
               }
            }
         }
      }
   }

/* Copy the linear truncation of the PolyMap - if it has been found. */
   if( in->lintrunc ) out->lintrunc = astCopy( in->lintrunc );

/* If an error has occurred, free all the resources allocated above. */
   if( !astOK ) {
      FreeArrays( out, 1, status );
      FreeArrays( out, 0, status );
   }

   return;

}

/* Destructor. */
/* ----------- */
static void Delete( AstObject *obj, int *status ) {
/*
*  Name:
*     Delete

*  Purpose:
*     Destructor for PolyMap objects.

*  Type:
*     Private function.

*  Synopsis:
*     void Delete( AstObject *obj, int *status )

*  Description:
*     This function implements the destructor for PolyMap objects.

*  Parameters:
*     obj
*        Pointer to the object to be deleted.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     void

*  Notes:
*     This function attempts to execute even if the global error status is
*     set.
*/

/* Local Variables: */
   AstPolyMap *this;
   int nc;
   int ic;
   int lstat;
   int error;

/* Obtain a pointer to the PolyMap structure. */
   this = (AstPolyMap *) obj;

/* Free the arrays. */
   FreeArrays( this, 1, status );
   FreeArrays( this, 0, status );

/* Free the resources used to store the Jacobian of the forward
   transformation. */
   if( this->jacobian ) {

/* Get the number of PolyMap inputs. We need to clear any error status
   first since astGetNin returns zero if an error has occurred. The
   Jacobian will only be non-NULL if the number of inputs and outputs
   are equal. */
      error = !astOK;
      if( error ) {
         lstat = astStatus;
         astClearStatus;
      }
      nc = astGetNin( this );
      if( error ) astSetStatus( lstat );

      for( ic = 0; ic < nc; ic++ ) {
         (this->jacobian)[ ic ] = astAnnul( (this->jacobian)[ ic ] );
      }
      this->jacobian = astFree( this->jacobian );
   }

/* Free the linear truncation. */
   if( this->lintrunc ) this->lintrunc = astAnnul( this->lintrunc );
}

/* Dump function. */
/* -------------- */
static void Dump( AstObject *this_object, AstChannel *channel, int *status ) {
/*
*  Name:
*     Dump

*  Purpose:
*     Dump function for PolyMap objects.

*  Type:
*     Private function.

*  Synopsis:
*     void Dump( AstObject *this, AstChannel *channel, int *status )

*  Description:
*     This function implements the Dump function which writes out data
*     for the PolyMap class to an output Channel.

*  Parameters:
*     this
*        Pointer to the PolyMap whose data are being written.
*     channel
*        Pointer to the Channel to which the data are being written.
*     status
*        Pointer to the inherited status variable.
*/

#define KEY_LEN 50               /* Maximum length of a keyword */

/* Local Variables: */
   AstPolyMap *this;             /* Pointer to the PolyMap structure */
   char buff[ KEY_LEN + 1 ];     /* Buffer for keyword string */
   char comm[ 100 ];             /* Buffer for comment string */
   double dval;                  /* Floating point attribute value */
   int i;                        /* Loop index */
   int iv;                       /* Vectorised keyword index */
   int ival;                     /* Integer value */
   int j;                        /* Loop index */
   int k;                        /* Loop index */
   int nin;                      /* No. of input coords */
   int nout;                     /* No. of output coords */
   int set;                      /* Attribute value set? */

/* Check the global error status. */
   if ( !astOK ) return;

/* Obtain a pointer to the PolyMap structure. */
   this = (AstPolyMap *) this_object;

/* Find the number of inputs and outputs of the uninverted Mapping. */
   nin = ( (AstMapping *) this )->nin;
   nout = ( (AstMapping *) this )->nout;

/* Write out values representing the instance variables for the
   PolyMap class.  */

/* First do the forward transformation arrays. Check they are used. */
   if( this->ncoeff_f ) {

/* Store the maximum power of each input axis value used by the forward
   transformation. */
      for( i = 0; i < nin; i++ ){
         (void) sprintf( buff, "MPF%d", i + 1 );
         (void) sprintf( comm, "Max. power of input %d in any forward polynomial", i + 1 );
          astWriteInt( channel, buff, 1, 1, (this->mxpow_f)[ i ], comm );
      }

/* Store the number of coefficients associated with each output of the forward
   transformation. */
      for( i = 0; i < nout; i++ ){
         (void) sprintf( buff, "NCF%d", i + 1 );
         (void) sprintf( comm, "No. of coeff.s for forward polynomial %d", i + 1 );
         astWriteInt( channel, buff, 1, 1, (this->ncoeff_f)[ i ], comm );
      }

/* Store the coefficient values used by the forward transformation. */
      iv = 1;
      for( i = 0; i < nout; i++ ){
         for( j = 0; j < this->ncoeff_f[ i ]; j++, iv++ ){
            if( (this->coeff_f)[ i ][ j ] != AST__BAD ) {
               (void) sprintf( buff, "CF%d", iv );
               (void) sprintf( comm, "Coeff %d of forward polynomial %d", j + 1, i + 1 );
               astWriteDouble( channel, buff, 1, 1, (this->coeff_f)[ i ][ j ], comm );
            }
         }
      }

/* Store the input axis powers associated with each coefficient of the forward
   transformation. */
      iv = 1;
      for( i = 0; i < nout; i++ ){
         for( j = 0; j < this->ncoeff_f[ i ]; j++ ){
            for( k = 0; k < nin; k++, iv++ ){
               if( (this->power_f)[ i ][ j ][ k ] > 0 ) {
                  (void) sprintf( buff, "PF%d", iv );
                  (void) sprintf( comm, "Power of i/p %d for coeff %d of fwd poly %d", k + 1, j + 1, i + 1 );
                  astWriteDouble( channel, buff, 1, 1, (this->power_f)[ i ][ j ][ k ], comm );
               }
            }
         }
      }
   }

/* Now do the inverse transformation arrays. Check they are used. */
   if( this->ncoeff_i ) {

/* Store the maximum power of each output axis value used by the inverse
   transformation. */
      for( i = 0; i < nout; i++ ){
         (void) sprintf( buff, "MPI%d", i + 1 );
         (void) sprintf( comm, "Max. power of output %d in any inverse polynomial", i + 1 );
          astWriteInt( channel, buff, 1, 1, (this->mxpow_i)[ i ], comm );
      }

/* Store the number of coefficients associated with each input of the inverse
   transformation. */
      for( i = 0; i < nin; i++ ){
         (void) sprintf( buff, "NCI%d", i + 1 );
         (void) sprintf( comm, "No. of coeff.s for inverse polynomial %d", i + 1 );
         astWriteInt( channel, buff, 1, 1, (this->ncoeff_i)[ i ], comm );
      }

/* Store the coefficient values used by the inverse transformation. */
      iv = 1;
      for( i = 0; i < nin; i++ ){
         for( j = 0; j < this->ncoeff_i[ i ]; j++, iv++ ){
            if( (this->coeff_i)[ i ][ j ] != AST__BAD ) {
               (void) sprintf( buff, "CI%d", iv );
               (void) sprintf( comm, "Coeff %d of inverse polynomial %d", j + 1, i + 1 );
               astWriteDouble( channel, buff, 1, 1, (this->coeff_i)[ i ][ j ], comm );
            }
         }
      }

/* Store the output axis powers associated with each coefficient of the inverse
   transformation. */
      iv = 1;
      for( i = 0; i < nin; i++ ){
         for( j = 0; j < this->ncoeff_i[ i ]; j++ ){
            for( k = 0; k < nout; k++, iv++ ){
               if( (this->power_i)[ i ][ j ][ k ] > 0 ) {
                  (void) sprintf( buff, "PI%d", iv );
                  (void) sprintf( comm, "Power of o/p %d for coeff %d of inv poly %d", k + 1, j + 1, i + 1 );
                  astWriteDouble( channel, buff, 1, 1, (this->power_i)[ i ][ j ][ k ], comm );
               }
            }
         }
      }
   }

/* Use an iterative inverse? */
   set = TestIterInverse( this, status );
   ival = set ? GetIterInverse( this, status ) : astGetIterInverse( this );
   astWriteInt( channel, "IterInv", set, 0, ival, ival ? "Use an iterative inverse" : "Do not use an iterative inverse" );

/* Max number of iterations for iterative inverse. */
   set = TestNiterInverse( this, status );
   ival = set ? GetNiterInverse( this, status ) : astGetNiterInverse( this );
   astWriteInt( channel, "NiterInv", set, 0, ival, "Max number of iterations for iterative inverse" );

/* Target relative error for iterative inverse. */
   set = TestTolInverse( this, status );
   dval = set ? GetTolInverse( this, status ) : astGetTolInverse( this );
   astWriteDouble( channel, "TolInv", set, 0, dval, "Target relative error for iterative inverse" );

/* Undefine macros local to this function. */
#undef KEY_LEN
}

/* Standard class functions. */
/* ========================= */
/* Implement the astIsAPolyMap and astCheckPolyMap functions using the macros
   defined for this purpose in the "object.h" header file. */
astMAKE_ISA(PolyMap,Mapping)
astMAKE_CHECK(PolyMap)

AstPolyMap *astPolyMap_( int nin, int nout, int ncoeff_f, const double coeff_f[],
                         int ncoeff_i, const double coeff_i[], const char *options, int *status, ...){
/*
*++
*  Name:
c     astPolyMap
f     AST_POLYMAP

*  Purpose:
*     Create a PolyMap.

*  Type:
*     Public function.

*  Synopsis:
c     #include "polymap.h"
c     AstPolyMap *astPolyMap( int nin, int nout, int ncoeff_f, const double coeff_f[],
c                             int ncoeff_i, const double coeff_i[],
c                             const char *options, ... )
f     RESULT = AST_POLYMAP( NIN, NOUT, NCOEFF_F, COEFF_F, NCOEFF_I, COEFF_I,
f                           OPTIONS, STATUS )

*  Class Membership:
*     PolyMap constructor.

*  Description:
*     This function creates a new PolyMap and optionally initialises
*     its attributes.
*
*     A PolyMap is a form of Mapping which performs a general polynomial
*     transformation.  Each output coordinate is a polynomial function of
*     all the input coordinates. The coefficients are specified separately
*     for each output coordinate. The forward and inverse transformations
*     are defined independantly by separate sets of coefficients. If no
*     inverse transformation is supplied, the default behaviour is to use
*     an iterative method to evaluate the inverse based only on the forward
*     transformation (see attribute IterInverse).

*  Parameters:
c     nin
f     NIN = INTEGER (Given)
*        The number of input coordinates.
c     nout
f     NOUT = INTEGER (Given)
*        The number of output coordinates.
c     ncoeff_f
f     NCOEFF_F = INTEGER (Given)
*        The number of non-zero coefficients necessary to define the
*        forward transformation of the PolyMap. If zero is supplied, the
*        forward transformation will be undefined.
c     coeff_f
f     COEFF_F( * ) = DOUBLE PRECISION (Given)
*        An array containing
c        "ncoeff_f*( 2 + nin )" elements. Each group of "2 + nin"
f        "NCOEFF_F*( 2 + NIN )" elements. Each group of "2 + NIN"
*        adjacent elements describe a single coefficient of the forward
*        transformation. Within each such group, the first element is the
*        coefficient value; the next element is the integer index of the
*        PolyMap output which uses the coefficient within its defining
*        polynomial (the first output has index 1); the remaining elements
*        of the group give the integer powers to use with each input
*        coordinate value (powers must not be negative, and floating
*        point values are rounded to the nearest integer).
c        If "ncoeff_f" is zero, a NULL pointer may be supplied for "coeff_f".
*
*        For instance, if the PolyMap has 3 inputs and 2 outputs, each group
*        consisting of 5 elements, A groups such as "(1.2, 2.0, 1.0, 3.0, 0.0)"
*        describes a coefficient with value 1.2 which is used within the
*        definition of output 2. The output value is incremented by the
*        product of the coefficient value, the value of input coordinate
*        1 raised to the power 1, and the value of input coordinate 2 raised
*        to the power 3. Input coordinate 3 is not used since its power is
*        specified as zero. As another example, the group "(-1.0, 1.0,
*        0.0, 0.0, 0.0 )" describes adds a constant value -1.0 onto
*        output 1 (it is a constant value since the power for every input
*        axis is given as zero).
*
c        Each final output coordinate value is the sum of the "ncoeff_f" terms
c        described by the "ncoeff_f" groups within the supplied array.
f        Each final output coordinate value is the sum of the "NCOEFF_F" terms
f        described by the "NCOEFF_F" groups within the supplied array.
c     ncoeff_i
f     NCOEFF_I = INTEGER (Given)
*        The number of non-zero coefficients necessary to define the
*        inverse transformation of the PolyMap. If zero is supplied, the
*        default behaviour is to use an iterative method to evaluate the
*        inverse based only on the forward transformation (see attribute
*        IterInverse).
c     coeff_i
f     COEFF_I( * ) = DOUBLE PRECISION (Given)
*        An array containing
c        "ncoeff_i*( 2 + nout )" elements. Each group of "2 + nout"
f        "NCOEFF_I*( 2 + NOUT )" elements. Each group of "2 + NOUT"
*        adjacent elements describe a single coefficient of the inverse
c        transformation, using the same schame as "coeff_f",
f        transformation, using the same schame as "COEFF_F",
*        except that "inputs" and "outputs" are transposed.
c        If "ncoeff_i" is zero, a NULL pointer may be supplied for "coeff_i".
c     options
f     OPTIONS = CHARACTER * ( * ) (Given)
c        Pointer to a null-terminated string containing an optional
c        comma-separated list of attribute assignments to be used for
c        initialising the new PolyMap. The syntax used is identical to
c        that for the astSet function and may include "printf" format
c        specifiers identified by "%" symbols in the normal way.
f        A character string containing an optional comma-separated
f        list of attribute assignments to be used for initialising the
f        new PolyMap. The syntax used is identical to that for the
f        AST_SET routine.
c     ...
c        If the "options" string contains "%" format specifiers, then
c        an optional list of additional arguments may follow it in
c        order to supply values to be substituted for these
c        specifiers. The rules for supplying these are identical to
c        those for the astSet function (and for the C "printf"
c        function).
f     STATUS = INTEGER (Given and Returned)
f        The global status.

*  Returned Value:
c     astPolyMap()
f     AST_POLYMAP = INTEGER
*        A pointer to the new PolyMap.

*  Notes:
*     - A null Object pointer (AST__NULL) will be returned if this
c     function is invoked with the AST error status set, or if it
f     function is invoked with STATUS set to an error value, or if it
*     should fail for any reason.
*--
*/

/* Local Variables: */
   astDECLARE_GLOBALS          /* Pointer to thread-specific global data */
   AstPolyMap *new;            /* Pointer to new PolyMap */
   va_list args;               /* Variable argument list */

/* Check the global status. */
   if ( !astOK ) return NULL;

/* Get a pointer to the thread specific global data structure. */
   astGET_GLOBALS(NULL);

/* Initialise the PolyMap, allocating memory and initialising the
   virtual function table as well if necessary. */
   new = astInitPolyMap( NULL, sizeof( AstPolyMap ), !class_init,
                         &class_vtab, "PolyMap", nin, nout,
                         ncoeff_f, coeff_f, ncoeff_i, coeff_i );

/* If successful, note that the virtual function table has been
   initialised. */
   if ( astOK ) {
      class_init = 1;

/* Obtain the variable argument list and pass it along with the options string
   to the astVSet method to initialise the new PolyMap's attributes. */
      va_start( args, status );
      astVSet( new, options, NULL, args );
      va_end( args );

/* If an error occurred, clean up by deleting the new object. */
      if ( !astOK ) new = astDelete( new );
   }

/* Return a pointer to the new PolyMap. */
   return new;
}

AstPolyMap *astPolyMapId_( int nin, int nout, int ncoeff_f, const double coeff_f[],
                           int ncoeff_i, const double coeff_i[], const char *options, ... ){
/*
*  Name:
*     astPolyMapId_

*  Purpose:
*     Create a PolyMap.

*  Type:
*     Private function.

*  Synopsis:
*     #include "polymap.h"
*     AstPolyMap *astPolyMap( int nin, int nout, int ncoeff_f, const double coeff_f[],
*                             int ncoeff_i, const double coeff_i[],
*                             const char *options, ... )

*  Class Membership:
*     PolyMap constructor.

*  Description:
*     This function implements the external (public) interface to the
*     astPolyMap constructor function. It returns an ID value (instead
*     of a true C pointer) to external users, and must be provided
*     because astPolyMap_ has a variable argument list which cannot be
*     encapsulated in a macro (where this conversion would otherwise
*     occur).
*
*     The variable argument list also prevents this function from
*     invoking astPolyMap_ directly, so it must be a re-implementation
*     of it in all respects, except for the final conversion of the
*     result to an ID value.

*  Parameters:
*     As for astPolyMap_.

*  Returned Value:
*     The ID value associated with the new PolyMap.
*/

/* Local Variables: */
   astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
   AstPolyMap *new;              /* Pointer to new PolyMap */
   va_list args;                 /* Variable argument list */
   int *status;                  /* Pointer to inherited status value */

/* Get a pointer to the inherited status value. */
   status = astGetStatusPtr;

/* Get a pointer to the thread specific global data structure. */
   astGET_GLOBALS(NULL);

/* Check the global status. */
   if ( !astOK ) return NULL;

/* Initialise the PolyMap, allocating memory and initialising the
   virtual function table as well if necessary. */
   new = astInitPolyMap( NULL, sizeof( AstPolyMap ), !class_init,
                         &class_vtab, "PolyMap", nin, nout,
                         ncoeff_f, coeff_f, ncoeff_i, coeff_i );

/* If successful, note that the virtual function table has been
   initialised. */
   if ( astOK ) {
      class_init = 1;

/* Obtain the variable argument list and pass it along with the options string
   to the astVSet method to initialise the new PolyMap's attributes. */
      va_start( args, options );
      astVSet( new, options, NULL, args );
      va_end( args );

/* If an error occurred, clean up by deleting the new object. */
      if ( !astOK ) new = astDelete( new );
   }

/* Return an ID value for the new PolyMap. */
   return astMakeId( new );
}

AstPolyMap *astInitPolyMap_( void *mem, size_t size, int init,
                             AstPolyMapVtab *vtab, const char *name,
                             int nin, int nout, int ncoeff_f, const double coeff_f[],
                             int ncoeff_i, const double coeff_i[], int *status ){
/*
*+
*  Name:
*     astInitPolyMap

*  Purpose:
*     Initialise a PolyMap.

*  Type:
*     Protected function.

*  Synopsis:
*     #include "polymap.h"
*     AstPolyMap *astInitPolyMap( void *mem, size_t size, int init,
*                                 AstPolyMapVtab *vtab, const char *name,
*                                 int nin, int nout, int ncoeff_f,
*                                 const double coeff_f[], int ncoeff_i,
*                                 const double coeff_i[] )

*  Class Membership:
*     PolyMap initialiser.

*  Description:
*     This function is provided for use by class implementations to initialise
*     a new PolyMap object. It allocates memory (if necessary) to accommodate
*     the PolyMap plus any additional data associated with the derived class.
*     It then initialises a PolyMap structure at the start of this memory. If
*     the "init" flag is set, it also initialises the contents of a virtual
*     function table for a PolyMap at the start of the memory passed via the
*     "vtab" parameter.

*  Parameters:
*     mem
*        A pointer to the memory in which the PolyMap is to be initialised.
*        This must be of sufficient size to accommodate the PolyMap data
*        (sizeof(PolyMap)) plus any data used by the derived class. If a value
*        of NULL is given, this function will allocate the memory itself using
*        the "size" parameter to determine its size.
*     size
*        The amount of memory used by the PolyMap (plus derived class data).
*        This will be used to allocate memory if a value of NULL is given for
*        the "mem" parameter. This value is also stored in the PolyMap
*        structure, so a valid value must be supplied even if not required for
*        allocating memory.
*     init
*        A logical flag indicating if the PolyMap's virtual function table is
*        to be initialised. If this value is non-zero, the virtual function
*        table will be initialised by this function.
*     vtab
*        Pointer to the start of the virtual function table to be associated
*        with the new PolyMap.
*     name
*        Pointer to a constant null-terminated character string which contains
*        the name of the class to which the new object belongs (it is this
*        pointer value that will subsequently be returned by the astGetClass
*        method).
*     nin
*        The number of input coordinate values per point. This is the
*        same as the number of columns in the matrix.
*     nout
*        The number of output coordinate values per point. This is the
*        same as the number of rows in the matrix.
*     ncoeff_f
*        The number of non-zero coefficients necessary to define the
*        forward transformation of the PolyMap. If zero is supplied, the
*        forward transformation will be undefined.
*     coeff_f
*        An array containing "ncoeff_f*( 2 + nin )" elements. Each group
*	 of "2 + nin" adjacent elements describe a single coefficient of
*	 the forward transformation. Within each such group, the first
*	 element is the coefficient value; the next element is the
*	 integer index of the PolyMap output which uses the coefficient
*	 within its defining polynomial (the first output has index 1);
*	 the remaining elements of the group give the integer powers to
*	 use with each input coordinate value (powers must not be
*	 negative)
*
*        For instance, if the PolyMap has 3 inputs and 2 outputs, each group
*        consisting of 5 elements, A groups such as "(1.2, 2.0, 1.0, 3.0, 0.0)"
*        describes a coefficient with value 1.2 which is used within the
*        definition of output 2. The output value is incremented by the
*        product of the coefficient value, the value of input coordinate
*        1 raised to the power 1, and the value of input coordinate 2 raised
*        to the power 3. Input coordinate 3 is not used since its power is
*        specified as zero. As another example, the group "(-1.0, 1.0,
*        0.0, 0.0, 0.0 )" describes adds a constant value -1.0 onto
*        output 1 (it is a constant value since the power for every input
*        axis is given as zero).
*
*        Each final output coordinate value is the sum of the "ncoeff_f" terms
*        described by the "ncoeff_f" groups within the supplied array.
*     ncoeff_i
*        The number of non-zero coefficients necessary to define the
*        inverse transformation of the PolyMap. If zero is supplied, the
*        inverse transformation will be undefined.
*     coeff_i
*        An array containing
*        "ncoeff_i*( 2 + nout )" elements. Each group of "2 + nout"
*        adjacent elements describe a single coefficient of the inverse
*        transformation, using the same schame as "coeff_f", except that
*        "inputs" and "outputs" are transposed.

*  Returned Value:
*     A pointer to the new PolyMap.

*  Notes:
*     -  A null pointer will be returned if this function is invoked with the
*     global error status set, or if it should fail for any reason.
*-
*/

/* Local Variables: */
   AstPolyMap *new;              /* Pointer to new PolyMap */

/* Check the global status. */
   if ( !astOK ) return NULL;

/* If necessary, initialise the virtual function table. */
   if ( init ) astInitPolyMapVtab( vtab, name );

/* Initialise a Mapping structure (the parent class) as the first component
   within the PolyMap structure, allocating memory if necessary. Specify that
   the Mapping should be defined in both the forward and inverse directions. */
   new = (AstPolyMap *) astInitMapping( mem, size, 0,
                                        (AstMappingVtab *) vtab, name,
                                        nin, nout, 1, 1 );
   if ( astOK ) {

/* Initialise the PolyMap data. */
/* ---------------------------- */

/* First initialise the pointers in case of errors. */
      new->ncoeff_f = NULL;
      new->power_f = NULL;
      new->coeff_f = NULL;
      new->mxpow_f = NULL;

      new->ncoeff_i = NULL;
      new->power_i = NULL;
      new->coeff_i = NULL;
      new->mxpow_i = NULL;

/* Store the forward transformation. */
      StoreArrays( new, 1, ncoeff_f, coeff_f, status );

/* Store the inverse transformation. */
      StoreArrays( new, 0, ncoeff_i, coeff_i, status );

/* Other class attributes. */
      new->iterinverse = -INT_MAX;
      new->niterinverse = -INT_MAX;
      new->tolinverse = AST__BAD;
      new->jacobian = NULL;
      new->lintrunc = NULL;

/* If an error occurred, clean up by deleting the new PolyMap. */
      if ( !astOK ) new = astDelete( new );
   }

/* Return a pointer to the new PolyMap. */
   return new;
}

AstPolyMap *astLoadPolyMap_( void *mem, size_t size,
                             AstPolyMapVtab *vtab, const char *name,
                             AstChannel *channel, int *status ) {
/*
*+
*  Name:
*     astLoadPolyMap

*  Purpose:
*     Load a PolyMap.

*  Type:
*     Protected function.

*  Synopsis:
*     #include "polymap.h"
*     AstPolyMap *astLoadPolyMap( void *mem, size_t size,
*                                 AstPolyMapVtab *vtab, const char *name,
*                                 AstChannel *channel )

*  Class Membership:
*     PolyMap loader.

*  Description:
*     This function is provided to load a new PolyMap using data read
*     from a Channel. It first loads the data used by the parent class
*     (which allocates memory if necessary) and then initialises a
*     PolyMap structure in this memory, using data read from the input
*     Channel.
*
*     If the "init" flag is set, it also initialises the contents of a
*     virtual function table for a PolyMap at the start of the memory
*     passed via the "vtab" parameter.


*  Parameters:
*     mem
*        A pointer to the memory into which the PolyMap is to be
*        loaded.  This must be of sufficient size to accommodate the
*        PolyMap data (sizeof(PolyMap)) plus any data used by derived
*        classes. If a value of NULL is given, this function will
*        allocate the memory itself using the "size" parameter to
*        determine its size.
*     size
*        The amount of memory used by the PolyMap (plus derived class
*        data).  This will be used to allocate memory if a value of
*        NULL is given for the "mem" parameter. This value is also
*        stored in the PolyMap structure, so a valid value must be
*        supplied even if not required for allocating memory.
*
*        If the "vtab" parameter is NULL, the "size" value is ignored
*        and sizeof(AstPolyMap) is used instead.
*     vtab
*        Pointer to the start of the virtual function table to be
*        associated with the new PolyMap. If this is NULL, a pointer
*        to the (static) virtual function table for the PolyMap class
*        is used instead.
*     name
*        Pointer to a constant null-terminated character string which
*        contains the name of the class to which the new object
*        belongs (it is this pointer value that will subsequently be
*        returned by the astGetClass method).
*
*        If the "vtab" parameter is NULL, the "name" value is ignored
*        and a pointer to the string "PolyMap" is used instead.

*  Returned Value:
*     A pointer to the new PolyMap.

*  Notes:
*     - A null pointer will be returned if this function is invoked
*     with the global error status set, or if it should fail for any
*     reason.
*-
*/

#define KEY_LEN 50               /* Maximum length of a keyword */

   astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
/* Local Variables: */
   AstPolyMap *new;              /* Pointer to the new PolyMap */
   char buff[ KEY_LEN + 1 ];     /* Buffer for keyword string */
   int i;                        /* Loop index */
   int iv;                       /* Vectorised keyword index */
   int j;                        /* Loop index */
   int k;                        /* Loop index */
   int nin;                      /* No. of input coords */
   int nout;                     /* No. of output coords */
   int undef;                    /* Is the transformation undefined? */

/* Get a pointer to the thread specific global data structure. */
   astGET_GLOBALS(channel);

/* Initialise. */
   new = NULL;

/* Check the global error status. */
   if ( !astOK ) return new;

/* If a NULL virtual function table has been supplied, then this is
   the first loader to be invoked for this PolyMap. In this case the
   PolyMap belongs to this class, so supply appropriate values to be
   passed to the parent class loader (and its parent, etc.). */
   if ( !vtab ) {
      size = sizeof( AstPolyMap );
      vtab = &class_vtab;
      name = "PolyMap";

/* If required, initialise the virtual function table for this class. */
      if ( !class_init ) {
         astInitPolyMapVtab( vtab, name );
         class_init = 1;
      }
   }

/* Invoke the parent class loader to load data for all the ancestral
   classes of the current one, returning a pointer to the resulting
   partly-built PolyMap. */
   new = astLoadMapping( mem, size, (AstMappingVtab *) vtab, name,
                         channel );

   if ( astOK ) {

/* Get the number of inputs and outputs for the uninverted Mapping. */
   nin = ( (AstMapping *) new )->nin;
   nout = ( (AstMapping *) new )->nout;

/* Read input data. */
/* ================ */
/* Request the input Channel to read all the input data appropriate to
   this class into the internal "values list". */
      astReadClassData( channel, "PolyMap" );

/* Allocate memory to hold the forward arrays. */
      new->ncoeff_f = astMalloc( sizeof( int )*(size_t) nout );
      new->mxpow_f = astMalloc( sizeof( int )*(size_t) nin );
      new->power_f = astMalloc( sizeof( int ** )*(size_t) nout );
      new->coeff_f = astMalloc( sizeof( double * )*(size_t) nout );
      if( astOK ) {

/* Assume the forward transformation is defined. */
         undef = 0;

/* Get the maximum power of each input axis value used by the forward
   transformation. Set a flag "undef" if no values relating to the
   forward transformation are found (this indicates that the forward
   transformation is not defined). */
         for( i = 0; i < nin && !undef; i++ ){
            (void) sprintf( buff, "mpf%d", i + 1 );
            (new->mxpow_f)[ i ] = astReadInt( channel, buff, INT_MAX );
            if( (new->mxpow_f)[ i ] == INT_MAX ) undef = 1;
         }

/* Get the number of coefficients associated with each output of the forward
   transformation. */
         for( i = 0; i < nout && !undef; i++ ){
            (void) sprintf( buff, "ncf%d", i + 1 );
            (new->ncoeff_f)[ i ] = astReadInt( channel, buff, INT_MAX );
            if( (new->ncoeff_f)[ i ] == INT_MAX ) undef = 1;
         }

/* Get the coefficient values used by the forward transformation. This
   uses new style vectorised key names if available. Otherwise it uses
   old style indexed names (which were superceded by vectorised names
   because they are shorter and so work better with FitsChans). */
         iv = 0;
         for( i = 0; i < nout && !undef; i++ ){

            (new->coeff_f)[ i ] = astMalloc( sizeof( double )*
                                           (size_t) new->ncoeff_f[ i ] );
            if( astOK ) {
               for( j = 0; j < new->ncoeff_f[ i ]; j++ ){
                  (void) sprintf( buff, "cf%d", ++iv );
                  (new->coeff_f)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD );
                  if( (new->coeff_f)[ i ][ j ] == AST__BAD ) {
                     (void) sprintf( buff, "cf%d_%d", i + 1, j + 1 );
                     (new->coeff_f)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD );
                  }
               }
            }
         }

/* Get the input axis powers associated with each coefficient of the forward
   transformation. */
         iv = 0;
         for( i = 0; i < nout && !undef; i++ ){
            (new->power_f)[ i ] = astMalloc( sizeof( int * )*
                                             (size_t) new->ncoeff_f[ i ] );
            if( astOK ) {
               for( j = 0; j < new->ncoeff_f[ i ]; j++ ){
                  (new->power_f)[ i ][ j ] = astMalloc( sizeof( int )* (size_t) nin );
                  if( astOK ) {
                     for( k = 0; k < nin; k++ ){
                        (void) sprintf( buff, "pf%d", ++iv );
                        (new->power_f)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 );
                        if( (new->power_f)[ i ][ j ][ k ] == 0 ) {
                           (void) sprintf( buff, "pf%d_%d_%d", i + 1, j + 1, k + 1 );
                           (new->power_f)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 );
                        }
                     }
                  }
               }
            }
         }

/* Free the arrays if the forward transformation is undefined. */
         if( undef ) {
            new->ncoeff_f = astFree( new->ncoeff_f );
            new->mxpow_f = astFree( new->mxpow_f );
            new->power_f = astFree( new->power_f );
            new->coeff_f = astFree( new->coeff_f );
         }
      }

/* Allocate memory to hold the inverse arrays. */
      new->ncoeff_i = astMalloc( sizeof( int )*(size_t) nin );
      new->mxpow_i = astMalloc( sizeof( int )*(size_t) nout );
      new->power_i = astMalloc( sizeof( int ** )*(size_t) nin );
      new->coeff_i = astMalloc( sizeof( double * )*(size_t) nin );
      if( astOK ) {

/* Assume the inverse transformation is defined. */
         undef = 0;

/* Get the maximum power of each output axis value used by the inverse
   transformation. Set a flag "undef" if no values relating to the
   inverse transformation are found (this indicates that the inverse
   transformation is not defined). */
         for( i = 0; i < nout && !undef; i++ ){
            (void) sprintf( buff, "mpi%d", i + 1 );
            (new->mxpow_i)[ i ] = astReadInt( channel, buff, INT_MAX );
            if( (new->mxpow_i)[ i ] == INT_MAX ) undef = 1;
         }

/* Get the number of coefficients associated with each input of the inverse
   transformation. */
         for( i = 0; i < nin && !undef; i++ ){
            (void) sprintf( buff, "nci%d", i + 1 );
            (new->ncoeff_i)[ i ] = astReadInt( channel, buff, INT_MAX );
            if( (new->ncoeff_i)[ i ] == INT_MAX ) undef = 1;
         }

/* Get the coefficient values used by the inverse transformation. */
         iv = 0;
         for( i = 0; i < nin && !undef; i++ ){

            (new->coeff_i)[ i ] = astMalloc( sizeof( double )*
                                           (size_t) new->ncoeff_i[ i ] );
            if( astOK ) {
               for( j = 0; j < new->ncoeff_i[ i ]; j++ ){
                  (void) sprintf( buff, "ci%d", ++iv );
                  (new->coeff_i)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD );
                  if( (new->coeff_i)[ i ][ j ] == AST__BAD ) {
                     (void) sprintf( buff, "ci%d_%d", i + 1, j + 1 );
                     (new->coeff_i)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD );
                  }
               }
            }
         }

/* Get the output axis powers associated with each coefficient of the inverse
   transformation. */
         iv = 0;
         for( i = 0; i < nin && !undef; i++ ){
            (new->power_i)[ i ] = astMalloc( sizeof( int * )*
                                             (size_t) new->ncoeff_i[ i ] );
            if( astOK ) {
               for( j = 0; j < new->ncoeff_i[ i ]; j++ ){
                  (new->power_i)[ i ][ j ] = astMalloc( sizeof( int )* (size_t) nout );
                  if( astOK ) {
                     for( k = 0; k < nout; k++ ){
                        (void) sprintf( buff, "pi%d", ++iv );
                        (new->power_i)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 );
                        if( (new->power_i)[ i ][ j ][ k ] == 0 ) {
                           (void) sprintf( buff, "pi%d_%d_%d", i + 1, j + 1, k + 1 );
                           (new->power_i)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 );
                        }
                     }
                  }
               }
            }
         }

/* Free the arrays if the inverse transformation is undefined. */
         if( undef ) {
            new->ncoeff_i = astFree( new->ncoeff_i );
            new->mxpow_i = astFree( new->mxpow_i );
            new->power_i = astFree( new->power_i );
            new->coeff_i = astFree( new->coeff_i );
         }
      }

/* Whether to use an iterative inverse transformation. */
      new->iterinverse = astReadInt( channel, "iterinv", -INT_MAX );
      if ( TestIterInverse( new, status ) ) SetIterInverse( new, new->iterinverse, status );

/* Max number of iterations for iterative inverse transformation. */
      new->niterinverse = astReadInt( channel, "niterinv", -INT_MAX );
      if ( TestNiterInverse( new, status ) ) SetNiterInverse( new, new->niterinverse, status );

/* Target relative error for iterative inverse transformation. */
      new->tolinverse = astReadDouble( channel, "tolinv", AST__BAD );
      if ( TestTolInverse( new, status ) ) SetTolInverse( new, new->tolinverse, status );

/* The Jacobian of the PolyMap's forward transformation has not yet been
   found. */
      new->jacobian = NULL;

/* The linear truncation of the PolyMap has not yet been found. */
      new->lintrunc = NULL;

/* If an error occurred, clean up by deleting the new PolyMap. */
      if ( !astOK ) new = astDelete( new );
   }

/* Return the new PolyMap pointer. */
   return new;

/* Undefine macros local to this function. */
#undef KEY_LEN
}

/* Virtual function interfaces. */
/* ============================ */
/* These provide the external interface to the virtual functions defined by
   this class. Each simply checks the global error status and then locates and
   executes the appropriate member function, using the function pointer stored
   in the object's virtual function table (this pointer is located using the
   astMEMBER macro defined in "object.h").

   Note that the member function may not be the one defined here, as it may
   have been over-ridden by a derived class. However, it should still have the
   same interface. */

void astPolyPowers_( AstPolyMap *this, double **work, int ncoord,
                     const int *mxpow, double **ptr, int point, int fwd,
                     int *status ){
   if ( !astOK ) return;
   (**astMEMBER(this,PolyMap,PolyPowers))( this, work, ncoord, mxpow, ptr,
                                           point, fwd, status );
}

AstPolyMap *astPolyTran_( AstPolyMap *this, int forward, double acc,
                          double maxacc, int maxorder, const double *lbnd,
                          const double *ubnd, int *status ){
   if ( !astOK ) return NULL;
   return (**astMEMBER(this,PolyMap,PolyTran))( this, forward, acc,
                                                maxacc, maxorder, lbnd,
                                                ubnd, status );
}

void astPolyCoeffs_( AstPolyMap *this, int forward, int nel, double *array,
                     int *ncoeff, int *status ){
   if ( !astOK ) return;
   (**astMEMBER(this,PolyMap,PolyCoeffs))( this, forward, nel,
                                           array, ncoeff, status );
}

void astFitPoly1DInit_( AstPolyMap *this, int forward, double **table,
                        AstMinPackData *data, double *scales,
                        int *status ){
   if ( !astOK ) return;
   (**astMEMBER(this,PolyMap,FitPoly1DInit))( this, forward, table, data, scales,
                                              status );
}


void astFitPoly2DInit_( AstPolyMap *this, int forward, double **table,
                        AstMinPackData *data, double *scales,
                        int *status ){
   if ( !astOK ) return;
   (**astMEMBER(this,PolyMap,FitPoly2DInit))( this, forward, table, data, scales,
                                              status );
}