path: root/ast/lutmap.c

/*
*class++
*  Name:
*     LutMap

*  Purpose:
*     Transform 1-dimensional coordinates using a lookup table.

*  Constructor Function:
c     astLutMap
f     AST_LUTMAP

*  Description:
*     A LutMap is a specialised form of Mapping which transforms
*     1-dimensional coordinates by using linear interpolation in a
*     lookup table.
*
*     Each input coordinate value is first scaled to give the index of
*     an entry in the table by subtracting a starting value (the input
*     coordinate corresponding to the first table entry) and dividing
*     by an increment (the difference in input coordinate value
*     between adjacent table entries).
*
*     The resulting index will usually contain a fractional part, so
*     the output coordinate value is then generated by interpolating
*     linearly between the appropriate entries in the table. If the
*     index lies outside the range of the table, linear extrapolation
*     is used based on the two nearest entries (i.e. the two entries
*     at the start or end of the table, as appropriate). If either of the
*     entries used for the interplation has a value of AST__BAD, then the
*     interpolated value is returned as AST__BAD.
*
*     If the lookup table entries increase or decrease monotonically
*     (ignoring any flat sections), then the inverse transformation may
*     also be performed.

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

*  Attributes:
*     In addition to those attributes common to all Mappings, every
*     LutMap also has the following attributes:
*
*     - LutEpsilon: The relative error of the values in the table.
*     - LutInterp: The interpolation method to use between table entries.

*  Functions:
c     The LutMap class does not define any new functions beyond those
f     The LutMap class does not define any new routines beyond those
*     which are applicable to all Mappings.

*  Copyright:
*     Copyright (C) 1997-2006 Council for the Central Laboratory of the
*     Research Councils
*     Copyright (C) 2007-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:
*     RFWS: R.F. Warren-Smith (Starlink)
*     DSB: David S. Berry (JAC, UCLan)

*  History:
*     8-JUL-1997 (RFWS):
*        Original version.
*     10-JUL-1997 (RFWS):
*        Added the MapMerge function.
*     8-JAN-2003 (DSB):
*        Changed private InitVtab method to protected astInitLutMapVtab
*        method.
*     12-JAN-2004 (DSB):
*        Check for AST__BAD values in the supplied lut array.
*     17-MAR-2006 (DSB):
*        - MapMerge changed so that a LutMap will cancel with its own
*        inverse.
*        - Added attribute LutInterp
*     10-MAY-2006 (DSB):
*        Override astEqual.
*     4-OCT-2006 (DSB):
*        - Correct "mintick" to "lutinterp" in SetAttrib.
*        - Do not include bad values in the dumped LUT array.
*     8-NOV-2007 (DSB):
*        - Take account of the requested invert flag when comparing two
*        neighbouring LutMaps for equality in MapMerge.
*     19-NOV-2010 (DSB):
*        Added (protected) astGetLutMapInfo function.
*     24-JAN-2011 (DSB):
*        Implement an inverse transformation even if the coordinate
*        array contains sections of equal or bad values. The inverse
*        transformation will generate bad values if used within such
*        regions of the coordinate array.
*     6-JUL-2011 (DSB):
*        Avoid indexing the lut array beyond the end when doing an
*        inverse transform.
*     2-OCT-2012 (DSB):
*        Check for Infs as well as NaNs.
*     21-MAY-2015 (DSB):
*        Aded LutEpsilon
*     23-SEP-2015 (DSB):
*        The GetMonotonic function had a bug that caused all LutMaps
*        to be considered monotonic, and thus have an inverse
*        transformation.
*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 LutMap

#define LINEAR 0
#define NEAR 1

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

#include "globals.h"             /* Thread-safe global data access */
#include "error.h"               /* Error reporting facilities */
#include "memory.h"              /* Memory management facilities */
#include "object.h"              /* Base Object class */
#include "pointset.h"            /* Sets of points/coordinates */
#include "mapping.h"             /* Coordinate mappings (parent class) */
#include "winmap.h"              /* Linear mappings between windows */
#include "channel.h"             /* I/O channels */
#include "unitmap.h"             /* Unit mappings */
#include "lutmap.h"              /* Interface definition for this class */
#include "globals.h"             /* Thread-safe global data access */

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

/* C header files. */
/* --------------- */
#include <float.h>
#include <math.h>
#include <limits.h>
#include <stdarg.h>
#include <stddef.h>
#include <stdio.h>
#include <string.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 * );

/* Define macros for accessing each item of thread specific global data. */
#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(LutMap)

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



/* If thread safety is not needed, declare and initialise globals at static
   variables. */
#else

static char getattrib_buff[ 101 ];


/* Define the class virtual function table and its initialisation flag
   as static variables. */
static AstLutMapVtab 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. */
AstLutMap *astLutMapId_( int, const double [], double, double, const char *, ... );

/* Prototypes for Private Member Functions. */
/* ======================================== */
static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
static int GetLinear( AstMapping *, int * );
static int GetMonotonic( int, const double *, int *, double **, int **, int **, int * );
static int MapMerge( AstMapping *, int, int, int *, AstMapping ***, int **, int * );
static void Copy( const AstObject *, AstObject *, int * );
static void Delete( AstObject *, int * );
static void Dump( AstObject *, AstChannel *, int * );
static int Equal( AstObject *, AstObject *, int * );
static double *GetLutMapInfo( AstLutMap *, double *, double *, int *, int * );

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 GetLutInterp( AstLutMap *, int * );
static int TestLutInterp( AstLutMap *, int * );
static void ClearLutInterp( AstLutMap *, int * );
static void SetLutInterp( AstLutMap *, int, int * );

static double GetLutEpsilon( AstLutMap *, int * );
static int TestLutEpsilon( AstLutMap *, int * );
static void ClearLutEpsilon( AstLutMap *, int * );
static void SetLutEpsilon( AstLutMap *, double, int * );

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

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

*  Type:
*     Private function.

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

*  Class Membership:
*     LutMap 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
*     LutMap, so that the default value will subsequently be used.

*  Parameters:
*     this
*        Pointer to the LutMap.
*     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: */
   AstLutMap *this;             /* Pointer to the LutMap structure */

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

/* Obtain a pointer to the LutMap structure. */
   this = (AstLutMap *) this_object;

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

/* LutInterp. */
/* ---------- */
   if ( !strcmp( attrib, "lutinterp" ) ) {
      astClearLutInterp( this );

/* LutEpsilon. */
/* ------------- */
   } else if ( !strcmp( attrib, "lutepsilon" ) ) {
      astClearLutEpsilon( 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 LutMaps are equivalent.

*  Type:
*     Private function.

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

*  Class Membership:
*     LutMap 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 LutMaps are equivalent.

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

*  Returned Value:
*     One if the LutMaps 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: */
   AstLutMap *that;
   AstLutMap *this;
   int i;
   int nin;
   int nout;
   int result;

/* Initialise. */
   result = 0;

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

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

/* Check the second object is a LutMap. We know the first is a
   LutMap since we have arrived at this implementation of the virtual
   function. */
   if( astIsALutMap( 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 LutMaps differ, it may still be possible
   for them to be equivalent. First compare the LutMaps if their Invert
   flags are the same. In this case all the attributes of the two LutMaps
   must be identical. */
         if( astGetInvert( this ) == astGetInvert( that ) ) {

            if( astEQUAL( this->start, that->start ) &&
                astEQUAL( this->inc, that->inc ) &&
                this->nlut == that->nlut &&
                this->lutinterp == that->lutinterp ){

               result = 1;
               for( i = 0; i < this->nlut; i++ ) {
                  if( !astEQUAL( (this->lut)[ i ], (that->lut)[ i ] ) ) {
                     result = 0;
                     break;
                  }
               }
            }

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

/* In the specific case of a LutMap, 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 const char *GetAttrib( AstObject *this_object, const char *attrib, int *status ) {
/*
*  Name:
*     GetAttrib

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

*  Type:
*     Private function.

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

*  Class Membership:
*     LutMap 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 LutMap, formatted as a character string.

*  Parameters:
*     this
*        Pointer to the LutMap.
*     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 LutMap, 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 LutMap. 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: */
   AstLutMap *this;             /* Pointer to the LutMap structure */
   astDECLARE_GLOBALS           /* Pointer to thread-specific global data */
   const char *result;          /* Pointer value to return */
   double luteps;               /* LutEpsilon attribute value */
   int lutinterp;               /* LutInterp 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 LutMap structure. */
   this = (AstLutMap *) 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. */

/* LutInterp. */
/* ---------- */
   if ( !strcmp( attrib, "lutinterp" ) ) {
      lutinterp = astGetLutInterp( this );
      if ( astOK ) {
         (void) sprintf( getattrib_buff, "%d", lutinterp );
         result = getattrib_buff;
      }

/* LutEpsilon. */
/* ------------- */
   } else if ( !strcmp( attrib, "lutepsilon" ) ) {
      luteps = astGetLutEpsilon( this );
      if ( astOK ) {
         (void) sprintf( getattrib_buff, "%.*g", AST__DBL_DIG, luteps );
         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 int GetLinear( AstMapping *this_mapping, int *status ) {
/*
*  Name:
*     GetLinear

*  Purpose:
*     Determine if a LutMap implements a linear coordinate transformation.

*  Type:
*     Private function.

*  Synopsis:
*     #include "lutmap.h"
*     int GetLinear( AstMapping *this, int *status )

*  Class Membership:
*     LutMap member function.

*  Description:
*     This function returns a boolean value to indicate if the LutMap
*     supplied is equivalent to a linear coordinate
*     transformation. This will be the case if the lookup table
*     elements increase or decrease linearly.

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

*  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: */
   AstLutMap *this;              /* Pointer to the LutMap structure */
   double *lut;                  /* Pointer to the lookup table */
   double eps;                   /* Relative error on the table values */
   double fract;                 /* Fractional position within table */
   double hi;                    /* Largest value */
   double interp;                /* Interpolated value */
   double lo;                    /* Smallest value */
   double tol1;                  /* First tolerance estimate */
   double tol2;                  /* Second tolerance estimate */
   double tol;                   /* Tolerance value used */
   int ilut;                     /* Loop counter for table elements */
   int linear;                   /* Result to be returned */
   int nlut;                     /* Number of lookup table elements */

/* Initialise. */
   linear = 0;

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

/* Obtain a pointer to the LutMap structure. */
   this = (AstLutMap *) this_mapping;

/* Nearest neighbour LutMaps are not considered to be linear because of
   the discontinuities at the start and end of the table. */
   if( astGetLutInterp( this ) != NEAR ) {

/* Obtain the lookup table details. */
      lut = this->lut;
      nlut = this->nlut;

/* Loop to identify the largest and smallest values in the lookup
      table. */
      lo = DBL_MAX;
      hi = -DBL_MAX;
      for ( ilut = 0; ilut < nlut; ilut++ ) {
         if ( lut[ ilut ] > hi ) hi = lut[ ilut ];
         if ( lut[ ilut ] < lo ) lo = lut[ ilut ];
      }

/* Check if the values are all the same (this makes the LutMap
      linear, although it will have no inverse). */
      linear = ( hi == lo );
      if ( !linear ) {

/* Get the relative error associated with the table values. */
         eps = astGetLutEpsilon( this );

/* Form a tolerance estimate based on the overall range of values in
   the lookup table. */
         tol1 = fabs( hi - lo ) * eps;

/* Now loop to inspect all the lookup table elements except the first
      and last. */
         linear = 1;
         for ( ilut = 1; ilut < ( nlut - 1 ); ilut++ ) {

/* Calculate the fractional position of the current element within the
      table. */
            fract = ( (double) ilut ) / ( (double) ( nlut - 1 ) );

/* Calculate the value it should have if the table is linear by
      interpolating between the first and last values. */
            interp = lut[ 0 ] * ( 1.0 - fract ) + lut[ nlut - 1 ] * fract;

/* Form a second tolerance estimate from this interpolated
   value. Select whichever tolerance estimate is larger (this avoids
   problems when values are near zero). */
            tol2 = fabs( interp ) * eps;
            tol = ( tol1 > tol2 ) ? tol1 : tol2;

/* Test for linearity within a small multiple of the tolerance. */
            linear = ( fabs( lut[ ilut ] - interp ) <= ( 2.0 * tol ) );
            if ( !linear ) break;
         }
      }
   }

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

static double *GetLutMapInfo( AstLutMap *this, double *start, double *inc,
                              int *nlut, int *status ){
/*
*  Name:
*     GetLutMapInfo

*  Purpose:
*     Return information about a LutMap.

*  Type:
*     Protected virtual function.

*  Synopsis:
*     #include "permmap.h"
*     double *astGetLutMapInfo( AstLutMap *this, double *start, double *inc,
*                               int *nlut, int *status )

*  Class Membership:
*     LutMap method

*  Description:
*     This function returns information about the supplied LutMap.

*  Parameters:
*     this
*        Pointer to the LutMap.
*     start
*        Pointer to a double in which to return the "start" value
*        supplied when the LutMap was created.
*     inc
*        Pointer to a double in which to return the "inc" value
*        supplied when the LutMap was created.
*     nlut
*        Pointer to a double in which to return the number of values in
*        the LUT.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     A pointer to a dynamically allocated array holding a copy of the
*     look-up table. This is an array of "nlut" elements, giving the
*     output values for input values "start", "start+inc", "start+2*inc",
*     etc. The pointer should be freed using astFree when no longer
*     needed.

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

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

/* Store the scalar values. */
   *start = this->start;
   *inc = this->inc;
   *nlut = this->nlut;

/* Return a copy of the look up table. */
   return astStore( NULL, this->lut, sizeof( double )*this->nlut );
}

static int GetMonotonic( int nlut, const double *lut, int *nluti, double **luti,
                         int **flagsi, int **indexi, int *status ) {
/*
*  Name:
*     GetMonotonic

*  Purpose:
*     Determine if a array is monotonic increasing or decreasing.

*  Type:
*     Private function.

*  Synopsis:
*     #include "lutmap.h"
*     int GetMonotonic( int nlut, const double *lut, int *nluti, double **luti,
*                       int **flagsi, int **indexi, int *status )

*  Class Membership:
*     LutMap member function.

*  Description:
*     This function returns a flag indiciating the supplied array is
*     monotonic increasing, monotonic decreasing, or non-monotonic.
*     Sections of equal or bad values do not invalidate an otherwise
*     monotonic array.
*
*     It also returns information needed to implement the inverse
*     transformation of a LutMap.

*  Parameters:
*     nlut
*        The length of the array.
*     lut
*        The array to check.
*     nluti
*        Address of an int in which to store the length of the returned
*        "luti" and "flags" arrays.
*     luti
*        Address at which to return a pointer to a newly allocated array
*        containing "*nluti" elements. This is a copy of the supplied
*        "lut" array but with any bad or NaN values omitted. Subsequent
*        elements are shuffled down to fill the holes left by removing
*        these bad values. A NULL pointer is returned if there are no bad
*        values in "lut".
*     flagsi
*        Address at which to return a pointer to a newly allocated array
*        containing "*nluti" elements. Each element is non-zero if the
*        corresponding value stored in "luti" was adjacent to a bad value
*        in the supplied "lut" array. A NULL pointer is returned if there
*        are no bad values in "lut".
*     indexi
*        Address at which to return a pointer to a newly allocated array
*        containing "*nluti" elements. Each element is the index within
*        "lut" of the corresponding value stored in "luti". A NULL pointer
*        is returned if there are no bad values in "lut".
*     status
*        Pointer to the inherited status variable.

*  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: */
   const double *p3;
   double *p1;
   double lval;
   int *p2;
   int *p4;
   int ilut;
   int nbad;
   int result;

/* Initialise. */
   result = 0;
   *nluti = 0;
   *luti = NULL;
   *flagsi = NULL;
   *indexi = NULL;

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

/* As yet we have not found a good value. */
   lval = AST__BAD;

/* As yet we have not found a bad value. */
   nbad = 0;

/* Loop  round the supplied array, ignoring bad (AST__BAD or NaN) values. */
   for ( ilut = 0; ilut < nlut; ilut++ ) {
      if( !astISBAD( lut[ ilut ] ) ) {

/* If this is the first good value, record it. */
         if( lval == AST__BAD ) {
            lval = lut[ ilut ];

/* If this is not the first good value, ignore it if it is equal to the
   previous god value. */
         } else if( lut[ ilut ] != lval ) {

/* Set the returned flag on the basis of the first pair of good values. */
            if( result == 0 ) {
               result = ( lut[ ilut ] > lval ) ? 1 : -1;

/* For subsequent pairs of good values, check the pair increases or
   decreases in the same way as the first pair. Reset the returned value
   to zero and break if not. */
            } else if( result == 1 && lut[ ilut ] < lval ) {
               result = 0;
               break;

            } else if( result == -1 && lut[ ilut ] >lval ) {
               result = 0;
               break;
            }

/* Record the current good value. */
            lval = lut[ ilut ];
         }
      } else {
         nbad++;
      }
   }

/* If any bad values were found, we now allocate the required returned
   arrays. */
   if( nbad ) {
      *nluti = nlut - nbad;
      *luti = astMalloc( sizeof( double )*( *nluti ) );
      *flagsi = astMalloc( sizeof( double )*( *nluti ) );
      *indexi = astMalloc( sizeof( double )*( *nluti ) );

      if( astOK ) {

/* Into "luti" copy all good values from "lut", shuffling down values to
   fill holes left by bad values. Into "flagsi", store a flag indicating
   if the corresponding "luti" value was adjacent to a bad value in the
   full "lut" array. */
         p1 = *luti;
         p2 = *flagsi;
         p3 = lut;
         p4 = *indexi;

/* Do the first input point (it has no lower neighbour). */
         if( !astISBAD( *p3 ) ) {
            *(p1++) = *p3;
            *(p2++) = astISBAD( p3[ +1 ] );
            *(p4++) = 0;
         }

/* Do all remaining input points except for the last one. */
         for ( ilut = 1,p3++; ilut < nlut-1; ilut++,p3++ ) {
            if( !astISBAD( *p3 ) ) {
               *(p1++) = *p3;
               *(p2++) = astISBAD( p3[ -1 ] ) || astISBAD( p3[ +1 ] );
               *(p4++) = ilut;
            }
         }

/* Do the last input point (it has no upper neighbour). */
         if( !astISBAD( *p3 ) ) {
            *p1 = *p3;
            *p2 = astISBAD( p3[ -1 ] );
            *p4 = ilut;
         }
      }
   }


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

void astInitLutMapVtab_(  AstLutMapVtab *vtab, const char *name, int *status ) {
/*
*+
*  Name:
*     astInitLutMapVtab

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

*  Type:
*     Protected function.

*  Synopsis:
*     #include "lutmap.h"
*     void astInitLutMapVtab( AstLutMapVtab *vtab, const char *name )

*  Class Membership:
*     LutMap vtab initialiser.

*  Description:
*     This function initialises the component of a virtual function
*     table which is used by the LutMap 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 astIsALutMap) 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->ClearLutInterp = ClearLutInterp;
   vtab->GetLutInterp = GetLutInterp;
   vtab->SetLutInterp = SetLutInterp;
   vtab->TestLutInterp = TestLutInterp;
   vtab->ClearLutEpsilon = ClearLutEpsilon;
   vtab->GetLutEpsilon = GetLutEpsilon;
   vtab->SetLutEpsilon = SetLutEpsilon;
   vtab->TestLutEpsilon = TestLutEpsilon;
   vtab->GetLutMapInfo = GetLutMapInfo;

/* 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_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;

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

/* Declare the class dump, copy and delete functions.*/
   astSetDump( vtab, Dump, "LutMap",
               "Map 1-d coordinates using a lookup table" );
   astSetCopy( (AstObjectVtab *) vtab, Copy );
   astSetDelete( (AstObjectVtab *) vtab, Delete );

/* 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 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 LutMap.

*  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:
*     LutMap 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 LutMap 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 LutMap with one 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 LutMap which is to be merged with
*        its neighbours. This should be a cloned copy of the LutMap
*        pointer contained in the array element "(*map_list)[where]"
*        (see below). This pointer will not be annulled, and the
*        LutMap it identifies will not be modified by this function.
*     where
*        Index in the "*map_list" array (below) at which the pointer
*        to the nominated LutMap 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: */
   AstLutMap *map;               /* Pointer to LutMap */
   AstLutMap *neb;               /* Pointer to neighbouring LutMap */
   AstMapping *new;              /* Pointer to replacement Mapping */
   double a1;                    /* First input coordinate value */
   double a2;                    /* Second input coordinate value */
   double b1;                    /* First output coordinate value */
   double b2;                    /* Second output coordinate value */
   int equal;                    /* Are LutMaps equal? */
   int i;                        /* Mapping index */
   int ilo;                      /* Index of lower LutMap */
   int invneb;                   /* Should the neigbour be used inverted? */
   int old_inv;                  /* Original Invert value for neigbour */
   int result;                   /* Result value to return */
   int simpler;                  /* Mapping simplified? */

/* Initialise the returned result. */
   result = -1;

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

/* Obtain a pointer to the nominated LutMap. */
   map = (AstLutMap *) ( *map_list )[ where ];

/* See if the LutMap is linear. If so, it can probably be
   simplified. */
   simpler = GetLinear( (AstMapping *) map, status );
   if ( simpler ) {

/* Obtain the range of input values corresponding to the first and
   last lookup table elements. */
      a1 = map->start;
      a2 = map->start + map->inc * ( map->nlut - 1 );

/* Obtain the corresponding range of output values and check these
   values are not the same. */
      b1 = map->lut[ 0 ];
      b2 = map->lut[ map->nlut - 1 ];
      if ( b1 != b2 ) {

/* Create a new WinMap that implements an equivalent linear Mapping,
   allowing for the invert flag associated with the LutMap. */
         if ( !( *invert_list )[ where ] ) {
            new = (AstMapping *) astWinMap( 1, &a1, &a2, &b1, &b2, "", status );
         } else {
            new = (AstMapping *) astWinMap( 1, &b1, &b2, &a1, &a2, "", status );
         }

/* If OK, annul the original LutMap pointer and substitute the new
   one. Also clear the associated invert flag. */
         if ( astOK ) {
            (void) astAnnul( ( *map_list )[ where ] );
            ( *map_list )[ where ] = new;
            ( *invert_list )[ where ] = 0;

/* Assign the result value. */
            result = where;
         }
      }

/* Otherwise, see if the LutMap is in series with its own inverse. If so
   the pair of LutMaps can be replaced by a UnitMap. */
   } else if( series ) {

/* Is the higher neighbour a LutMap? If so get a pointer to it, and
   note the index of the lower of the two adjacent LutMaps. */
      if( where < ( *nmap - 1 ) &&
          astIsALutMap( ( *map_list )[ where + 1 ] ) ){
         neb = (AstLutMap *) ( *map_list )[ where + 1 ];
         invneb = ( *invert_list )[ where + 1 ];
         ilo = where;

/* If not, is the lower neighbour a LutMap? If so get a pointer to it,
   and note the index of the lower of the two adjacent LutMaps. */
      } else if( where > 0 &&
                 astIsALutMap( ( *map_list )[ where - 1 ] ) ){
         neb = (AstLutMap *) ( *map_list )[ where - 1 ];
         invneb = ( *invert_list )[ where - 1 ];
         ilo =  where - 1;

      } else {
         neb = NULL;
      }

/* If a neighbouring LutMap was found, we can replace the pair by a
   UnitMap if the two LutMaps are equal but have opposite values for
   their Invert flags. Temporarily invert the neighbour, then compare
   the two LutMaps for equality, then re-invert the neighbour. */
      if( neb ) {
         old_inv = astGetInvert( neb );
         astSetInvert( neb, invneb );
         astInvert( neb );
         equal = astEqual( map, neb );
         astSetInvert( neb, old_inv );

/* If the two LutMaps are equal but opposite, annul the first of the two
   Mappings, and replace it with a UnitMap. Also set the invert flag. */
         if( equal ) {
            new = (AstMapping *) astUnitMap( 1, "", status );
            (void) astAnnul( ( *map_list )[ ilo ] );
            ( *map_list )[ ilo ] = new;
            ( *invert_list )[ ilo ] = 0;

/* Annul the second of the two Mappings, and shuffle down the rest of the
   list to fill the gap. */
            (void) astAnnul( ( *map_list )[ ilo + 1 ] );
            for ( i = ilo + 2; i < *nmap; i++ ) {
               ( *map_list )[ i - 1 ] = ( *map_list )[ i ];
               ( *invert_list )[ i - 1 ] = ( *invert_list )[ i ];
            }

/* Clear the vacated element at the end. */
            ( *map_list )[ *nmap - 1 ] = NULL;
            ( *invert_list )[ *nmap - 1 ] = 0;

/* Decrement the Mapping count and return the index of the first
   modified element. */
            ( *nmap )--;
            result = where;
         }
      }
   }

/* If an error occurred, clear the returned result. */
   if ( !astOK ) result = -1;

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

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

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

*  Type:
*     Private function.

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

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

*  Description:
*     This function assigns an attribute value for a LutMap, 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 LutMap.
*     setting
*        Pointer to a null-terminated string specifying the new attribute
*        value.
*/

/* Local Variables: */
   AstLutMap *this;              /* Pointer to the LutMap structure */
   double luteps;                /* LutEpsilon attribute value */
   int lutinterp;                /* LutInterp 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 LutMap structure. */
   this = (AstLutMap *) 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. */

/* LutInterp. */
/* ---------- */
   if ( nc = 0,
        ( 1 == astSscanf( setting, "lutinterp= %d %n", &lutinterp, &nc ) )
        && ( nc >= len ) ) {
      astSetLutInterp( this, lutinterp );

/* LutEpsilon. */
/* ------------- */
   } else if ( nc = 0,
        ( 1 == astSscanf( setting, "lutepsilon= %lf %n", &luteps, &nc ) )
        && ( nc >= len ) ) {
      astSetLutEpsilon( this, luteps );

/* 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 int TestAttrib( AstObject *this_object, const char *attrib, int *status ) {
/*
*  Name:
*     TestAttrib

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

*  Type:
*     Private function.

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

*  Class Membership:
*     LutMap 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 LutMap's attributes.

*  Parameters:
*     this
*        Pointer to the LutMap.
*     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: */
   AstLutMap *this;             /* Pointer to the LutMap structure */
   int result;                   /* Result value to return */

/* Initialise. */
   result = 0;

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

/* Obtain a pointer to the LutMap structure. */
   this = (AstLutMap *) this_object;

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

/* LutInterp. */
/* ---------- */
   if ( !strcmp( attrib, "lutinterp" ) ) {
      result = astTestLutInterp( this );

/* LutEpsilon. */
/* ------------- */
   } else if ( !strcmp( attrib, "lutepsilon" ) ) {
      result = astTestLutEpsilon( 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 LutMap to transform a set of points.

*  Type:
*     Private function.

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

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

*  Description:
*     This function takes a LutMap and a set of points encapsulated
*     in a PointSet and transforms the points so as to apply the
*     lookup table transformation.

*  Parameters:
*     this
*        Pointer to the LutMap.
*     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 equal 1.
*     - If an output PointSet is supplied, it must have space for
*     sufficient number of points (with 1 coordinate value per point)
*     to accommodate the result. Any excess space will be ignored.
*/

/* Local Variables: */
   AstLutMap *map;               /* Pointer to LutMap to be applied */
   AstPointSet *result;          /* Pointer to output PointSet */
   double **ptr_in;              /* Pointer to input coordinate data */
   double **ptr_out;             /* Pointer to output coordinate data */
   double *lut;                  /* Pointer to LUT */
   double d1;                    /* Offset to I1 value */
   double d2;                    /* Offset to I2 value */
   double fract;                 /* Fractional interpolation distance */
   double scale;                 /* Normalising scale factor */
   double value_in;              /* Input coordinate value */
   double value_out;             /* Output coordinate value */
   double x;                     /* Value normalised to LUT increment */
   double xi;                    /* Integer value of "x" */
   int *flags;                   /* Flags indicating an adjacent bad value */
   int *index;                   /* Translates reduced to original indices */
   int i1;                       /* Lower adjacent LUT index */
   int i2;                       /* Upper adjacent LUT index */
   int i;                        /* New LUT index */
   int istart;                   /* Original LUT index at start of interval */
   int ix;                       /* "x" converted to an int */
   int near;                     /* Perform nearest neighbour interpolation? */
   int nlut;                     /* Number of LUT entries */
   int nlutm1;                   /* Number of LUT entries minus one */
   int npoint;                   /* Number of points */
   int ok;                       /* Lookup table is not flat */
   int point;                    /* Loop counter for points */
   int up;                       /* LUT values are increasing? */

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

/* Obtain a pointer to the LutMap. */
   map = (AstLutMap *) 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 );

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

/* Determine the numbers of points from the input PointSet and obtain
   pointers for accessing the input and output coordinate values. */
   npoint = astGetNpoint( in );
   ptr_in = astGetPoints( in );
   ptr_out = astGetPoints( result );

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

/* Forward transformation. */
/* ----------------------- */
   if( astOK ){
      if ( forward ) {

/* Obtain lookup table details. */
         lut = map->lut;
         nlut = map->nlut;
         near = ( astGetLutInterp( map ) == NEAR );
         nlutm1 = nlut - 1;

/* Calculate the scale factor required. */
         scale = 1.0 / map->inc;

/* Loop to transform each input point. */
         for ( point = 0; point < npoint; point++ ) {

/* Extract the input coordinate value. */
            value_in = ptr_in[ 0 ][ point ];

/* First check if this is the same value as we transformed last. If
   so, re-use the last result. */
            if ( value_in == map->last_fwd_in ) {
               value_out = map->last_fwd_out;

/* Check for bad input coordinates and generate a bad result if
   necessary. */
            } else if ( value_in == AST__BAD ) {
               value_out = AST__BAD;

/* For nearest-neighbour interpolation, return the value of the lookup table
   entry corresponding to the input coordinate. */
            } else if( near ){
               x = ( value_in - map->start ) * scale;
               xi = floor( x + 0.5 );
               ix = (int) xi;
               if ( ix < 0 || ix >= nlut ) {
                  value_out = AST__BAD;
               } else {
                  value_out = lut[ ix ];
               }

/* Otherwise, (for linear interpolation) identify the lookup table entry
   corresponding to the input coordinate. */
            } else {
               x = ( value_in - map->start ) * scale;
               xi = floor( x );
               ix = (int) xi;

/* If the input value lies below the first lookup table entry,
   extrapolate using the first two table values. */
               if ( ix < 0 ) {
                  if( lut[ 0 ] != AST__BAD && lut[ 1 ] != AST__BAD ) {
                     value_out = lut[ 0 ] + x * ( lut[ 1 ] - lut[ 0 ] );
                  } else {
                     value_out = AST__BAD;
                  }

/* If the input value lies above the last lookup table entry (or equals
   it), extrapolate using the last two table values. */
               } else if ( ix >= nlutm1 ) {
                  if( lut[ nlutm1 ] != AST__BAD &&
                      lut[ nlut - 2 ] != AST__BAD ) {
                     value_out = lut[ nlutm1 ] +
                                 ( x - (double) ( nlutm1 ) ) *
                                 ( lut[ nlutm1 ] - lut[ nlut - 2 ] );
                  } else {
                     value_out = AST__BAD;
                  }

/* Otherwise, interpolate between the adjacent entries. */
               } else {
                  if( lut[ ix ] != AST__BAD &&
                      lut[ ix + 1 ] != AST__BAD ) {
                     fract = x - xi;
                     value_out = lut[ ix ] * ( 1.0 - fract ) +
                                 lut[ ix + 1 ] * fract;
                  } else {
                     value_out = AST__BAD;
                  }
               }
            }

/* Assign the output coordinate value. */
            ptr_out[ 0 ][ point ] = value_out;

/* Retain the input and output coordinate values for possible re-use
   in future. */
            map->last_fwd_in = value_in;
            map->last_fwd_out = value_out;
         }

/* Inverse transformation. */
/* ----------------------- */
      } else {

/* Obtain details of the lookup table to be used by the inverse
   transformation. This is the same as the forward transformation lookup
   table, except that any bad values are omitted. Also, get a pointer to a
   array of flags that indicate if the corresponding lookup table entries
   were adjacent to a bad value or not in the full lookup table. */
         if( map->luti ) {
            lut = map->luti;
            flags = map->flagsi;
            nlut = map->nluti;
            index = map->indexi;
         } else {
            lut = map->lut;
            flags = NULL;
            nlut = map->nlut;
            index = NULL;
         }
         near = ( astGetLutInterp( map ) == NEAR );
         nlutm1 = nlut - 1;

/* Loop to transform each input point. */
         for ( point = 0; point < npoint; point++ ) {

/* Extract the input coordinate value. */
            value_in = ptr_in[ 0 ][ point ];

/* First check if this is the same value as we transformed last. If
   so, re-use the last result. */
            if ( value_in == map->last_inv_in ) {
               value_out = map->last_inv_out;

/* Check for bad input coordinates and generate a bad result if
   necessary. */
            } else if ( value_in == AST__BAD ) {
               value_out = AST__BAD;

/* Otherwise, we can determine an inverse. Note the inverse transformation
   will not be defined, so will not be attempted, unless all the table
   entries are monotonically increasing or decreasing, possibly with sections
   of equal or bad values. */
            } else {
               up = ( lut[ nlutm1 ] > lut[ 0 ] );

/* Perform a binary search to identify two adjacent lookup table
   elements whose values bracket the input coordinate value. */
               i1 = -1;
               i2 = nlutm1;
               while ( i2 > ( i1 + 1 ) ) {
                  i = ( i1 + i2 ) / 2;
                  *( ( ( value_in >= lut[ i ] ) == up ) ? &i1 : &i2 ) = i;
               }

/* If the lower table value is equal to the required value, and either of
   its neighbours is also equal to the required value, then we have been
   asked to find the inverse in a flat region of the table, so return
   a bad value. Likewise, if the upper table value is equal to the required
   value, and either of its neighbours is also equal to the required value,
   then we have been asked to find the inverse in a flat region of the table,
   so return a bad value. */
               ok = 1;
               if( lut[ i1 ] == value_in ) {
                  if( i1 > 0 && lut[ i1 - 1 ] == value_in ) ok = 0;
                  if( lut[ i2 ] == value_in ) ok = 0;
               } else if( lut[ i2 ] == value_in ) {
                  if( i2 < nlutm1 && lut[ i2 + 1 ] == value_in ) ok = 0;
                  if( lut[ i1 ] == value_in ) ok = 0;
               }

               if( !ok ) {
                  value_out = AST__BAD;

/* If both of the two table elements were adjacent to a bad value in the
   full lookup table, return a bad output value. */
               } else if( flags && ( flags[ i1 ] && flags[ i2 ] ) ) {
                  value_out = AST__BAD;

/* Nearest neighbour interpolation: return the closest of i1 or i2. Return
   AST__BAD if the supplied value is less than either or greater than
   either.  */
               } else if( near ) {
                  d1 = lut[ i1 ] - value_in;
                  d2 = lut[ i2 ] - value_in;
                  if( ( d1 > 0.0 && d2 > 0.0 ) ||
                      ( d1 < 0.0 && d2 < 0.0 ) ) {
                     value_out = AST__BAD;

                  } else {

                     if( fabs( d1 ) < fabs( d2 ) ){
                        istart = index ? index[ i1 ] : i1;
                     } else {
                        istart = index ? index[ i2 ] : i2;
                     }
                     value_out = map->start + map->inc * istart;

                  }

/* Linear interpolation... */
               } else {

/* We are interested in the lower bracketing table element. If
   necessary, restrict this element's index to lie within the
   table. This causes extrapolation to occur (instead of
   interpolation) if the input value actually lies outside the range
   of the lookup table. */
                  if ( i1 < 0 ) i1 = 0;
                  if ( i1 > ( nlut - 2 ) ) i1 = nlut - 2;

/* Interpolate (or extrapolate) to derive the output coordinate
   value. */
                  istart = index ? index[ i1 ] : i1;
                  value_out = map->start + map->inc * ( (double) istart +
                                             ( ( value_in - lut[ i1 ] ) /
                                             ( lut[ i1 + 1 ] - lut[ i1 ] ) ) );
               }
            }

/* Assign the output coordinate value. */
            ptr_out[ 0 ][ point ] = value_out;

/* Retain the input and output coordinate values for possible re-use
   in future. */
            map->last_inv_in = value_in;
            map->last_inv_out = value_out;
         }
      }
   }

/* 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:
*     LutInterp

*  Purpose:
*     Look-up table interpolation method.

*  Type:
*     Public attribute.

*  Synopsis:
*     Integer.

*  Description:
*     This attribute indicates the method to be used when finding the
*     output value of a LutMap for an input value part way between two
*     table entries. If it is set to 0 (the default) then linear
*     interpolation is used. Otherwise, nearest neighbour interpolation
*     is used.
*
*     Using nearest neighbour interpolation causes AST__BAD to be returned
*     for any point which falls outside the bounds of the table. Linear
*     interpolation results in an extrapolated value being returned based
*     on the two end entries in the table.
*
*     Note, the value of this attribute may changed only if the LutMap
*     has no more than one reference. That is, an error is reported if the
*     LutMap has been cloned, either by including it within another object
*     such as a CmpMap or FrameSet or by calling the
c     astClone
f     AST_CLONE
*     function.

*  Applicability:
*     LutMap
*        All LutMaps have this attribute.

*att--
*/
astMAKE_CLEAR1(LutMap,LutInterp,lutinterp,-INT_MAX)
astMAKE_GET(LutMap,LutInterp,int,LINEAR,( ( this->lutinterp == -INT_MAX ) ?
                                          LINEAR : this->lutinterp ))
astMAKE_SET1(LutMap,LutInterp,int,lutinterp,(( value == LINEAR ) ? LINEAR : NEAR ))
astMAKE_TEST(LutMap,LutInterp,( this->lutinterp != -INT_MAX ))

/*
*att++
*  Name:
*     LutEpsilon

*  Purpose:
*     The relative error of the values held in the took-up table.

*  Type:
*     Public attribute.

*  Synopsis:
*     Double precision.

*  Description:
*     This attribute holds the relative error of the values held in the
*     took-up table. It is used when simplifying a LutMap, to determine
*     if the LutMap should be considered linear. Setting a larger value
*     makes it more likely that a LutMap will be replaced by a WinMap
*     (i.e. a linear Mapping) when simplified.
*
*     The default value is the value of the system constant DBL_EPSILON
*     (typically around 1e-16 or 2E-16). If the values in the look-up
*     table were derived from single precision data, it may be appropriate
*     to set this attribute to a value around 1E-7.
*
*     Note, the value of this attribute may changed only if the LutMap
*     has no more than one reference. That is, an error is reported if the
*     LutMap has been cloned, either by including it within another object
*     such as a CmpMap or FrameSet or by calling the
c     astClone
f     AST_CLONE
*     function.

*  Applicability:
*     LutMap
*        All LutMaps have this attribute.

*att--
*/
astMAKE_CLEAR1(LutMap,LutEpsilon,lutepsilon,AST__BAD)
astMAKE_GET(LutMap,LutEpsilon,double,DBL_EPSILON,( ( this->lutepsilon == AST__BAD ) ?
                                          DBL_EPSILON : this->lutepsilon ))
astMAKE_SET1(LutMap,LutEpsilon,double,lutepsilon,(value))
astMAKE_TEST(LutMap,LutEpsilon,( this->lutepsilon != AST__BAD ))

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

*  Purpose:
*     Copy constructor for LutMap objects.

*  Type:
*     Private function.

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

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

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

/* Local Variables: */
   AstLutMap *out;               /* Pointer to output LutMap */
   AstLutMap *in;                /* Pointer to input LutMap */

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

/* Obtain a pointer to the input and output LutMaps. */
   in= (AstLutMap *) objin;
   out = (AstLutMap *) objout;

/* Nullify all output pointers. */
   out->lut = NULL;
   out->luti = NULL;
   out->flagsi = NULL;
   out->indexi = NULL;

/* Allocate memory and store a copy of the lookup table data. */
   out->lut = astStore( NULL, in->lut,
                        sizeof( double ) * (size_t) in->nlut );

/* Do the arrays used for the inverse transformation, if they exist. */
   if( in->luti ) out->luti = astStore( NULL, in->luti,
                                        sizeof( double ) * (size_t) in->nluti );
   if( in->flagsi ) out->flagsi = astStore( NULL, in->flagsi,
                                        sizeof( double ) * (size_t) in->nluti );
   if( in->indexi ) out->indexi = astStore( NULL, in->indexi,
                                        sizeof( double ) * (size_t) in->nluti );
}

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

*  Purpose:
*     Destructor for LutMap objects.

*  Type:
*     Private function.

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

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

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

/* Local Variables: */
   AstLutMap *this;              /* Pointer to LutMap */

/* Obtain a pointer to the LutMap structure. */
   this = (AstLutMap *) obj;

/* Free the memory holding the lookup tables, etc. */
   this->lut = astFree( this->lut );
   this->luti = astFree( this->luti );
   this->flagsi = astFree( this->flagsi );
   this->indexi = astFree( this->indexi );
}

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

*  Purpose:
*     Dump function for LutMap 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 LutMap class to an output Channel.

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

/* Local Constants: */
#define KEY_LEN 50               /* Maximum length of a keyword */

/* Local Variables: */
   AstLutMap *this;              /* Pointer to the LutMap structure */
   char buff[ KEY_LEN + 1 ];     /* Buffer for keyword string */
   double dval;                  /* Double value */
   int ilut;                     /* Loop counter for table elements */
   int ival;                     /* Integer value */
   int set;                      /* Attribute value set? */

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

/* Obtain a pointer to the LutMap structure. */
   this = (AstLutMap *) this_object;

/* Write out values representing the instance variables for the LutMap
   class.  Accompany these with appropriate comment strings, possibly
   depending on the values being written. */

/* Number of lookup table elements. */
   astWriteInt( channel, "Nlut", 1, 1, this->nlut,
                "Number of lookup table elements" );

/* Input coordinate at first element centre. */
   astWriteDouble( channel, "Start", ( this->start != 0.0 ), 1, this->start,
                   "Input value at first element" );

/* Element spacing. */
   astWriteDouble( channel, "Incr", ( this->inc != 1.0 ), 1, this->inc,
                   "Input value increment between elements" );

/* Interpolation method */
   set = TestLutInterp( this, status );
   ival = set ? GetLutInterp( this, status ) : astGetLutInterp( this );
   astWriteInt( channel, "LutInt", set, 1, ival, "Interpolation method" );

/* Precision */
   if( TestLutEpsilon( this, status ) ) {
      dval = GetLutEpsilon( this, status );
      astWriteDouble( channel, "LutEps", 1, 1, dval, "Table relative error" );
   }

/* Lookup table contents. */
   for ( ilut = 0; ilut < this->nlut; ilut++ ) {
      if( this->lut[ ilut ] != AST__BAD ) {
         (void) sprintf( buff, "L%d", ilut + 1 );
         astWriteDouble( channel, buff, 1, 1, this->lut[ ilut ],
                         ilut ? "" : "Lookup table elements..." );
      }
   }

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

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

AstLutMap *astLutMap_( int nlut, const double lut[],
                       double start, double inc,
                       const char *options, int *status, ...) {
/*
*++
*  Name:
c     astLutMap
f     AST_LUTMAP

*  Purpose:
*     Create a LutMap.

*  Type:
*     Public function.

*  Synopsis:
c     #include "lutmap.h"
c     AstLutMap *astLutMap( int nlut, const double lut[],
c                           double start, double inc,
c                           const char *options, ... )
f     RESULT = AST_LUTMAP( NLUT, LUT, START, INC, OPTIONS, STATUS )

*  Class Membership:
*     LutMap constructor.

*  Description:
*     This function creates a new LutMap and optionally initialises
*     its attributes.
*
*     A LutMap is a specialised form of Mapping which transforms
*     1-dimensional coordinates by using linear interpolation in a
*     lookup table.  Each input coordinate value is first scaled to
*     give the index of an entry in the table by subtracting a
*     starting value (the input coordinate corresponding to the first
*     table entry) and dividing by an increment (the difference in
*     input coordinate value between adjacent table entries).
*
*     The resulting index will usually contain a fractional part, so
*     the output coordinate value is then generated by interpolating
*     linearly between the appropriate entries in the table. If the
*     index lies outside the range of the table, linear extrapolation
*     is used based on the two nearest entries (i.e. the two entries
*     at the start or end of the table, as appropriate).
*
*     If the lookup table entries increase or decrease monotonically,
*     then the inverse transformation may also be performed.

*  Parameters:
c     nlut
f     NLUT = INTEGER (Given)
*        The number of entries in the lookup table. This value must be
*        at least 2.
c     lut
f     LUT( NLUT ) = DOUBLE PRECISION (Given)
c        An array containing the "nlut"
f        An array containing the
*        lookup table entries.
c     start
f     START = DOUBLE PRECISION (Given)
*        The input coordinate value which corresponds to the first lookup
*        table entry.
c     inc
f     INC = DOUBLE PRECISION (Given)
*        The lookup table spacing (the increment in input coordinate
*        value between successive lookup table entries). This value
*        may be positive or negative, but must not be zero.
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 LutMap. 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 LutMap. 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     astLutMap()
f     AST_LUTMAP = INTEGER
*        A pointer to the new LutMap.

*  Notes:
*     - If the entries in the lookup table either increase or decrease
*     monotonically, then the new LutMap's TranInverse attribute will
*     have a value of one, indicating that the inverse transformation
*     can be performed. Otherwise, it will have a value of zero, so
*     that any attempt to use the inverse transformation will result
*     in an error.
*     - 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.

*  Status Handling:
*     The protected interface to this function includes an extra
*     parameter at the end of the parameter list descirbed above. This
*     parameter is a pointer to the integer inherited status
*     variable: "int *status".

*--
*/

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

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

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

/* Initialise the LutMap, allocating memory and initialising the
   virtual function table as well if necessary. */
   new = astInitLutMap( NULL, sizeof( AstLutMap ), !class_init, &class_vtab,
                        "LutMap", nlut, lut, start, inc );

/* 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
   LutMap'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 LutMap. */
   return new;
}

AstLutMap *astLutMapId_( int nlut, const double lut[],
                         double start, double inc,
                         const char *options, ... ) {
/*
*  Name:
*     astLutMapId_

*  Purpose:
*     Create a LutMap.

*  Type:
*     Private function.

*  Synopsis:
*     #include "lutmap.h"
*     AstLutMap *astLutMapId( int nlut, const double lut[],
*                             double start, double inc,
*                             const char *options, ... )

*  Class Membership:
*     LutMap constructor.

*  Description:
*     This function implements the external (public) interface to the
*     astLutMap constructor function. It returns an ID value (instead
*     of a true C pointer) to external users, and must be provided
*     because astLutMap_ 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 astLutMap_ 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 astLutMap_.

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

/* Local Variables: */
   astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
   AstLutMap *new;               /* Pointer to new LutMap */
   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 LutMap, allocating memory and initialising the
   virtual function table as well if necessary. */
   new = astInitLutMap( NULL, sizeof( AstLutMap ), !class_init, &class_vtab,
                        "LutMap", nlut, lut, start, inc );

/* 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 LutMap'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 LutMap. */
   return astMakeId( new );
}

AstLutMap *astInitLutMap_( void *mem, size_t size, int init,
                           AstLutMapVtab *vtab, const char *name,
                           int nlut, const double lut[],
                           double start, double inc, int *status ) {
/*
*+
*  Name:
*     astInitLutMap

*  Purpose:
*     Initialise a LutMap.

*  Type:
*     Protected function.

*  Synopsis:
*     #include "lutmap.h"
*     AstLutMap *astInitLutMap( void *mem, size_t size, int init,
*                               AstLutMapVtab *vtab, const char *name,
*                               int nlut, const double lut[],
*                               double start, double inc )

*  Class Membership:
*     LutMap initialiser.

*  Description:
*     This function is provided for use by class implementations to
*     initialise a new LutMap object. It allocates memory (if
*     necessary) to accommodate the LutMap plus any additional data
*     associated with the derived class.  It then initialises a LutMap
*     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 LutMap at the start of the memory passed via the
*     "vtab" parameter.

*  Parameters:
*     mem
*        A pointer to the memory in which the LutMap is to be
*        initialised.  This must be of sufficient size to accommodate
*        the LutMap data (sizeof(LutMap)) 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 LutMap (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 LutMap structure, so a valid value must be
*        supplied even if not required for allocating memory.
*     init
*        A logical flag indicating if the LutMap'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 LutMap.
*     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).
*     nlut
*        The number of elements in the lookup table. This value must
*        be at least 2.
*     lut
*        An array containing the "nlut" lookup table elements.
*     start
*        The input coordinate value which corresponds with the first
*        lookup table element.
*     inc
*        The lookup table element spacing (i.e. the increment in input
*        coordinate value between successive lookup table elements).

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

*  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: */
   AstLutMap *new;       /* Pointer to new LutMap */
   double *luti;         /* Pointer to table for inverse transformation */
   double *p;            /* Pointer to next lut element */
   int *flagsi;          /* Pointer to flags for inverse transformation */
   int *indexi;          /* Pointer to translation from original to reduced */
   int dirn;             /* +1 => values increasing, -1 => values decreasing */
   int ilut;             /* Loop counter for LUT elements */
   int nluti;            /* Length of "luti" array */

/* Initialise. */
   new = NULL;

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

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

/* Check that the number of lookup table elements is valid. */
   if ( nlut < 2 ) {
      astError( AST__LUTIN, "astInitLutMap(%s): Invalid number of lookup "
                "table elements (%d).", status, name, nlut );
      astError( AST__LUTIN, "This value should be at least 2." , status);

/* Also check that the input value increment is not zero. */
   } else if ( inc == 0.0 ) {
      astError( AST__LUTII, "astInitLutMap(%s): An input value increment of "
                "zero between lookup table elements is not allowed.", status, name );

/* Determine if the element values increase or decrease monotonically (except
   that adjacent entries can be equal). We can only implement the inverse
   transformation if this is so. The inverse transformation will generate
   AST__BAD output values for sections of the table that contain equal
   adjacent values, or hold AST__BAD values. */
   } else {
      dirn = GetMonotonic( nlut, lut, &nluti, &luti, &flagsi, &indexi,
                           status );

/* Initialise a Mapping structure (the parent class) as the first
   component within the LutMap structure, allocating memory if
   necessary. Specify that the Mapping should be defined in the
   forward direction, and conditionally in the inverse direction. */
      new = (AstLutMap *) astInitMapping( mem, size, 0,
                                         (AstMappingVtab *) vtab, name,
                                          1, 1, 1, ( dirn != 0 ) );

      if ( astOK ) {

/* Initialise the LutMap data. */
/* ---------------------------- */
         new->nlut = nlut;
         new->start = start;
         new->inc = inc;
         new->lutinterp = LINEAR;
         new->lutepsilon = AST__BAD;
         new->nluti = nluti;
         new->luti = luti;
         new->flagsi = flagsi;
         new->indexi = indexi;

/* Allocate memory and store the lookup table. */
         new->lut = astStore( NULL, lut, sizeof( double ) * (size_t) nlut );

/* Replace an NaN values by AST__BAD */
         p = new->lut;
         for ( ilut = 0; ilut < nlut; ilut++, p++ ) {
            if( !astISFINITE(*p) ) *p = AST__BAD;
         }

/* Initialise the retained input and output coordinate values. */
         new->last_fwd_in = AST__BAD;
         new->last_fwd_out = AST__BAD;
         new->last_inv_in = AST__BAD;
         new->last_inv_out = AST__BAD;
      }

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

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

AstLutMap *astLoadLutMap_( void *mem, size_t size,
                           AstLutMapVtab *vtab, const char *name,
                           AstChannel *channel, int *status ) {
/*
*+
*  Name:
*     astLoadLutMap

*  Purpose:
*     Load a LutMap.

*  Type:
*     Protected function.

*  Synopsis:
*     #include "lutmap.h"
*     AstLutMap *astLoadLutMap( void *mem, size_t size,
*                               AstLutMapVtab *vtab, const char *name,
*                               AstChannel *channel )

*  Class Membership:
*     LutMap loader.

*  Description:
*     This function is provided to load a new LutMap 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
*     LutMap 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 LutMap at the start of the memory
*     passed via the "vtab" parameter.


*  Parameters:
*     mem
*        A pointer to the memory into which the LutMap is to be
*        loaded.  This must be of sufficient size to accommodate the
*        LutMap data (sizeof(LutMap)) 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 LutMap (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 LutMap 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(AstLutMap) is used instead.
*     vtab
*        Pointer to the start of the virtual function table to be
*        associated with the new LutMap. If this is NULL, a pointer
*        to the (static) virtual function table for the LutMap 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 "LutMap" is used instead.

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

*  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 Constants: */
   astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
#define KEY_LEN 50               /* Maximum length of a keyword */

/* Local Variables: */
   AstLutMap *new;               /* Pointer to the new LutMap */
   char buff[ KEY_LEN + 1 ];     /* Buffer for keyword string */
   int ilut;                     /* Get a pointer to the thread specific global data structure. */
   astGET_GLOBALS(channel);

/* Loop counter for table elements */

/* 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 LutMap. In this case the
   LutMap belongs to this class, so supply appropriate values to be
   passed to the parent class loader (and its parent, etc.). */
   if ( !vtab ) {
      size = sizeof( AstLutMap );
      vtab = &class_vtab;
      name = "LutMap";

/* If required, initialise the virtual function table for this class. */
      if ( !class_init ) {
         astInitLutMapVtab( 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 LutMap. */
   new = astLoadMapping( mem, size, (AstMappingVtab *) vtab, name,
                         channel );

   if ( astOK ) {

/* Read input data. */
/* ================ */
/* Request the input Channel to read all the input data appropriate to
   this class into the internal "values list". */
      astReadClassData( channel, "LutMap" );

/* Now read each individual data item from this list and use it to
   initialise the appropriate instance variable(s) for this class. */

/* Number of lookup table elements. */
      new->nlut = astReadInt( channel, "nlut", 2 );

/* Starting input coordinate value. */
      new->start = astReadDouble( channel, "start", 0.0 );

/* Input coordinate value increment. */
      new->inc = astReadDouble( channel, "incr", 1.0 );

/* Interpolation method */
      new->lutinterp = astReadInt( channel, "lutint", LINEAR );
      if ( TestLutInterp( new, status ) ) SetLutInterp( new, new->lutinterp, status );

/* Precision */
      new->lutepsilon = astReadDouble( channel, "luteps", AST__BAD );
      if ( TestLutEpsilon( new, status ) ) SetLutEpsilon( new, new->lutepsilon, status );

/* Allocate memory to hold the lookup table elements. */
      new->lut = astMalloc( sizeof( double ) * (size_t) new->nlut );

/* If OK, loop to read each element. */
      if ( astOK ) {
         for ( ilut = 0; ilut < new->nlut; ilut++ ) {
            (void) sprintf( buff, "l%d", ilut + 1 );
            new->lut[ ilut ] = astReadDouble( channel, buff, AST__BAD );
         }

/* Initialise the retained input and output coordinate values. */
         new->last_fwd_in = AST__BAD;
         new->last_fwd_out = AST__BAD;
         new->last_inv_in = AST__BAD;
         new->last_inv_out = AST__BAD;

/* See if the array is monotonic increasing or decreasing. */
         (void) GetMonotonic( new->nlut, new->lut, &(new->nluti),
                              &(new->luti), &(new->flagsi), &(new->indexi),
                              status );
      }
   }

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

/* Return the new LutMap 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. */

double *astGetLutMapInfo_( AstLutMap *this, double *start, double *inc,
                           int *nlut, int *status ){
   if( !astOK ) return NULL;
   return (**astMEMBER(this,LutMap,GetLutMapInfo))( this, start, inc, nlut,
                                                    status );
}