path: root/ast/mathmap.c

/*
*class++
*  Name:
*     MathMap

*  Purpose:
*     Transform coordinates using mathematical expressions.

*  Constructor Function:
c     astMathMap
f     AST_MATHMAP

*  Description:
c     A MathMap is a Mapping which allows you to specify a set of forward
c     and/or inverse transformation functions using arithmetic operations
c     and mathematical functions similar to those available in C. The
c     MathMap interprets these functions at run-time, whenever its forward
c     or inverse transformation is required. Because the functions are not
c     compiled in the normal sense (unlike an IntraMap), they may be used to
c     describe coordinate transformations in a transportable manner. A
c     MathMap therefore provides a flexible way of defining new types of
c     Mapping whose descriptions may be stored as part of a dataset and
c     interpreted by other programs.
f     A MathMap is a Mapping which allows you to specify a set of forward
f     and/or inverse transformation functions using arithmetic operations
f     and mathematical functions similar to those available in Fortran. The
f     MathMap interprets these functions at run-time, whenever its forward
f     or inverse transformation is required. Because the functions are not
f     compiled in the normal sense (unlike an IntraMap), they may be used to
f     describe coordinate transformations in a transportable manner. A
f     MathMap therefore provides a flexible way of defining new types of
f     Mapping whose descriptions may be stored as part of a dataset and
f     interpreted by other programs.

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

*  Attributes:
*     In addition to those attributes common to all Mappings, every
*     MathMap also has the following attributes:
*     - Seed: Random number seed
*     - SimpFI: Forward-inverse MathMap pairs simplify?
*     - SimpIF: Inverse-forward MathMap pairs simplify?

*  Functions:
c     The MathMap class does not define any new functions beyond those
f     The MathMap 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

*  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)

*  History:
*     3-SEP-1999 (RFWS):
*        Original version.
*     8-JAN-2003 (DSB):
*        Changed private InitVtab method to protected astInitMathMapVtab
*        method.
*     14-FEB-2006 (DSB):
*        Override astGetObjSize.
*     14-MAR-2006 (DSB):
*        - Add QIF function.
*        - Override astEqual method.
*     20-NOV-2006 (DSB):
*        Re-implement the Equal method to avoid use of astSimplify.
*     30-AUG-2012 (DSB):
*        Fix bug in undocumented Gaussian noise function.
*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 MathMap

/* Allocate pointer array. */
/* ----------------------- */
/* This macro allocates an array of pointers. If successful, each element
   of the array is initialised to NULL. */
#define MALLOC_POINTER_ARRAY(array_name,array_type,array_size) \
\
/* Allocate the array. */ \
   (array_name) = astMalloc( sizeof(array_type) * (size_t) (array_size) ); \
   if ( astOK ) { \
\
/* If successful, loop to initialise each element. */ \
      int array_index_; \
      for ( array_index_ = 0; array_index_ < (array_size); array_index_++ ) { \
         (array_name)[ array_index_ ] = NULL; \
      } \
   }

/* Free pointer array. */
/* ------------------- */
/* This macro frees a dynamically allocated array of pointers, each of
   whose elements may point at a further dynamically allocated array
   (which is also to be freed). It also allows for the possibility of any
   of the pointers being NULL. */
#define FREE_POINTER_ARRAY(array_name,array_size) \
\
/* Check that the main array pointer is not NULL. */ \
   if ( (array_name) ) { \
\
/* If OK, loop to free each of the sub-arrays. */ \
      int array_index_; \
      for ( array_index_ = 0; array_index_ < (array_size); array_index_++ ) { \
\
/* Check that each sub-array pointer is not NULL before freeing it. */ \
         if ( (array_name)[ array_index_ ] ) { \
            (array_name)[ array_index_ ] = \
               astFree( (array_name)[ array_index_ ] ); \
         } \
      } \
\
/* Free the main pointer array. */ \
      (array_name) = astFree( (array_name) ); \
   }

/* SizeOf pointer array. */
/* --------------------- */
/* This macro increments "result" by the number of bytes allocated for an
   array of pointers, each of whose elements may point at a further
   dynamically allocated array (which is also to be included). It also
   allows for the possibility of any of the pointers being NULL. */
#define SIZEOF_POINTER_ARRAY(array_name,array_size) \
\
/* Check that the main array pointer is not NULL. */ \
   if ( (array_name) ) { \
\
/* If OK, loop to measure each of the sub-arrays. */ \
      int array_index_; \
      for ( array_index_ = 0; array_index_ < (array_size); array_index_++ ) { \
\
/* Check that each sub-array pointer is not NULL before measuring it. */ \
         if ( (array_name)[ array_index_ ] ) { \
            result += astTSizeOf( (array_name)[ array_index_ ] ); \
         } \
      } \
\
/* Include the main pointer array. */ \
      result += astTSizeOf( (array_name) ); \
   }

/* Header files. */
/* ============= */
/* Interface definitions. */
/* ---------------------- */
#include "channel.h"             /* I/O channels */

#include "globals.h"             /* Thread-safe global data access */
#include "error.h"               /* Error reporting facilities */
#include "mapping.h"             /* Coordinate mappings (parent class) */
#include "cmpmap.h"              /* Compound Mappings */
#include "mathmap.h"             /* Interface definition for this class */
#include "memory.h"              /* Memory allocation facilities */
#include "globals.h"             /* Thread-safe global data access */
#include "object.h"              /* Base Object class */
#include "pointset.h"            /* Sets of points */
#include "unitmap.h"             /* Unit Mapping */

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

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

/* Module Variables. */
/* ================= */
/* This type is made obscure since it is publicly accessible (but not
   useful). Provide shorthand for use within this module. */
typedef AstMathMapRandContext_ Rcontext;



/* 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 int (* parent_getobjsize)( AstObject *, int * );
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 * );

/* This declaration enumerates the operation codes recognised by the
   EvaluateFunction function which evaluates arithmetic expressions. */
typedef enum {

/* User-supplied constants and variables. */
   OP_LDCON,                     /* Load constant */
   OP_LDVAR,                     /* Load variable */

/* System constants. */
   OP_LDBAD,                     /* Load bad value (AST__BAD) */
   OP_LDDIG,                     /* Load # decimal digits (DBL_DIG) */
   OP_LDEPS,                     /* Load relative precision (DBL_EPSILON) */
   OP_LDMAX,                     /* Load largest value (DBL_MAX) */
   OP_LDMAX10E,                  /* Max. decimal exponent (DBL_MAX_10_EXP) */
   OP_LDMAXE,                    /* Load maximum exponent (DBL_MAX_EXP) */
   OP_LDMDIG,                    /* Load # mantissa digits (DBL_MANT_DIG) */
   OP_LDMIN,                     /* Load smallest value (DBL_MIN) */
   OP_LDMIN10E,                  /* Min. decimal exponent (DBL_MIN_10_EXP) */
   OP_LDMINE,                    /* Load minimum exponent (DBL_MIN_EXP) */
   OP_LDRAD,                     /* Load floating radix (FLT_RADIX) */
   OP_LDRND,                     /* Load rounding mode (FLT_ROUNDS) */

/* Mathematical constants. */
   OP_LDE,                       /* Load e (base of natural logarithms) */
   OP_LDPI,                      /* Load pi */

/* Functions with one argument. */
   OP_ABS,                       /* Absolute value (sign removal) */
   OP_ACOS,                      /* Inverse cosine (radians) */
   OP_ACOSD,                     /* Inverse cosine (degrees) */
   OP_ACOSH,                     /* Inverse hyperbolic cosine */
   OP_ACOTH,                     /* Inverse hyperbolic cotangent */
   OP_ACSCH,                     /* Inverse hyperbolic cosecant */
   OP_ASECH,                     /* Inverse hyperbolic secant */
   OP_ASIN,                      /* Inverse sine (radians) */
   OP_ASIND,                     /* Inverse sine (degrees) */
   OP_ASINH,                     /* Inverse hyperbolic sine */
   OP_ATAN,                      /* Inverse tangent (radians) */
   OP_ATAND,                     /* Inverse tangent (degrees) */
   OP_ATANH,                     /* Inverse hyperbolic tangent */
   OP_CEIL,                      /* C ceil function (round up) */
   OP_COS,                       /* Cosine (radians) */
   OP_COSD,                      /* Cosine (degrees) */
   OP_COSH,                      /* Hyperbolic cosine */
   OP_COTH,                      /* Hyperbolic cotangent */
   OP_CSCH,                      /* Hyperbolic cosecant */
   OP_EXP,                       /* Exponential function */
   OP_FLOOR,                     /* C floor function (round down) */
   OP_INT,                       /* Integer value (round towards zero) */
   OP_ISBAD,                     /* Test for bad value */
   OP_LOG,                       /* Natural logarithm */
   OP_LOG10,                     /* Base 10 logarithm */
   OP_NINT,                      /* Fortran NINT function (round to nearest) */
   OP_POISS,                     /* Poisson random number */
   OP_SECH,                      /* Hyperbolic secant */
   OP_SIN,                       /* Sine (radians) */
   OP_SINC,                      /* Sinc function [= sin(x)/x] */
   OP_SIND,                      /* Sine (degrees) */
   OP_SINH,                      /* Hyperbolic sine */
   OP_SQR,                       /* Square */
   OP_SQRT,                      /* Square root */
   OP_TAN,                       /* Tangent (radians) */
   OP_TAND,                      /* Tangent (degrees) */
   OP_TANH,                      /* Hyperbolic tangent */

/* Functions with two arguments. */
   OP_ATAN2,                     /* Inverse tangent (2 arguments, radians) */
   OP_ATAN2D,                    /* Inverse tangent (2 arguments, degrees) */
   OP_DIM,                       /* Fortran DIM (positive difference) fn. */
   OP_GAUSS,                     /* Gaussian random number */
   OP_MOD,                       /* Modulus function */
   OP_POW,                       /* Raise to power */
   OP_RAND,                      /* Uniformly distributed random number */
   OP_SIGN,                      /* Transfer of sign function */

/* Functions with three arguments. */
   OP_QIF,                       /* C "question mark" operator "a?b:c" */

/* Functions with variable numbers of arguments. */
   OP_MAX,                       /* Maximum of 2 or more values */
   OP_MIN,                       /* Minimum of 2 or more values */

/* Unary arithmetic operators. */
   OP_NEG,                       /* Negate (change sign) */

/* Unary boolean operators. */
   OP_NOT,                       /* Boolean NOT */

/* Binary arithmetic operators. */
   OP_ADD,                       /* Add */
   OP_DIV,                       /* Divide */
   OP_MUL,                       /* Multiply */
   OP_SUB,                       /* Subtract */

/* Bit-shift operators. */
   OP_SHFTL,                     /* Shift bits left */
   OP_SHFTR,                     /* Shift bits right */

/* Relational operators. */
   OP_EQ,                        /* Relational equal */
   OP_GE,                        /* Greater than or equal */
   OP_GT,                        /* Greater than */
   OP_LE,                        /* Less than or equal */
   OP_LT,                        /* Less than */
   OP_NE,                        /* Not equal */

/* Bit-wise operators. */
   OP_BITAND,                    /* Bit-wise AND */
   OP_BITOR,                     /* Bit-wise OR */
   OP_BITXOR,                    /* Bit-wise exclusive OR */

/* Binary boolean operators. */
   OP_AND,                       /* Boolean AND */
   OP_EQV,                       /* Fortran logical .EQV. operation */
   OP_OR,                        /* Boolean OR */
   OP_XOR,                       /* Boolean exclusive OR */

/* Null operation. */
   OP_NULL                       /* Null operation */
} Oper;

/* This structure holds a description of each symbol which may appear
   in an expression. */
typedef struct {
   const char *text;             /* Symbol text as it appears in expressions */
   const int size;               /* Size of symbol text */
   const int operleft;           /* An operator when seen from the left? */
   const int operright;          /* An operator when seen from the right? */
   const int unarynext;          /* May be followed by a unary +/- ? */
   const int unaryoper;          /* Is a unary +/- ? */
   const int leftpriority;       /* Priority when seen from the left */
   const int rightpriority;      /* Priority when seen from the right */
   const int parincrement;       /* Change in parenthesis level */
   const int stackincrement;     /* Change in evaluation stack size */
   const int nargs;              /* Number of function arguments */
   const Oper opcode;            /* Resulting operation code */
} Symbol;

/* This initialises an array of Symbol structures to hold data on all
   the supported symbols. The order is not important, but symbols are
   arranged here in approximate order of descending evaluation
   priority. The end of the array is indicated by an element with a NULL
   "text" component. */
static const Symbol symbol[] = {

/* User-supplied constants and variables. */
   { ""            ,  0,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDCON    },
   { ""            ,  0,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDVAR    },

/* System constants. */
   { "<bad>"       ,  5,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDBAD    },
   { "<dig>"       ,  5,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDDIG    },
   { "<epsilon>"   ,  9,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDEPS    },
   { "<mant_dig>"  , 10,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDMDIG   },
   { "<max>"       ,  5,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDMAX    },
   { "<max_10_exp>", 12,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDMAX10E },
   { "<max_exp>"   ,  9,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDMAXE   },
   { "<min>"       ,  5,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDMIN    },
   { "<min_10_exp>", 12,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDMIN10E },
   { "<min_exp>"   ,  9,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDMINE   },
   { "<radix>"     ,  7,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDRAD    },
   { "<rounds>"    ,  8,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDRND    },

/* Mathematical constants. */
   { "<e>"         ,  3,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDE      },
   { "<pi>"        ,  4,  0,  0,  0,  0, 19, 19,  0,  1,  0,  OP_LDPI     },

/* Functions with one argument. */
   { "abs("        ,  4,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ABS      },
   { "acos("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ACOS     },
   { "acosd("      ,  6,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ACOSD    },
   { "acosh("      ,  6,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ACOSH    },
   { "acoth("      ,  6,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ACOTH    },
   { "acsch("      ,  6,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ACSCH    },
   { "aint("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_INT      },
   { "asech("      ,  6,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ASECH    },
   { "asin("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ASIN     },
   { "asind("      ,  6,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ASIND    },
   { "asinh("      ,  6,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ASINH    },
   { "atan("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ATAN     },
   { "atand("      ,  6,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ATAND    },
   { "atanh("      ,  6,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ATANH    },
   { "ceil("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_CEIL     },
   { "cos("        ,  4,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_COS      },
   { "cosd("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_COSD     },
   { "cosh("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_COSH     },
   { "coth("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_COTH     },
   { "csch("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_CSCH     },
   { "exp("        ,  4,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_EXP      },
   { "fabs("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ABS      },
   { "floor("      ,  6,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_FLOOR    },
   { "int("        ,  4,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_INT      },
   { "isbad("      ,  6,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_ISBAD    },
   { "log("        ,  4,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_LOG      },
   { "log10("      ,  6,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_LOG10    },
   { "nint("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_NINT     },
   { "poisson("    ,  8,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_POISS    },
   { "sech("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_SECH     },
   { "sin("        ,  4,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_SIN      },
   { "sinc("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_SINC     },
   { "sind("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_SIND     },
   { "sinh("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_SINH     },
   { "sqr("        ,  4,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_SQR      },
   { "sqrt("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_SQRT     },
   { "tan("        ,  4,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_TAN      },
   { "tand("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_TAND     },
   { "tanh("       ,  5,  0,  1,  1,  0, 19,  1,  1,  0,  1,  OP_TANH     },

/* Functions with two arguments. */
   { "atan2("      ,  6,  0,  1,  1,  0, 19,  1,  1, -1,  2,  OP_ATAN2    },
   { "atan2d("     ,  7,  0,  1,  1,  0, 19,  1,  1, -1,  2,  OP_ATAN2D   },
   { "dim("        ,  4,  0,  1,  1,  0, 19,  1,  1, -1,  2,  OP_DIM      },
   { "fmod("       ,  5,  0,  1,  1,  0, 19,  1,  1, -1,  2,  OP_MOD      },
   { "gauss("      ,  6,  0,  1,  1,  0, 19,  1,  1, -1,  2,  OP_GAUSS    },
   { "mod("        ,  4,  0,  1,  1,  0, 19,  1,  1, -1,  2,  OP_MOD      },
   { "pow("        ,  4,  0,  1,  1,  0, 19,  1,  1, -1,  2,  OP_POW      },
   { "rand("       ,  5,  0,  1,  1,  0, 19,  1,  1, -1,  2,  OP_RAND     },
   { "sign("       ,  5,  0,  1,  1,  0, 19,  1,  1, -1,  2,  OP_SIGN     },

/* Functions with two arguments. */
   { "qif("        ,  4,  0,  1,  1,  0, 19,  1,  1, -2,  3,  OP_QIF      },

/* Functions with variable numbers of arguments. */
   { "max("        ,  4,  0,  1,  1,  0, 19,  1,  1, -1, -2,  OP_MAX      },
   { "min("        ,  4,  0,  1,  1,  0, 19,  1,  1, -1, -2,  OP_MIN      },

/* Parenthesised expressions. */
   { ")"           ,  1,  1,  0,  0,  0,  2, 19, -1,  0,  0,  OP_NULL     },
   { "("           ,  1,  0,  1,  1,  0, 19,  1,  1,  0,  0,  OP_NULL     },

/* Unary arithmetic operators. */
   { "+"           ,  1,  0,  1,  1,  1, 17, 16,  0,  0,  0,  OP_NULL     },
   { "-"           ,  1,  0,  1,  1,  1, 17, 16,  0,  0,  0,  OP_NEG      },

/* Unary boolean operators. */
   { "!"           ,  1,  0,  1,  1,  0, 17, 16,  0,  0,  0,  OP_NOT      },
   { ".not."       ,  5,  0,  1,  1,  0, 17, 16,  0,  0,  0,  OP_NOT      },

/* Binary arithmetic operators. */
   { "**"          ,  2,  1,  1,  1,  0, 18, 15,  0, -1,  0,  OP_POW      },
   { "*"           ,  1,  1,  1,  1,  0, 14, 14,  0, -1,  0,  OP_MUL      },
   { "/"           ,  1,  1,  1,  1,  0, 14, 14,  0, -1,  0,  OP_DIV      },
   { "+"           ,  1,  1,  1,  1,  0, 13, 13,  0, -1,  0,  OP_ADD      },
   { "-"           ,  1,  1,  1,  1,  0, 13, 13,  0, -1,  0,  OP_SUB      },

/* Bit-shift operators. */
   { "<<"          ,  2,  1,  1,  1,  0, 12, 12,  0, -1,  0,  OP_SHFTL    },
   { ">>"          ,  2,  1,  1,  1,  0, 12, 12,  0, -1,  0,  OP_SHFTR    },

/* Relational operators. */
   { "<"           ,  1,  1,  1,  1,  0, 11, 11,  0, -1,  0,  OP_LT       },
   { ".lt."        ,  4,  1,  1,  1,  0, 11, 11,  0, -1,  0,  OP_LT       },
   { "<="          ,  2,  1,  1,  1,  0, 11, 11,  0, -1,  0,  OP_LE       },
   { ".le."        ,  4,  1,  1,  1,  0, 11, 11,  0, -1,  0,  OP_LE       },
   { ">"           ,  1,  1,  1,  1,  0, 11, 11,  0, -1,  0,  OP_GT       },
   { ".gt."        ,  4,  1,  1,  1,  0, 11, 11,  0, -1,  0,  OP_GT       },
   { ">="          ,  2,  1,  1,  1,  0, 11, 11,  0, -1,  0,  OP_GE       },
   { ".ge."        ,  4,  1,  1,  1,  0, 11, 11,  0, -1,  0,  OP_GE       },
   { "=="          ,  2,  1,  1,  1,  0, 10, 10,  0, -1,  0,  OP_EQ       },
   { ".eq."        ,  4,  1,  1,  1,  0, 10, 10,  0, -1,  0,  OP_EQ       },
   { "!="          ,  2,  1,  1,  1,  0, 10, 10,  0, -1,  0,  OP_NE       },
   { ".ne."        ,  4,  1,  1,  1,  0, 10, 10,  0, -1,  0,  OP_NE       },

/* Bit-wise operators. */
   { "&"           ,  1,  1,  1,  1,  0,  9,  9,  0, -1,  0,  OP_BITAND   },
   { "^"           ,  1,  1,  1,  1,  0,  8,  8,  0, -1,  0,  OP_BITXOR   },
   { "|"           ,  1,  1,  1,  1,  0,  7,  7,  0, -1,  0,  OP_BITOR    },

/* Binary boolean operators. */
   { "&&"          ,  2,  1,  1,  1,  0,  6,  6,  0, -1,  0,  OP_AND      },
   { ".and."       ,  5,  1,  1,  1,  0,  6,  6,  0, -1,  0,  OP_AND      },
   { "^^"          ,  2,  1,  1,  1,  0,  5,  5,  0, -1,  0,  OP_XOR      },
   { "||"          ,  2,  1,  1,  1,  0,  4,  4,  0, -1,  0,  OP_OR       },
   { ".or."        ,  4,  1,  1,  1,  0,  4,  4,  0, -1,  0,  OP_OR       },
   { ".eqv."       ,  5,  1,  1,  1,  0,  3,  3,  0, -1,  0,  OP_EQV      },
   { ".neqv."      ,  6,  1,  1,  1,  0,  3,  3,  0, -1,  0,  OP_XOR      },
   { ".xor."       ,  5,  1,  1,  1,  0,  3,  3,  0, -1,  0,  OP_XOR      },

/* Separators. */
   { ","           ,  1,  1,  1,  1,  0,  2,  2,  0,  0,  0,  OP_NULL     },

/* End of symbol data. */
   { NULL          ,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  OP_NULL     }
};

/* These variables identify indices in the above array which hold
   special symbols used explicitly in the code. */
static const int symbol_ldcon = 0; /* Load a constant */
static const int symbol_ldvar = 1; /* Load a variable */

/* 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(MathMap)

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



static pthread_mutex_t mutex2 = PTHREAD_MUTEX_INITIALIZER;
#define LOCK_MUTEX2 pthread_mutex_lock( &mutex2 );
#define UNLOCK_MUTEX2 pthread_mutex_unlock( &mutex2 );

static pthread_mutex_t mutex3 = PTHREAD_MUTEX_INITIALIZER;
#define LOCK_MUTEX3 pthread_mutex_lock( &mutex3 );
#define UNLOCK_MUTEX3 pthread_mutex_unlock( &mutex3 );

static pthread_mutex_t mutex4 = PTHREAD_MUTEX_INITIALIZER;
#define LOCK_MUTEX4 pthread_mutex_lock( &mutex4 );
#define UNLOCK_MUTEX4 pthread_mutex_unlock( &mutex4 );

static pthread_mutex_t mutex5 = PTHREAD_MUTEX_INITIALIZER;
#define LOCK_MUTEX5 pthread_mutex_lock( &mutex5 );
#define UNLOCK_MUTEX5 pthread_mutex_unlock( &mutex5 );

static pthread_mutex_t mutex6 = PTHREAD_MUTEX_INITIALIZER;
#define LOCK_MUTEX6 pthread_mutex_lock( &mutex6 );
#define UNLOCK_MUTEX6 pthread_mutex_unlock( &mutex6 );

static pthread_mutex_t mutex7 = PTHREAD_MUTEX_INITIALIZER;
#define LOCK_MUTEX7 pthread_mutex_lock( &mutex7 );
#define UNLOCK_MUTEX7 pthread_mutex_unlock( &mutex7 );

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

static char getattrib_buff[ 51 ];


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

#define LOCK_MUTEX2
#define UNLOCK_MUTEX2

#define LOCK_MUTEX3
#define UNLOCK_MUTEX3

#define LOCK_MUTEX4
#define UNLOCK_MUTEX4

#define LOCK_MUTEX5
#define UNLOCK_MUTEX5

#define LOCK_MUTEX6
#define UNLOCK_MUTEX6

#define LOCK_MUTEX7
#define UNLOCK_MUTEX7

#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. */
AstMathMap *astMathMapId_( int, int, int, const char *[], int, const char *[], const char *, ... );

/* Prototypes for Private Member Functions. */
/* ======================================== */
static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
static int GetObjSize( AstObject *, int * );
static const char *GetAttrib( AstObject *, const char *, int * );
static double Gauss( Rcontext *, int * );
static double LogGamma( double, int * );
static double Poisson( Rcontext *, double, int * );
static double Rand( Rcontext *, int * );
static int DefaultSeed( const Rcontext *, int * );
static int Equal( AstObject *, AstObject *, int * );
static int GetSeed( AstMathMap *, int * );
static int GetSimpFI( AstMathMap *, int * );
static int GetSimpIF( AstMathMap *, int * );
static int MapMerge( AstMapping *, int, int, int *, AstMapping ***, int **, int * );
static int TestAttrib( AstObject *, const char *, int * );
static int TestSeed( AstMathMap *, int * );
static int TestSimpFI( AstMathMap *, int * );
static int TestSimpIF( AstMathMap *, int * );
static void CleanFunctions( int, const char *[], char ***, int * );
static void ClearAttrib( AstObject *, const char *, int * );
static void ClearSeed( AstMathMap *, int * );
static void ClearSimpFI( AstMathMap *, int * );
static void ClearSimpIF( AstMathMap *, int * );
static void CompileExpression( const char *, const char *, const char *, int, const char *[], int **, double **, int *, int * );
static void CompileMapping( const char *, const char *, int, int, int, const char *[], int, const char *[], int ***, int ***, double ***, double ***, int *, int *, int * );
static void Copy( const AstObject *, AstObject *, int * );
static void Delete( AstObject *, int * );
static void Dump( AstObject *, AstChannel *, int * );
static void EvaluateFunction( Rcontext *, int, const double **, const int *, const double *, int, double *, int * );
static void EvaluationSort( const double [], int, int [], int **, int *, int * );
static void ExtractExpressions( const char *, const char *, int, const char *[], int, char ***, int * );
static void ExtractVariables( const char *, const char *, int, const char *[], int, int, int, int, int, char ***, int * );
static void ParseConstant( const char *, const char *, const char *, int, int *, double *, int * );
static void ParseName( const char *, int, int *, int * );
static void ParseVariable( const char *, const char *, const char *, int, int, const char *[], int *, int *, int * );
static void SetAttrib( AstObject *, const char *, int * );
static void SetSeed( AstMathMap *, int, int * );
static void SetSimpFI( AstMathMap *, int, int * );
static void SetSimpIF( AstMathMap *, int, int * );
static void ValidateSymbol( const char *, const char *, const char *, int, int, int *, int **, int **, int *, double **, int * );

/* Member functions. */
/* ================= */
static void CleanFunctions( int nfun, const char *fun[], char ***clean, int *status ) {
/*
*  Name:
*     CleanFunctions

*  Purpose:
*     Make a clean copy of a set of functions.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     void CleanFunctions( int nfun, const char *fun[], char ***clean, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     This function copies an array of strings, eliminating any white space
*     characters and converting to lower case. It is intended for cleaning
*     up arrays of function definitions prior to compilation. The returned
*     copy is stored in dynamically allocated memory.

*  Parameters:
*     nfun
*        The number of functions to be cleaned.
*     fun
*        Pointer to an array, with "nfun" elements, of pointers to null
*        terminated strings which contain each of the functions.
*     clean
*        Address in which to return a pointer to an array (with "nfun"
*        elements) of pointers to null terminated strings containing the
*        cleaned functions (i.e. this returns an array of strings).
*
*        Both the returned array of pointers, and the strings to which they
*        point, will be dynamically allocated and should be freed by the
*        caller (using astFree) when no longer required.
*     status
*        Pointer to the inherited status variable.

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

/* Local Variables: */
   char c;                       /* Character from function string */
   int i;                        /* Loop counter for characters */
   int ifun;                     /* Loop counter for functions */
   int nc;                       /* Count of non-blank characters */

/* Initialise. */
   *clean = NULL;

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

/* Allocate and initialise an array to hold the returned pointers. */
   MALLOC_POINTER_ARRAY( *clean, char *, nfun )

/* Loop through all the input functions. */
   if ( astOK ) {
      for ( ifun = 0; ifun < nfun; ifun++ ) {

/* Count the number of non-blank characters in each function string. */
         nc = 0;
         for ( i = 0; ( c = fun[ ifun ][ i ] ); i++ ) nc += !isspace( c );

/* Allocate a string long enough to hold the function with all the
   white space removed, storing its pointer in the array allocated
   earlier. Check for errors. */
         ( *clean )[ ifun ] = astMalloc( sizeof( char ) *
                                         (size_t) ( nc + 1 ) );
         if ( !astOK ) break;

/* Loop to copy the non-blank function characters into the new
   string. */
         nc = 0;
         for ( i = 0; ( c = fun[ ifun ][ i ] ); i++ ) {
            if ( !isspace( c ) ) ( *clean )[ ifun ][ nc++ ] = tolower( c );
         }

/* Null-terminate the result. */
         ( *clean )[ ifun ][ nc ] = '\0';
      }

/* If an error occurred, then free the main pointer array together
   with any strings that have been allocated, resetting the output
   value. */
      if ( !astOK ) {
         FREE_POINTER_ARRAY( *clean, nfun )
      }
   }
}

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

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

*  Type:
*     Private function.

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

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

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

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

/* Obtain a pointer to the MathMap structure. */
   this = (AstMathMap *) this_object;

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

/* Seed. */
/* ----- */
   if ( !strcmp( attrib, "seed" ) ) {
      astClearSeed( this );

/* SimpFI. */
/* ------- */
   } else if ( !strcmp( attrib, "simpfi" ) ) {
      astClearSimpFI( this );

/* SimpIF. */
/* ------- */
   } else if ( !strcmp( attrib, "simpif" ) ) {
      astClearSimpIF( this );

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

static void CompileExpression( const char *method, const char *class,
                               const char *exprs, int nvar, const char *var[],
                               int **code, double **con, int *stacksize, int *status ) {
/*
*  Name:
*     CompileExpression

*  Purpose:
*     Compile a mathematical expression.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     void CompileExpression( const char *method, const char *class,
*                             const char *exprs, int nvar, const char *var[],
*                             int **code, double **con, int *stacksize )

*  Class Membership:
*     MathMap member function.

*  Description:
*     This function checks and compiles a mathematical expression. It
*     produces a sequence of operation codes (opcodes) and a set of
*     numerical constants which may subsequently be used to evaluate the
*     expression on a push-down stack.

*  Parameters:
*     method
*        Pointer to a constant null-terminated character string
*        containing the name of the method that invoked this function.
*        This method name is used solely for constructing error messages.
*     class
*        Pointer to a constant null-terminated character string containing the
*        class name of the Object being processed. This name is used solely
*        for constructing error messages.
*     exprs
*        Pointer to a null-terminated string containing the expression
*        to be compiled. This is case sensitive and should contain no white
*        space.
*     nvar
*        The number of variable names defined for use in the expression.
*     var
*        An array of pointers (with "nvar" elements) to null-terminated
*        strings. Each of these should contain a variable name which may
*        appear in the expression. These strings are case sensitive and
*        should contain no white space.
*     code
*        Address of a pointer which will be set to point at a dynamically
*        allocated array of int containing the set of opcodes (cast to int)
*        produced by this function. The first element of this array will
*        contain a count of the number of opcodes which follow.
*
*        The allocated space must be freed by the caller (using astFree) when
*        no longer required.
*     con
*        Address of a pointer which will be set to point at a dynamically
*        allocated array of double containing the set of constants
*        produced by this function (this may be NULL if no constants are
*        produced).
*
*        The allocated space must be freed by the caller (using astFree) when
*        no longer required.
*     stacksize
*        Pointer to an int in which to return the size of the push-down stack
*        required to evaluate the expression using the returned opcodes and
*        constants.

*  Algorithm:
*     The function passes through the input expression searching for
*     symbols. It looks for standard symbols (arithmetic operators,
*     parentheses, function calls and delimiters) in the next part of the
*     expression to be parsed, using identification information stored in
*     the static "symbol" array. It ignores certain symbols, according to
*     whether they appear to be operators or operands. The choice depends on
*     what the previous symbol was; for instance, two operators may not
*     occur in succession. Unary +/- operators are also ignored in
*     situations where they are not permitted.
*
*     If a standard symbol is found, it is passed to the ValidateSymbol
*     function, which keeps track of the current level of parenthesis in the
*     expression and of the number of arguments supplied to any (possibly
*     nested) function calls. This function then accepts or rejects the
*     symbol according to whether it is valid within the current context. An
*     error is reported if it is rejected.
*
*     If the part of the expression currently being parsed did not contain a
*     standard symbol, an attempt is made to parse it first as a constant,
*     then as a variable name. If either of these succeeds, an appropriate
*     symbol number is added to the list of symbols identified so far, and a
*     value is added to the list of constants - this is either the value of
*     the constant itself, or the identification number of the variable. If
*     the expression cannot be parsed, an error is reported.
*
*     When the entire expression has been analysed as a sequence of symbols
*     (and associated constants), the EvaluationSort function is
*     invoked. This sorts the symbols into evaluation order, which is the
*     order in which the associated operations must be performed on a
*     push-down arithmetic stack to evaluate the expression. This routine
*     also substitutes operation codes (defined in the "Oper" enum) for the
*     symbol numbers and calculates the size of evaluation stack which will
*     be required.

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

/* Local Variables: */
   double c;                     /* Value of parsed constant */
   int *argcount;                /* Array of argument count information */
   int *opensym;                 /* Array of opening parenthesis information */
   int *symlist;                 /* Array of symbol indices */
   int found;                    /* Standard symbol identified? */
   int iend;                     /* Ending index in the expression string */
   int istart;                   /* Staring index in the expression string */
   int isym;                     /* Loop counter for symbols */
   int ivar;                     /* Index of variable name */
   int lpar;                     /* Parenthesis level */
   int ncon;                     /* Number of constants generated */
   int nsym;                     /* Number of symbols identified */
   int opernext;                 /* Next symbol an operator (from left)? */
   int size;                     /* Size of symbol matched */
   int sym;                      /* Index of symbol in static "symbol" array */
   int unarynext;                /* Next symbol may be unary +/- ? */

/* Initialise. */
   *code = NULL;
   *con = NULL;
   *stacksize = 0;

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

/* Further initialisation. */
   argcount = NULL;
   lpar = 0;
   ncon = 0;
   nsym = 0;
   opensym = NULL;
   symlist = NULL;
   sym = 0;
   ivar = 0;

/* The first symbol to be encountered must not look like an operator
   from the left. It may be a unary + or - operator. */
   opernext = 0;
   unarynext = 1;

/* Search through the expression to classify each symbol which appears
   in it. Stop when there are no more input characters or an error is
   detected. */
   istart = 0;
   for ( istart = 0; astOK && exprs[ istart ]; istart = iend + 1 ) {

/* Compare each of the symbols in the symbol data with the next
   section of the expression, looking for the longest symbol text which
   will match. Stop if a NULL "text" value is found, which acts as the
   end flag. */
      found = 0;
      size = 0;
      for ( isym = 0; symbol[ isym ].text; isym++ ) {

/* Only consider symbols which have text associated with them and
   which look like operators or operands from the left, according to the
   setting of the "opernext" flag. Thus, if an operator or operand is
   missing from the input expression, the next symbol will not be
   identified, because it will be of the wrong type. Also exclude unary
   +/- operators if they are out of context. */
         if ( symbol[ isym ].size &&
              ( symbol[ isym ].operleft == opernext ) &&
              ( !symbol[ isym ].unaryoper || unarynext ) ) {

/* Test if the text of the symbol matches the expression at the
   current position. If so, note that a match has been found. */
            if ( !strncmp( exprs + istart, symbol[ isym ].text,
                           (size_t) symbol[ isym ].size ) ) {
               found = 1;

/* If this symbol matches more characters than any previous symbol,
   then store the symbol's index and note its size. */
               if ( symbol[ isym ].size > size ) {
                  sym = isym;
                  size = symbol[ isym ].size;

/* Calculate the index of the last symbol character in the expression
   string. */
                  iend = istart + size - 1;
               }
            }
         }
      }

/* If the symbol was identified as one of the standard symbols, then
   validate it, updating the parenthesis level and argument count
   information at the same time. */
      if ( found ) {
         ValidateSymbol( method, class, exprs, iend, sym, &lpar, &argcount,
                         &opensym, &ncon, con, status );

/* If it was not one of the standard symbols, then check if the next
   symbol was expected to be an operator. If so, then there is a missing
   operator, so report an error. */
      } else {
         if ( opernext ) {
            astError( AST__MIOPR,
                      "%s(%s): Missing or invalid operator in the expression "
                      "\"%.*s\".", status,
                      method, class, istart + 1, exprs );

/* If the next symbol was expected to be an operand, then it may be a
   constant, so try to parse it as one. */
         } else {
            ParseConstant( method, class, exprs, istart, &iend, &c, status );
            if ( astOK ) {

/* If successful, set the symbol number to "symbol_ldcon" (load
   constant) and extend the "*con" array to accommodate a new
   constant. Check for errors. */
               if ( iend >= istart ) {
                  sym = symbol_ldcon;
                  *con = astGrow( *con, ncon + 1, sizeof( double ) );
                  if ( astOK ) {

/* Append the constant to the "*con" array. */
                     ( *con )[ ncon++ ] = c;
                  }

/* If the symbol did not parse as a constant, then it may be a
   variable name, so try to parse it as one. */
               } else {
                  ParseVariable( method, class, exprs, istart, nvar, var,
                                 &ivar, &iend, status );
                  if ( astOK ) {

/* If successful, set the symbol to "symbol_ldvar" (load variable) and
   extend the "*con" array to accommodate a new constant. Check for
   errors. */
                     if ( ivar != -1 ) {
                        sym = symbol_ldvar;
                        *con = astGrow( *con, ncon + 1, sizeof( double ) );
                        if ( astOK ) {

/* Append the variable identification number as a constant to the
   "*con" array. */
                           ( *con )[ ncon++ ] = (double) ivar;
                        }

/* If the expression did not parse as a variable name, then there is a
   missing operand in the expression, so report an error. */
                     } else {
                        astError( AST__MIOPA,
                                  "%s(%s): Missing or invalid operand in the "
                                  "expression \"%.*s\".", status,
                                  method, class, istart + 1, exprs );
                     }
                  }
               }
            }
         }
      }

/* If there has been no error, then the next symbol in the input
   expression has been identified and is valid. */
      if ( astOK ) {

/* Decide whether the next symbol should look like an operator or an
   operand from the left. This is determined by the nature of the symbol
   just identified (seen from the right) - two operands or two operators
   cannot be adjacent. */
         opernext = !symbol[ sym ].operright;

/* Also decide whether the next symbol may be a unary +/- operator,
   according to the "unarynext" symbol data entry for the symbol just
   identified. */
         unarynext = symbol[ sym ].unarynext;

/* Extend the "symlist" array to accommodate the symbol just
   identified. Check for errors. */
         symlist = astGrow( symlist, nsym + 1, sizeof( int ) );
         if ( astOK ) {

/* Append the symbol's index to the end of this list. */
            symlist[ nsym++ ] = sym;
         }
      }
   }

/* If there has been no error, check the final context after
   identifying all the symbols... */
   if ( astOK ) {

/* If an operand is still expected, then there is an unsatisfied
   operator on the end of the expression, so report an error. */
      if ( !opernext ) {
         astError( AST__MIOPA,
                   "%s(%s): Missing or invalid operand in the expression "
                   "\"%s\".", status,
                   method, class, exprs );

/* If the final parenthesis level is positive, then there is a missing
   right parenthesis, so report an error. */
      } else if ( lpar > 0 ) {
         astError( AST__MRPAR,
                   "%s(%s): Missing right parenthesis in the expression "
                   "\"%s\".", status,
                   method, class, exprs );
      }
   }

/* Sort the symbols into evaluation order to produce output opcodes. */
   EvaluationSort( *con, nsym, symlist, code, stacksize, status );

/* Free any memory used as workspace. */
   if ( argcount ) argcount = astFree( argcount );
   if ( opensym ) opensym = astFree( opensym );
   if ( symlist ) symlist = astFree( symlist );

/* If OK, re-allocate the "*con" array to have the correct size (since
   astGrow may have over-allocated space). */
   if ( astOK && *con ) {
      *con = astRealloc( *con, sizeof( double ) * (size_t) ncon );
   }

/* If an error occurred, free any allocated memory and reset the
   output values. */
   if ( !astOK ) {
      *code = astFree( *code );
      *con = astFree( *con );
      *stacksize = 0;
   }
}

static void CompileMapping( const char *method, const char *class,
                            int nin, int nout,
                            int nfwd, const char *fwdfun[],
                            int ninv, const char *invfun[],
                            int ***fwdcode, int ***invcode,
                            double ***fwdcon, double ***invcon,
                            int *fwdstack, int *invstack, int *status ) {
/*
*  Name:
*     CompileMapping

*  Purpose:
*     Compile the transformation functions for a MathMap.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     void CompileMapping( const char *method, const char *class,
*                          int nin, int nout,
*                          int nfwd, const char *fwdfun[],
*                          int ninv, const char *invfun[],
*                          int ***fwdcode, int ***invcode,
*                          double ***fwdcon, double ***invcon,
*                          int *fwdstack, int *invstack, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     This function checks and compiles the transformation functions required
*     to create a MathMap. It produces sequences of operation codes (opcodes)
*     and numerical constants which may subsequently be used to evaluate the
*     functions on a push-down stack.

*  Parameters:
*     method
*        Pointer to a constant null-terminated character string
*        containing the name of the method that invoked this function.
*        This method name is used solely for constructing error messages.
*     class
*        Pointer to a constant null-terminated character string containing the
*        class name of the Object being processed. This name is used solely
*        for constructing error messages.
*     nin
*        Number of input variables for the MathMap.
*     nout
*        Number of output variables for the MathMap.
*     nfwd
*        The number of forward transformation functions being supplied.
*        This must be at least equal to "nout".
*     fwdfun
*        Pointer to an array, with "nfwd" elements, of pointers to null
*        terminated strings which contain each of the forward transformation
*        functions. These must be in lower case and should contain no white
*        space.
*     ninv
*        The number of inverse transformation functions being supplied.
*        This must be at least equal to "nin".
*     invfun
*        Pointer to an array, with "ninv" elements, of pointers to null
*        terminated strings which contain each of the inverse transformation
*        functions. These must be in lower case and should contain no white
*        space.
*     fwdcode
*        Address in which to return a pointer to an array (with "nfwd"
*        elements) of pointers to arrays of int containing the set of opcodes
*        (cast to int) for each forward transformation function. The number
*        of opcodes produced for each function is given by the first element
*        of the opcode array.
*
*        Both the returned array of pointers, and the arrays of int to which
*        they point, will be stored in dynamically allocated memory and should
*        be freed by the caller (using astFree) when no longer required.
*
*        If the right hand sides (including the "=" sign) of all the supplied
*        functions are absent, then this indicates an undefined transformation
*        and the returned pointer value will be NULL. An error results if
*        an "=" sign is present but no expression follows it.
*     invcode
*        Address in which to return a pointer to an array (with "ninv"
*        elements) of pointers to arrays of int containing the set of opcodes
*        (cast to int) for each inverse transformation function. The number
*        of opcodes produced for each function is given by the first element
*        of the opcode array.
*
*        Both the returned array of pointers, and the arrays of int to which
*        they point, will be stored in dynamically allocated memory and should
*        be freed by the caller (using astFree) when no longer required.
*
*        If the right hand sides (including the "=" sign) of all the supplied
*        functions are absent, then this indicates an undefined transformation
*        and the returned pointer value will be NULL. An error results if
*        an "=" sign is present but no expression follows it.
*     fwdcon
*        Address in which to return a pointer to an array (with "nfwd"
*        elements) of pointers to arrays of double containing the set of
*        constants for each forward transformation function.
*
*        Both the returned array of pointers, and the arrays of double to which
*        they point, will be stored in dynamically allocated memory and should
*        be freed by the caller (using astFree) when no longer required. Note
*        that any of the pointers to the arrays of double may be NULL if no
*        constants are associated with a particular function.
*
*        If the forward transformation is undefined, then the returned pointer
*        value will be NULL.
*     invcon
*        Address in which to return a pointer to an array (with "ninv"
*        elements) of pointers to arrays of double containing the set of
*        constants for each inverse transformation function.
*
*        Both the returned array of pointers, and the arrays of double to which
*        they point, will be stored in dynamically allocated memory and should
*        be freed by the caller (using astFree) when no longer required. Note
*        that any of the pointers to the arrays of double may be NULL if no
*        constants are associated with a particular function.
*
*        If the inverse transformation is undefined, then the returned pointer
*        value will be NULL.
*     fwdstack
*        Pointer to an int in which to return the size of the push-down stack
*        required to evaluate the forward transformation functions.
*     invstack
*        Pointer to an int in which to return the size of the push-down stack
*        required to evaluate the inverse transformation functions.
*     status
*        Pointer to the inherited status variable.

*  Notes:
*     - A value of NULL will be returned for the "*fwdcode", "*invcode",
*     "*fwdcon" and "*invcon" pointers and a value of zero will be returned
*     for the "*fwdstack" and "*invstack" values if this function is invoked
*     with the global error status set, or if it should fail for any reason.
*/

/* Local Variables: */
   char **exprs;                 /* Pointer to array of expressions */
   char **var;                   /* Pointer to array of variable names */
   const char **strings;         /* Pointer to temporary array of strings */
   int ifun;                     /* Loop counter for functions */
   int nvar;                     /* Number of variables to extract */
   int stacksize;                /* Required stack size */

/* Initialise. */
   *fwdcode = NULL;
   *invcode = NULL;
   *fwdcon = NULL;
   *invcon = NULL;
   *fwdstack = 0;
   *invstack = 0;
   nvar = 0;

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

/* Further initialisation. */
   exprs = NULL;
   var = NULL;

/* Compile the forward transformation. */
/* ----------------------------------- */
/* Allocate space for an array of pointers to the functions from which
   we will extract variable names. */
   strings = astMalloc( sizeof( char * ) * (size_t) ( nin + nfwd ) );

/* Fill the first elements of this array with pointers to the inverse
   transformation functions ("nin" in number) which yield the final input
   values. These will have the names of the input variables on their left
   hand sides. */
   if ( astOK ) {
      nvar = 0;
      for ( ifun = ninv - nin; ifun < ninv; ifun++ ) {
         strings[ nvar++ ] = invfun[ ifun ];
      }

/* Fill the remaining elements of the array with pointers to the
   forward transformation functions. These will have the names of any
   intermediate variables plus the final output variables on their left
   hand sides. */
      for ( ifun = 0; ifun < nfwd; ifun++ ) strings[ nvar++ ] = fwdfun[ ifun ];

/* Extract the variable names from the left hand sides of these
   functions and check them for validity and absence of duplication. */
      ExtractVariables( method, class, nvar, strings, nin, nout, nfwd, ninv, 1,
                        &var, status );
   }

/* Free the temporary array of string pointers. */
   strings = astFree( strings );

/* Extract the expressions from the right hand sides of the forward
   transformation functions. */
   ExtractExpressions( method, class, nfwd, fwdfun, 1, &exprs, status );

/* If OK, and the forward transformation is defined, then allocate and
   initialise space for an array of pointers to the opcodes for each
   expression and, similarly, for the constants for each expression. */
   if ( astOK && exprs ) {
      MALLOC_POINTER_ARRAY( *fwdcode, int *, nfwd )
      MALLOC_POINTER_ARRAY( *fwdcon, double *, nfwd )

/* If OK, loop to compile each of the expressions, storing pointers to
   the resulting opcodes and constants in the arrays allocated above. On
   each loop, we make progressively more of the variable names in "var"
   visible to the compilation function. This ensures that each expression
   can only use variables which have been defined earlier. */
      if ( astOK ) {
         for ( ifun = 0; ifun < nfwd; ifun++ ) {
            CompileExpression( method, class, exprs[ ifun ],
                               nin + ifun, (const char **) var,
                               &( *fwdcode )[ ifun ], &( *fwdcon )[ ifun ],
                               &stacksize, status );

/* If an error occurs, then report contextual information and quit. */
            if ( !astOK ) {
               astError( astStatus,
                         "Error in forward transformation function %d.", status,
                         ifun + 1 );
               break;
            }

/* If OK, calculate the maximum evaluation stack size required by any
   of the expressions. */
            *fwdstack = ( *fwdstack > stacksize ) ? *fwdstack : stacksize;
         }
      }
   }

/* Free the memory containing the extracted expressions and variables. */
   FREE_POINTER_ARRAY( exprs, nfwd )
   FREE_POINTER_ARRAY( var, nvar )

/* Compile the inverse transformation. */
/* ----------------------------------- */
/* Allocate space for an array of pointers to the functions from which
   we will extract variable names. */
   strings = astMalloc( sizeof( char * ) * (size_t) ( nout + ninv ) );

/* Fill the first elements of this array with pointers to the forward
   transformation functions ("nout" in number) which yield the final
   output values. These will have the names of the output variables on
   their left hand sides. */
   if ( astOK ) {
      nvar = 0;
      for ( ifun = nfwd - nout; ifun < nfwd; ifun++ ) {
         strings[ nvar++ ] = fwdfun[ ifun ];
      }

/* Fill the remaining elements of the array with pointers to the
   inverse transformation functions. These will have the names of any
   intermediate variables plus the final input variables on their left
   hand sides. */
      for ( ifun = 0; ifun < ninv; ifun++ ) strings[ nvar++ ] = invfun[ ifun ];

/* Extract the variable names from the left hand sides of these
   functions and check them for validity and absence of duplication. */
      ExtractVariables( method, class, nvar, strings, nin, nout, nfwd, ninv, 0,
                        &var, status );
   }

/* Free the temporary array of string pointers. */
   strings = astFree( strings );

/* Extract the expressions from the right hand sides of the inverse
   transformation functions. */
   ExtractExpressions( method, class, ninv, invfun, 0, &exprs, status );

/* If OK, and the forward transformation is defined, then allocate and
   initialise space for an array of pointers to the opcodes for each
   expression and, similarly, for the constants for each expression. */
   if ( astOK && exprs ) {
      MALLOC_POINTER_ARRAY( *invcode, int *, ninv )
      MALLOC_POINTER_ARRAY( *invcon, double *, ninv )

/* If OK, loop to compile each of the expressions, storing pointers to
   the resulting opcodes and constants in the arrays allocated above. On
   each loop, we make progressively more of the variable names in "var"
   visible to the compilation function. This ensures that each expression
   can only use variables which have been defined earlier. */
      if ( astOK ) {
         for ( ifun = 0; ifun < ninv; ifun++ ) {
            CompileExpression( method, class, exprs[ ifun ],
                               nout + ifun, (const char **) var,
                               &( *invcode )[ ifun ], &( *invcon )[ ifun ],
                               &stacksize, status );

/* If an error occurs, then report contextual information and quit. */
            if ( !astOK ) {
               astError( astStatus,
                         "Error in inverse transformation function %d.", status,
                         ifun + 1 );
               break;
            }

/* If OK, calculate the maximum evaluation stack size required by any
   of the expressions. */
            *invstack = ( *invstack > stacksize ) ? *invstack : stacksize;
         }
      }
   }

/* Free the memory containing the extracted expressions and variables. */
   FREE_POINTER_ARRAY( exprs, ninv )
   FREE_POINTER_ARRAY( var, nvar )

/* If an error occurred, then free all remaining allocated memory and
   reset the output values. */
   if ( !astOK ) {
      FREE_POINTER_ARRAY( *fwdcode, nfwd )
      FREE_POINTER_ARRAY( *invcode, ninv )
      FREE_POINTER_ARRAY( *fwdcon, nfwd )
      FREE_POINTER_ARRAY( *invcon, ninv )
      *fwdstack = 0;
      *invstack = 0;
   }
}

static int DefaultSeed( const Rcontext *context, int *status ) {
/*
*  Name:
*     DefaultSeed

*  Purpose:
*     Generate an unpredictable seed for a random number generator.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     int DefaultSeed( Rcontext *context, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     On each invocation this function returns an integer value which is
*     highly unpredictable. This value may be used as a default seed for the
*     random number generator associated with a MathMap, so that it
*     generates a different sequence on each occasion.

*  Parameters:
*     context
*        Pointer to the random number generator context associated with
*        the MathMap.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     The unpredictable integer.

*  Notes:
*     - This function does not perform error checking and will execute even
*     if the global error status is set.
*/

/* Local Constants: */
   const int nwarm = 5;          /* Number of warm-up iterations */
   const long int a = 8121L;     /* Constants for random number generator... */
   const long int c = 28411L;
   const long int m = 134456L;

/* Local Variables; */
   int iwarm;                    /* Loop counter for warm-up iterations */
   static long init = 0;         /* Local initialisation performed? */
   static long int rand;         /* Local random integer */
   unsigned long int bits;       /* Bit pattern for producing result */

/* On the first invocation, initialise a local random number generator
   to a value derived by combining bit patterns obtained from the system
   clock and the processor time used. The result needs to be positive and
   lie in the range 0 to "m-1". */
   LOCK_MUTEX5
   if ( !init ) {
      rand = (long int) ( ( (unsigned long int) time( NULL ) ^
                            (unsigned long int) clock() ) %
                          (unsigned long int) m );

/* These values will typically only change in their least significant
   bits between programs run successively, but by using the bit pattern
   as a seed, we ensure that these differences are rapidly propagated to
   other bits. To hasten this process, we "warm up" the local generator
   with a few iterations. This is a quick and dirty generator using
   constants from Press et al. (Numerical recipes). */
      for ( iwarm = 0; iwarm < nwarm; iwarm++ ) {
         rand = ( rand * a + c ) % m;
      }

/* Note that this initialisation has been performed. */
      init = 1;
   }
   UNLOCK_MUTEX5

/* Generate a new bit pattern from the system time. Apart from the
   first invocation, this will be a different time to that used above. */
   bits = (unsigned long int) time( NULL );

/* Mask in a pattern derived from the CPU time used. */
   bits ^= (unsigned long int) clock();

/* The system time may change quite slowly (e.g. every second), so
   also mask in the address of the random number generator context
   supplied. This makes the seed depend on which MathMap is in use. */
   bits ^= (unsigned long int) context;

/* Now mask in the last random integer produced by the random number
   generator whose context has been supplied. This makes the seed depend
   on the MathMap's past use of random numbers. */
   bits ^= (unsigned long int) context->random_int;

/* Finally, in order to produce different seeds when this function is
   invoked twice in rapid succession on the same object (with no
   intermediate processing), we also mask in a pseudo-random value
   generated here. Generate the next local random integer. */
   rand = ( rand * a + c ) % m;

/* We then scale this value to give an integer in the range 0 to
   ULONG_MAX and mask the corresponding bit pattern into our seed. */
   bits ^= (unsigned long int) ( ( (double) rand / (double) ( m - 1UL ) ) *
                                 ( ( (double) ULONG_MAX + 1.0 ) *
                                   ( 1.0 - DBL_EPSILON ) ) );

/* Return the integer value of the seed (which may involve discarding
   some unwanted bits). */
   return (int) bits;
}

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

*  Purpose:
*     Test if two MathMaps are equivalent.

*  Type:
*     Private function.

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

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

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

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

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

*  Notes:
*     - The two MathMaps are considered equivalent if the combination of
*     the first in series with the inverse of the second simplifies to a
*     UnitMap.
*     - 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: */
   AstMathMap *that;          /* Pointer to the second MathMap structure */
   AstMathMap *this;          /* Pointer to the first MathMap structure */
   double **that_con;         /* Lists of constants from "that" */
   double **this_con;         /* Lists of constants from "this" */
   int **that_code;           /* Lists of opcodes from "that" */
   int **this_code;           /* Lists of opcodes from "this" */
   int code;                  /* Opcode value */
   int icode;                 /* Opcode index */
   int icon;                  /* Constant index */
   int ifun;                  /* Function index */
   int ncode;                 /* No. of opcodes for current "this" function */
   int ncode_that;            /* No. of opcodes for current "that" function */
   int nin;                   /* Number of inputs */
   int nout;                  /* Number of outputs */
   int pass;                  /* Check fwd or inv */
   int result;                /* Result value to return */
   int that_nfun;             /* Number of functions from "that" */
   int this_nfun;             /* Number of functions from "this" */

/* Initialise. */
   result = 0;

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

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

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

/* Check they have the same number of inputs and outputs */
      nin = astGetNin( this );
      nout = astGetNout( this );
      if( astGetNout( that ) == nout && astGetNin( that ) == nin ) {

/* Assume equality. */
         result = 1;

/* The first pass through this next loop compares forward functions, and
   the second pass compares inverse functions. */
         for( pass = 0; pass < 2 && result; pass++ ) {

/* On the first pass, get pointers to the lists of opcodes and constants for
   the effective forward transformations (taking into account the value
   of the Invert attribute), together with the number of such functions. */
            if( pass == 0 ) {
               if( !astGetInvert( this ) ) {
                  this_code = this->fwdcode;
                  this_con = this->fwdcon;
                  this_nfun = this->nfwd;
               } else {
                  this_code = this->invcode;
                  this_con = this->invcon;
                  this_nfun = this->ninv;
               }

               if( !astGetInvert( that ) ) {
                  that_code = that->fwdcode;
                  that_con = that->fwdcon;
                  that_nfun = that->nfwd;
               } else {
                  that_code = that->invcode;
                  that_con = that->invcon;
                  that_nfun = that->ninv;
               }

/* On the second pass, get pointers to the lists of opcodes and constants for
   the effective inverse transformations, together with the number of such
   functions. */
            } else {

               if( astGetInvert( this ) ) {
                  this_code = this->fwdcode;
                  this_con = this->fwdcon;
                  this_nfun = this->nfwd;
               } else {
                  this_code = this->invcode;
                  this_con = this->invcon;
                  this_nfun = this->ninv;
               }

               if( astGetInvert( that ) ) {
                  that_code = that->fwdcode;
                  that_con = that->fwdcon;
                  that_nfun = that->nfwd;
               } else {
                  that_code = that->invcode;
                  that_con = that->invcon;
                  that_nfun = that->ninv;
               }
            }

/* Check that "this" and "that" have the same number of functions */
            if( that_nfun != this_nfun ) result = 0;

/* Loop round each function. */
            for( ifun = 0; ifun < this_nfun && result; ifun++ ) {

/* The first element in the opcode array is the number of subsequent
   opcodes. Obtain and compare these counts. */
               ncode = this_code ? this_code[ ifun ][ 0 ] : 0;
               ncode_that = that_code ? that_code[ ifun ][ 0 ] : 0;
               if( ncode != ncode_that ) result = 0;

/* Compare the following opcodes. Some opcodes consume constants from the
   list of constants associated with the MathMap. Compare the constants
   for such opcodes. */
               icon = 0;
               for( icode = 0; icode < ncode && result; icode++ ){
                  code = this_code[ ifun ][ icode ];
                  if( that_code[ ifun ][ icode ] != code ) {
                     result = 0;

                  } else if( code == OP_LDCON ||
                             code == OP_LDVAR ||
                             code == OP_MAX ||
                             code == OP_MIN ) {

                     if( this_con[ ifun ][ icon ] !=
                         that_con[ ifun ][ icon ] ) {
                        result = 0;
                     } else {
                        icon++;
                     }
                  }
               }
            }
         }
      }
   }

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

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

static void EvaluateFunction( Rcontext *rcontext, int npoint,
                              const double **ptr_in, const int *code,
                              const double *con, int stacksize, double *out, int *status ) {
/*
*  Name:
*     EvaluateFunction

*  Purpose:
*     Evaluate a compiled function.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     void EvaluateFunction( Rcontext *rcontext, int npoint,
*                            const double **ptr_in, const int *code,
*                            const double *con, int stacksize, double *out, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     This function implements a "virtual machine" which executes operations
*     on an arithmetic stack in order to evaluate transformation functions.
*     Each operation is specified by an input operation code (opcode) and
*     results in the execution of a vector operation on a stack. The final
*     result, after executing all the supplied opcodes, is returned as a
*     vector.
*
*     This function detects arithmetic errors (such as overflow and division
*     by zero) and propagates any "bad" coordinate values, including those
*     present in the input, to the output.

*  Parameters:
*     npoint
*        The number of points to be transformd (i.e. the size of the vector
*        of values on which operations are to be performed).
*     ptr_in
*        Pointer to an array of pointers to arrays of double (with "npoint"
*        elements). These arrays should contain the input coordinate values,
*        such that coordinate number "coord" for point number "point" can be
*        found in "ptr_in[coord][point]".
*     code
*        Pointer to an array of int containing the set of opcodes (cast to int)
*        for the operations to be performed. The first element of this array
*        should contain a count of the number of opcodes which follow.
*     con
*        Pointer to an array of double containing the set of constants required
*        to evaluate the function (this may be NULL if no constants are
*        required).
*     stacksize
*        The size of the stack required to evaluate the expression using the
*        opcodes and constants supplied. This value should be calculated during
*        expression compilation.
*     out
*        Pointer to an array of double (with "npoint" elements) in which to
*        return the vector of result values.
*     status
*        Pointer to the inherited status variable.
*/

/* Local Constants: */
   const int bits =              /* Number of bits in an unsigned long */
      sizeof( unsigned long ) * CHAR_BIT;
   const double eps =            /* Smallest number subtractable from 2.0 */
      2.0 * DBL_EPSILON;
   const double scale =          /* 2.0 raised to the power "bits" */
      ldexp( 1.0, bits );
   const double scale1 =         /* 2.0 raised to the power "bits-1" */
      scale * 0.5;
   const double rscale =         /* Reciprocal scale factor */
      1.0 / scale;
   const double rscale1 =        /* Reciprocal initial scale factor */
      1.0 / scale1;
   const int nblock =            /* Number of blocks of bits to process */
      ( sizeof( double ) + sizeof( unsigned long ) - 1 ) /
      sizeof( unsigned long );
   const unsigned long signbit = /* Mask for extracting sign bit */
      1UL << ( bits - 1 );

/* Local Variables: */
   double **stack;               /* Array of pointers to stack elements */
   double *work;                 /* Pointer to stack workspace */
   double *xv1;                  /* Pointer to first argument vector */
   double *xv2;                  /* Pointer to second argument vector */
   double *xv3;                  /* Pointer to third argument vector */
   double *xv;                   /* Pointer to sole argument vector */
   double *y;                    /* Pointer to result */
   double *yv;                   /* Pointer to result vector */
   double abs1;                  /* Absolute value (temporary variable) */
   double abs2;                  /* Absolute value (temporary variable) */
   double frac1;                 /* First (maybe normalised) fraction */
   double frac2;                 /* Second (maybe normalised) fraction */
   double frac;                  /* Sole normalised fraction */
   double newexp;                /* New power of 2 exponent value */
   double ran;                   /* Random number */
   double result;                /* Function result value */
   double unscale;               /* Factor for removing scaling */
   double value;                 /* Value to be assigned to stack vector */
   double x1;                    /* First argument value */
   double x2;                    /* Second argument value */
   double x3;                    /* Third argument value */
   double x;                     /* Sole argument value */
   int expon1;                   /* First power of 2 exponent */
   int expon2;                   /* Second power of 2 exponent */
   int expon;                    /* Sole power of 2 exponent */
   int iarg;                     /* Loop counter for arguments */
   int iblock;                   /* Loop counter for blocks of bits */
   int icode;                    /* Opcode value */
   int icon;                     /* Counter for number of constants used */
   int istk;                     /* Loop counter for stack elements */
   int ivar;                     /* Input variable number */
   int narg;                     /* Number of function arguments */
   int ncode;                    /* Number of opcodes to process */
   int point;                    /* Loop counter for stack vector elements */
   int sign;                     /* Argument is non-negative? */
   int tos;                      /* Top of stack index */
   static double d2r;            /* Degrees to radians conversion factor */
   static double log2;           /* Natural logarithm of 2.0 */
   static double pi;             /* Value of PI */
   static double r2d;            /* Radians to degrees conversion factor */
   static double rsafe_sq;       /* Reciprocal of "safe_sq" */
   static double safe_sq;        /* Huge value that can safely be squared */
   static int init = 0;          /* Initialisation performed? */
   unsigned long b1;             /* Block of bits from first argument */
   unsigned long b2;             /* Block of bits from second argument */
   unsigned long b;              /* Block of bits for result */
   unsigned long neg;            /* Result is negative? (sign bit) */

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

/* If this is the first invocation of this function, then initialise
   constant values. */
   LOCK_MUTEX2
   if ( !init ) {

/* Trigonometrical conversion factors. */
      pi = acos( -1.0 );
      r2d = 180.0 / pi;
      d2r = pi / 180.0;

/* Natural logarithm of 2.0. */
      log2 = log( 2.0 );

/* This value must be safe to square without producing overflow, yet
   large enough that adding or subtracting 1.0 from the square makes no
   difference. We also need its reciprocal. */
      safe_sq = 0.9 * sqrt( DBL_MAX );
      rsafe_sq = 1.0 / safe_sq;

/* Note that initialisation has been performed. */
      init = 1;
   }
   UNLOCK_MUTEX2

/* Allocate space for an array of pointers to elements of the
   workspace stack (each stack element being an array of double). */
   stack = astMalloc( sizeof( double * ) * (size_t) stacksize );

/* Allocate space for the stack itself. */
   work = astMalloc( sizeof( double ) *
                     (size_t) ( npoint * ( stacksize - 1 ) ) );

/* If OK, then initialise the stack pointer array to identify the
   start of each vector on the stack. The first element points at the
   output array (in which the result will be accumulated), while other
   elements point at successive vectors within the workspace allocated
   above. */
   if ( astOK ) {
      stack[ 0 ] = out;
      for ( istk = 1; istk < stacksize; istk++ ) {
         stack[ istk ] = work + ( istk - 1 ) * npoint;
      }

/* Define stack operations. */
/* ======================== */
/* We now define a set of macros for performing vector operations on
   elements of the stack. Each is in the form of a "case" block for
   execution in response to the appropriate operation code (opcode). */

/* Zero-argument operation. */
/* ------------------------ */
/* This macro performs a zero-argument operation, which results in the
   insertion of a new vector on to the stack. */
#define ARG_0(oper,setup,function) \
\
/* Test for the required opcode value. */ \
   case oper: \
\
/* Perform any required initialisation. */ \
      {setup;} \
\
/* Increment the top of stack index and obtain a pointer to the new stack \
   element (vector). */ \
      yv = stack[ ++tos ]; \
\
/* Loop to access each vector element, obtaining a pointer to it. */ \
      for ( point = 0; point < npoint; point++ ) { \
         y = yv + point; \
\
/* Perform the processing, which results in assignment to this element. */ \
         {function;} \
      } \
\
/* Break out of the "case" block. */ \
      break;

/* One-argument operation. */
/* ----------------------- */
/* This macro performs a one-argument operation, which processes the
   top stack element without changing the stack size. */
#define ARG_1(oper,function) \
\
/* Test for the required opcode value. */ \
   case oper: \
\
/* Obtain a pointer to the top stack element (vector). */ \
      xv = stack[ tos ]; \
\
/* Loop to access each vector element, obtaining its value and \
   checking that it is not bad. */ \
      for ( point = 0; point < npoint; point++ ) { \
         if ( ( x = xv[ point ] ) != AST__BAD ) { \
\
/* Also obtain a pointer to the element. */ \
            y = xv + point; \
\
/* Perform the processing, which uses the element's value and then \
   assigns the result to this element. */ \
            {function;} \
         } \
      } \
\
/* Break out of the "case" block. */ \
      break;

/* One-argument boolean operation. */
/* ------------------------------- */
/* This macro is similar in function to ARG_1 above, except that no
   checks are made for bad argument values. It is intended for use with
   boolean functions where bad values are handled explicitly. */
#define ARG_1B(oper,function) \
\
/* Test for the required opcode value. */ \
   case oper: \
\
/* Obtain a pointer to the top stack element (vector). */ \
      xv = stack[ tos ]; \
\
/* Loop to access each vector element, obtaining the argument value \
   and a pointer to the element. */ \
      for ( point = 0; point < npoint; point++ ) { \
         x = xv[ point ]; \
         y = xv + point; \
\
/* Perform the processing, which uses the element's value and then \
   assigns the result to this element. */ \
         {function;} \
      } \
\
/* Break out of the "case" block. */ \
      break;

/* Two-argument operation. */
/* ----------------------- */
/* This macro performs a two-argument operation, which processes the
   top two stack elements and produces a single result, resulting in the
   stack size decreasing by one. In this case, we first define a macro
   without the "case" block statements present. */
#define DO_ARG_2(function) \
\
/* Obtain pointers to the top two stack elements (vectors), decreasing \
   the top of stack index by one. */ \
      xv2 = stack[ tos-- ]; \
      xv1 = stack[ tos ]; \
\
/* Loop to access each vector element, obtaining the value of the \
   first argument and checking that it is not bad. */ \
      for ( point = 0; point < npoint; point++ ) { \
         if ( ( x1 = xv1[ point ] ) != AST__BAD ) { \
\
/* Also obtain a pointer to the element which is to receive the \
   result. */ \
            y = xv1 + point; \
\
/* Obtain the value of the second argument, again checking that it is \
   not bad. */ \
            if ( ( x2 = xv2[ point ] ) != AST__BAD ) { \
\
/* Perform the processing, which uses the two argument values and then \
   assigns the result to the appropriate top of stack element. */ \
               {function;} \
\
/* If the second argument was bad, so is the result. */ \
            } else { \
               *y = AST__BAD; \
            } \
         } \
      }

/* This macro simply wraps the one above up in a "case" block. */
#define ARG_2(oper,function) \
   case oper: \
      DO_ARG_2(function) \
      break;

/* Two-argument boolean operation. */
/* ------------------------------- */
/* This macro is similar in function to ARG_2 above, except that no
   checks are made for bad argument values. It is intended for use with
   boolean functions where bad values are handled explicitly. */
#define ARG_2B(oper,function) \
\
/* Test for the required opcode value. */ \
   case oper: \
\
/* Obtain pointers to the top two stack elements (vectors), decreasing \
   the top of stack index by one. */ \
      xv2 = stack[ tos-- ]; \
      xv1 = stack[ tos ]; \
\
/* Loop to access each vector element, obtaining the value of both \
   arguments and a pointer to the element which is to receive the \
   result. */ \
      for ( point = 0; point < npoint; point++ ) { \
         x1 = xv1[ point ]; \
         x2 = xv2[ point ]; \
         y = xv1 + point; \
\
/* Perform the processing, which uses the two argument values and then \
   assigns the result to the appropriate top of stack element. */ \
         {function;} \
      } \
\
/* Break out of the "case" block. */ \
      break;

/* Three-argument boolean operation. */
/* --------------------------------- */
/* This macro is similar in function to ARG_2B above, except that it
   takes three values of the stack and puts one back. It performs no
   checks for bad values. */
#define ARG_3B(oper,function) \
\
/* Test for the required opcode value. */ \
   case oper: \
\
/* Obtain pointers to the top three stack elements (vectors), decreasing \
   the top of stack index by two. */ \
      xv3 = stack[ tos-- ]; \
      xv2 = stack[ tos-- ]; \
      xv1 = stack[ tos ]; \
\
/* Loop to access each vector element, obtaining the value of all 3 \
   arguments and a pointer to the element which is to receive the \
   result. */ \
      for ( point = 0; point < npoint; point++ ) { \
         x1 = xv1[ point ]; \
         x2 = xv2[ point ]; \
         x3 = xv3[ point ]; \
         y = xv1 + point; \
\
/* Perform the processing, which uses the three argument values and then \
   assigns the result to the appropriate top of stack element. */ \
         {function;} \
      } \
\
/* Break out of the "case" block. */ \
      break;

/* Define arithmetic operations. */
/* ============================= */
/* We now define macros for performing some of the arithmetic
   operations we will require in a "safe" way - i.e. trapping numerical
   problems such as overflow and invalid arguments and translating them
   into the AST__BAD value. */

/* Absolute value. */
/* --------------- */
/* This is just shorthand. */
#define ABS(x) ( ( (x) >= 0.0 ) ? (x) : -(x) )

/* Integer part. */
/* ------------- */
/* This implements rounding towards zero without involving conversion
   to an integer (which could overflow). */
#define INT(x) ( ( (x) >= 0.0 ) ? floor( (x) ) : ceil( (x) ) )

/* Trap maths overflow. */
/* -------------------- */
/* This macro calls a C maths library function and checks for overflow
   in the result. */
#define CATCH_MATHS_OVERFLOW(function) \
   ( \
\
/* Clear the "errno" value. */ \
      errno = 0, \
\
/* Evaluate the function. */ \
      result = (function), \
\
/* Check if "errno" and the returned result indicate overflow and \
   return the appropriate result. */ \
      ( ( errno == ERANGE ) && ( ABS( result ) == HUGE_VAL ) ) ? AST__BAD : \
                                                                 result \
   )

/* Trap maths errors. */
/* ------------------ */
/* This macro is similar to the one above, except that it also checks
   for domain errors (i.e. invalid argument values). */
#define CATCH_MATHS_ERROR(function) \
   ( \
\
/* Clear the "errno" value. */ \
      errno = 0, \
\
/* Evaluate the function. */ \
      result = (function), \
\
/* Check if "errno" and the returned result indicate a domain error or \
   overflow and return the appropriate result. */ \
      ( ( errno == EDOM ) || \
        ( ( errno == ERANGE ) && ( ABS( result ) == HUGE_VAL ) ) ) ? \
                                 AST__BAD : result \
   )

/* Tri-state boolean OR. */
/* --------------------- */
/* This evaluates a boolean OR using tri-state logic. For example,
   "a||b" may evaluate to 1 if "a" is bad but "b" is non-zero, so that
   the normal rules of bad value propagation do not apply. */
#define TRISTATE_OR(x1,x2) \
\
/* Test if the first argument is bad. */ \
   ( (x1) == AST__BAD ) ? ( \
\
/* If so, test the second argument. */ \
      ( ( (x2) == 0.0 ) || ( (x2) == AST__BAD ) ) ? AST__BAD : 1.0 \
   ) : ( \
\
/* Test if the second argument is bad. */ \
      ( (x2) == AST__BAD ) ? ( \
\
/* If so, test the first argument. */ \
         ( (x1) == 0.0 ) ? AST__BAD : 1.0 \
\
/* If neither argument is bad, use the normal OR operator. */ \
      ) : ( \
         ( (x1) != 0.0 ) || ( (x2) != 0.0 ) \
      ) \
   )

/* Tri-state boolean AND. */
/* ---------------------- */
/* This evaluates a boolean AND using tri-state logic. */
#define TRISTATE_AND(x1,x2) \
\
/* Test if the first argument is bad. */ \
   ( (x1) == AST__BAD ) ? ( \
\
/* If so, test the second argument. */ \
      ( (x2) != 0.0 ) ? AST__BAD : 0.0 \
   ) : ( \
\
/* Test if the second argument is bad. */ \
      ( (x2) == AST__BAD ) ? ( \
\
/* If so, test the first argument. */ \
         ( (x1) != 0.0 ) ? AST__BAD : 0.0 \
\
/* If neither argument is bad, use the normal AND operator. */ \
      ) : ( \
         ( (x1) != 0.0 ) && ( (x2) != 0.0 ) \
      ) \
   )

/* Safe addition. */
/* -------------- */
/* This macro performs addition while avoiding possible overflow. */
#define SAFE_ADD(x1,x2) ( \
\
/* Test if the first argument is non-negative. */ \
   ( (x1) >= 0.0 ) ? ( \
\
/* If so, then we can perform addition if the second argument is \
   non-positive. Otherwise, we must calculate the most positive safe \
   second argument value that can be added and test for this (the test \
   itself is safe against overflow). */ \
      ( ( (x2) <= 0.0 ) || ( ( (DBL_MAX) - (x1) ) >= (x2) ) ) ? ( \
\
/* Perform addition if it is safe, otherwise return AST__BAD. */ \
         (x1) + (x2) \
      ) : ( \
         AST__BAD \
      ) \
\
/* If the first argument is negative, then we can perform addition if \
   the second argument is non-negative. Otherwise, we must calculate the \
   most negative second argument value that can be added and test for \
   this (the test itself is safe against overflow). */ \
   ) : ( \
      ( ( (x2) >= 0.0 ) || ( ( (DBL_MAX) + (x1) ) >= -(x2) ) ) ? ( \
\
/* Perform addition if it is safe, otherwise return AST__BAD. */ \
         (x1) + (x2) \
      ) : ( \
         AST__BAD \
      ) \
   ) \
)

/* Safe subtraction. */
/* ----------------- */
/* This macro performs subtraction while avoiding possible overflow. */
#define SAFE_SUB(x1,x2) ( \
\
/* Test if the first argument is non-negative. */ \
   ( (x1) >= 0.0 ) ? ( \
\
/* If so, then we can perform subtraction if the second argument is \
   also non-negative. Otherwise, we must calculate the most negative safe \
   second argument value that can be subtracted and test for this (the \
   test itself is safe against overflow). */ \
      ( ( (x2) >= 0.0 ) || ( ( (DBL_MAX) - (x1) ) >= -(x2) ) ) ? ( \
\
/* Perform subtraction if it is safe, otherwise return AST__BAD. */ \
         (x1) - (x2) \
      ) : ( \
         AST__BAD \
      ) \
\
/* If the first argument is negative, then we can perform subtraction \
   if the second argument is non-positive. Otherwise, we must calculate \
   the most positive second argument value that can be subtracted and \
   test for this (the test itself is safe against overflow). */ \
   ) : ( \
      ( ( (x2) <= 0.0 ) || ( ( (DBL_MAX) + (x1) ) >= (x2) ) ) ? ( \
\
/* Perform subtraction if it is safe, otherwise return AST__BAD. */ \
         (x1) - (x2) \
      ) : ( \
         AST__BAD \
      ) \
   ) \
)

/* Safe multiplication. */
/* -------------------- */
/* This macro performs multiplication while avoiding possible overflow. */
#define SAFE_MUL(x1,x2) ( \
\
/* Multiplication is safe if the absolute value of either argument is \
   unity or less. Otherwise, we must use the first argument to calculate \
   the maximum absolute value that the second argument may have and test \
   for this (the test itself is safe against overflow). */ \
   ( ( ( abs1 = ABS( (x1) ) ) <= 1.0 ) || \
     ( ( abs2 = ABS( (x2) ) ) <= 1.0 ) || \
     ( ( (DBL_MAX) / abs1 ) >= abs2 ) ) ? ( \
\
/* Perform multiplication if it is safe, otherwise return AST__BAD. */ \
      (x1) * (x2) \
   ) : ( \
      AST__BAD \
   ) \
)

/* Safe division. */
/* -------------- */
/* This macro performs division while avoiding possible overflow. */
#define SAFE_DIV(x1,x2) ( \
\
/* Division is unsafe if the second argument is zero. Otherwise, it is \
   safe if the abolute value of the second argument is unity or \
   more. Otherwise, we must use the second argument to calculate the \
   maximum absolute value that the first argument may have and test for \
   this (the test itself is safe against overflow). */ \
   ( ( (x2) != 0.0 ) && \
     ( ( ( abs2 = ABS( (x2) ) ) >= 1.0 ) || \
       ( ( (DBL_MAX) * abs2 ) >= ABS( (x1) ) ) ) ) ? ( \
\
/* Perform division if it is safe, otherwise return AST__BAD. */ \
      (x1) / (x2) \
   ) : ( \
      AST__BAD \
   ) \
)

/* Bit-shift operation. */
/* -------------------- */
/* This macro shifts the bits in a double value a specified number of
   places to the left, which simply corresponds to multiplying by the
   appropriate power of two. */
#define SHIFT_BITS(x1,x2) ( \
\
/* Decompose the value into a normalised fraction and a power of 2. */ \
   frac = frexp( (x1), &expon ), \
\
/* Calculate the new power of 2 which should apply after the shift, \
   rounding towards zero to give an integer value. */ \
   newexp = INT( (x2) ) + (double) expon, \
\
/* If the new exponent is too negative to convert to an integer, then \
   the result must underflow to zero. */ \
   ( newexp < (double) -INT_MAX ) ? ( \
      0.0 \
\
/* Otherwise, if it is too positive to convert to an integer, then the \
   result must overflow, unless the normalised fraction is zero. */ \
   ) : ( ( newexp > (double) INT_MAX ) ? ( \
      ( frac == 0.0 ) ? 0.0 : AST__BAD \
\
/* Otherwise, convert the new exponent to an integer and apply \
   it. Trap any overflow which may still occur. */ \
   ) : ( \
      CATCH_MATHS_OVERFLOW( ldexp( frac, (int) newexp ) ) \
   ) ) \
)

/* Two-argument bit-wise boolean operation. */
/* ---------------------------------------- */
/* This macro expands to code which performs a bit-wise boolean
   operation on a pair of arguments and assigns the result to the
   variable "result". It operates on floating point (double) values,
   which are regarded as if they are fixed-point binary numbers with
   negative values expressed in twos-complement notation. This means that
   it delivers the same results for integer values as the normal
   (integer) C bit-wise operations. However, it will also operate on the
   fraction bits of floating point numbers. It also offers greater
   precision (the first 53 or so significant bits of the result being
   preserved for typical IEEE floating point implementations). */
#define BIT_OPER(oper,x1,x2) \
\
/* Convert each argument to a normalised fraction in the range \
   [0.5,1.0) and a power of two exponent, removing any sign \
   information. */ \
   frac1 = frexp( ABS( (x1) ), &expon1 ); \
   frac2 = frexp( ABS( (x2) ), &expon2 ); \
\
/* Set "expon" to be the larger of the two exponents. If the two \
   exponents are not equal, divide the fraction with the smaller exponent \
   by 2 to the power of the exponent difference. This gives both \
   fractions the same effective exponent (although one of them may no \
   longer be normalised). Note that overflow is avoided because all \
   numbers remain less than 1.0, but underflow may occur. */ \
   expon = expon1; \
   if ( expon2 > expon1 ) { \
      expon = expon2; \
      frac1 = ldexp( frac1, expon1 - expon ); \
   } else if ( expon1 > expon2 ) { \
      frac2 = ldexp( frac2, expon2 - expon ); \
   } \
\
/* If either of the original arguments is negative, we now subtract \
   the corresponding fraction from 2.0. If we think of the fraction as \
   represented in fixed-point binary notation, this corresponds to \
   converting negative numbers into the twos-complement form normally used \
   for integers (the sign bit being the bit with value 1) instead \
   of having a separate sign bit as for floating point numbers. \
\
   Note that one of the fractions may have underflowed during the \
   scaling above. In that case (if the original argument was negative), \
   we must subtract the value "eps" (= 2.0 * DBL_EPSILON) from 2.0 \
   instead, so that we produce the largest number less than 2.0. In \
   twos-complement notation this represents the smallest possible \
   negative number and corresponds to extending the sign bit of the \
   original number up into more significant bits. This causes all bits to \
   be set as we require (rather than all being clear if the underflow \
   is simply ignored). */ \
   if ( (x1) < 0.0 ) frac1 = 2.0 - ( ( frac1 > eps ) ? frac1 : eps ); \
   if ( (x2) < 0.0 ) frac2 = 2.0 - ( ( frac2 > eps ) ? frac2 : eps ); \
\
/* We now extract the bits from the fraction values into integer \
   variables so that we may perform bit-wise operations on them. However, \
   since a double may be longer than any available integer, we may \
   have to handle several successive blocks of bits individually. */ \
\
/* Extract the first block of bits by scaling by the required power of \
   2 to shift the required bits to the left of the binary point. Then \
   extract the integer part. Note that this initial shift is one bit less \
   than the number of bits in an unsigned long, because we have \
   introduced an extra sign bit. */ \
   frac1 *= scale1; \
   frac2 *= scale1; \
   b1 = (unsigned long) frac1; \
   b2 = (unsigned long) frac2; \
\
/* Perform the required bit-wise operation on the extracted blocks of \
   bits. */ \
   b = b1 oper b2; \
\
/* Extract the sign bit from this initial result. This determines \
   whether the final result bit pattern should represent a negative \
   floating point number. */ \
   neg = b & signbit; \
\
/* Initialise the floating point result by setting it to the integer \
   result multipled by the reciprocal of the scale factor used to shift \
   the bits above. This returns the result bits to their correct \
   significance. */ \
   unscale = rscale1; \
   result = (double) b * unscale; \
\
/* We now loop to extract and process further blocks of bits (if \
   present). The number of blocks is determined by the relative lengths \
   of a double and an unsigned long. In practice, some bits of the double \
   will be used by its exponent, so the last block may be incomplete and \
   will simply be padded with zeros. */ \
   for ( iblock = 1; iblock < nblock; iblock++ ) { \
\
/* Subtract the integer part (which has already been processed) from \
   each fraction, to leave the bits which remain to be processed. Then \
   multiply by a scale factor to shift the next set of bits to the left \
   of the binary point. This time, we use as many bits as will fit into \
   an unsigned long. */ \
      frac1 = ( frac1 - (double) b1 ) * scale; \
      frac2 = ( frac2 - (double) b2 ) * scale; \
\
/* Extract the integer part, which contains the required bits. */ \
      b1 = (unsigned long) frac1; \
      b2 = (unsigned long) frac2; \
\
/* Perform the required bit-wise operation on the extracted blocks of \
   bits. */ \
      b = b1 oper b2; \
\
/* Update the result floating point value by adding the new integer \
   result multiplied by a scale factor to return the bits to their \
   original significance. */ \
      unscale *= rscale; \
      result += (double) b * unscale; \
   } \
\
/* If the (normalised fraction) result represents a negative number, \
   then subtract 2.0 from it (equivalent to subtracting it from 2 and \
   negating the result). This converts back to using a separate sign bit \
   instead of twos-complement notation. */ \
   if ( neg ) result -= 2.0; \
\
/* Scale by the required power of 2 to remove the initial \
   normalisation applied and assign the result to the "result" \
   variable. */ \
   result = ldexp( result, expon )

/* Gaussian random number. */
/* ----------------------- */
/* This macro expands to code which assigns a pseudo-random value to
   the "result" variable. The value is drawn from a Gaussian distribution
   with mean "x1" and standard deviation "ABS(x2)". */
#define GAUSS(x1,x2) \
\
/* Loop until a satisfactory result is obtained. */ \
   do { \
\
/* Obtain a value drawn from a standard Gaussian distribution. */ \
      ran = Gauss( rcontext, status ); \
\
/* Multiply by "ABS(x2)", trapping possible overflow. */ \
      result = ABS( (x2) ); \
      result = SAFE_MUL( ran, result ); \
\
/* If OK, add "x1", again trapping possible overflow. */ \
      if ( result != AST__BAD ) result = SAFE_ADD( result, (x1) ); \
\
/* Continue generating values until one is found which does not cause \
   overflow. */ \
   } while ( result == AST__BAD );

/* Implement the stack-based arithmetic. */
/* ===================================== */
/* Initialise the top of stack index and constant counter. */
      tos = -1;
      icon = 0;

/* Determine the number of opcodes to be processed and loop to process
   them, executing the appropriate "case" block for each one. */
      ncode = code[ 0 ];
      for ( icode = 1; icode <= ncode; icode++ ) {
         switch ( (Oper) code[ icode ] ) {

/* Ignore any null opcodes (which shouldn't occur). */
            case OP_NULL: break;

/* Otherwise, perform the required vector operation on the stack... */

/* User-supplied constants and variables. */
/* -------------------------------------- */
/* Loading a constant involves incrementing the constant count and
   assigning the next constant's value to the top of stack element. */
            ARG_0( OP_LDCON,    value = con[ icon++ ], *y = value )

/* Loading a variable involves obtaining the variable's index by
   consuming a constant (as above), and then copying the variable's
   values into the top of stack element. */
            ARG_0( OP_LDVAR,    ivar = (int) ( con[ icon++ ] + 0.5 ),
                                *y = ptr_in[ ivar ][ point ] )

/* System constants. */
/* ----------------- */
/* Loading a "bad" value simply means assigning AST__BAD to the top of
   stack element. */
            ARG_0( OP_LDBAD,    ;, *y = AST__BAD )

/* The following load constants associated with the (double) floating
   point representation into the top of stack element. */
            ARG_0( OP_LDDIG,    ;, *y = (double) DBL_DIG )
            ARG_0( OP_LDEPS,    ;, *y = DBL_EPSILON )
            ARG_0( OP_LDMAX,    ;, *y = DBL_MAX )
            ARG_0( OP_LDMAX10E, ;, *y = (double) DBL_MAX_10_EXP )
            ARG_0( OP_LDMAXE,   ;, *y = (double) DBL_MAX_EXP )
            ARG_0( OP_LDMDIG,   ;, *y = (double) DBL_MANT_DIG )
            ARG_0( OP_LDMIN,    ;, *y = DBL_MIN )
            ARG_0( OP_LDMIN10E, ;, *y = (double) DBL_MIN_10_EXP )
            ARG_0( OP_LDMINE,   ;, *y = (double) DBL_MIN_EXP )
            ARG_0( OP_LDRAD,    ;, *y = (double) FLT_RADIX )
            ARG_0( OP_LDRND,    ;, *y = (double) FLT_ROUNDS )

/* Mathematical constants. */
/* ----------------------- */
/* The following load mathematical constants into the top of stack
   element. */
            ARG_0( OP_LDE,      value = exp( 1.0 ), *y = value )
            ARG_0( OP_LDPI,     ;, *y = pi )

/* Functions with one argument. */
/* ---------------------------- */
/* The following simply evaluate a function of the top of stack
   element and assign the result to the same element. */
            ARG_1( OP_ABS,      *y = ABS( x ) )
            ARG_1( OP_ACOS,     *y = ( ABS( x ) <= 1.0 ) ?
                                     acos( x ) : AST__BAD )
            ARG_1( OP_ACOSD,    *y = ( ABS( x ) <= 1.0 ) ?
                                     acos( x ) * r2d : AST__BAD )
            ARG_1( OP_ACOSH,    *y = ( x < 1.0 ) ? AST__BAD :
                                     ( ( x > safe_sq ) ? log( x ) + log2 :
                                       log( x + sqrt( x * x - 1.0 ) ) ) )
            ARG_1( OP_ACOTH,    *y = ( ABS( x ) <= 1.0 ) ? AST__BAD :
                                     0.5 * ( log( ( x + 1.0 ) /
                                                  ( x - 1.0 ) ) ) )
            ARG_1( OP_ACSCH,    *y = ( ( x == 0.0 ) ? AST__BAD :
                                       ( sign = ( x >= 0.0 ), x = ABS( x ),
                                       ( sign ? 1.0 : -1.0 ) *
                                       ( ( x < rsafe_sq ) ? log2 - log( x ) :
                                         ( x = 1.0 / x,
                                       log( x + sqrt( x * x + 1.0 ) ) ) ) ) ) )
            ARG_1( OP_ASECH,    *y = ( ( x <= 0 ) || ( x > 1.0 ) ) ? AST__BAD :
                                       ( ( x < rsafe_sq ) ? log2 - log( x ) :
                                         ( x = 1.0 / x,
                                           log( x + sqrt( x * x - 1.0 ) ) ) ) )
            ARG_1( OP_ASIN,     *y = ( ABS( x ) <= 1.0 ) ?
                                     asin( x ) : AST__BAD )
            ARG_1( OP_ASIND,    *y = ( ABS( x ) <= 1.0 ) ?
                                     asin( x ) * r2d : AST__BAD )
            ARG_1( OP_ASINH,    *y = ( sign = ( x >= 0.0 ), x = ABS( x ),
                                       ( sign ? 1.0 : -1.0 ) *
                                       ( ( x > safe_sq ) ? log( x ) + log2 :
                                         log( x + sqrt( x * x + 1.0 ) ) ) ) )
            ARG_1( OP_ATAN,     *y = atan( x ) )
            ARG_1( OP_ATAND,    *y = atan( x ) * r2d )
            ARG_1( OP_ATANH,    *y = ( ABS( x ) >= 1.0 ) ? AST__BAD :
                                     0.5 * ( log( ( 1.0 + x ) /
                                                  ( 1.0 - x ) ) ) )
            ARG_1( OP_CEIL,     *y = ceil( x ) )
            ARG_1( OP_COS,      *y = cos( x ) )
            ARG_1( OP_COSD,     *y = cos( x * d2r ) )
            ARG_1( OP_COSH,     *y = CATCH_MATHS_OVERFLOW( cosh( x ) ) )
            ARG_1( OP_COTH,     *y = ( x = tanh( x ), SAFE_DIV( 1.0, x ) ) )
            ARG_1( OP_CSCH,     *y = ( x = CATCH_MATHS_OVERFLOW( sinh( x ) ),
                                       ( x == AST__BAD ) ?
                                       0.0 : SAFE_DIV( 1.0, x ) ) )
            ARG_1( OP_EXP,      *y = CATCH_MATHS_OVERFLOW( exp( x ) ) )
            ARG_1( OP_FLOOR,    *y = floor( x ) )
            ARG_1( OP_INT,      *y = INT( x ) )
            ARG_1B( OP_ISBAD,   *y = ( x == AST__BAD ) )
            ARG_1( OP_LOG,      *y = ( x > 0.0 ) ? log( x ) : AST__BAD )
            ARG_1( OP_LOG10,    *y = ( x > 0.0 ) ? log10( x ) : AST__BAD )
            ARG_1( OP_NINT,     *y = ( x >= 0 ) ?
                                     floor( x + 0.5 ) : ceil( x - 0.5 ) )
            ARG_1( OP_POISS,    *y = Poisson( rcontext, x, status ) )
            ARG_1( OP_SECH,     *y = ( x = CATCH_MATHS_OVERFLOW( cosh( x ) ),
                                       ( x == AST__BAD ) ? 0.0 : 1.0 / x ) )
            ARG_1( OP_SIN,      *y = sin( x ) )
            ARG_1( OP_SINC,     *y = ( x == 0.0 ) ? 1.0 : sin( x ) / x )
            ARG_1( OP_SIND,     *y = sin( x * d2r ) )
            ARG_1( OP_SINH,     *y = CATCH_MATHS_OVERFLOW( sinh( x ) ) )
            ARG_1( OP_SQR,      *y = SAFE_MUL( x, x ) )
            ARG_1( OP_SQRT,     *y = ( x >= 0.0 ) ? sqrt( x ) : AST__BAD )
            ARG_1( OP_TAN,      *y = CATCH_MATHS_OVERFLOW( tan( x ) ) )
            ARG_1( OP_TAND,     *y = tan( x * d2r ) )
            ARG_1( OP_TANH,     *y = tanh( x ) )

/* Functions with two arguments. */
/* ----------------------------- */
/* These evaluate a function of the top two entries on the stack. */
            ARG_2( OP_ATAN2,    *y = atan2( x1, x2 ) )
            ARG_2( OP_ATAN2D,   *y = atan2( x1, x2 ) * r2d )
            ARG_2( OP_DIM,      *y = ( x1 > x2 ) ? x1 - x2 : 0.0 )
            ARG_2( OP_GAUSS,    GAUSS( x1, x2 ); *y = result )
            ARG_2( OP_MOD,      *y = ( x2 != 0.0 ) ?
                                     fmod( x1, x2 ) : AST__BAD )
            ARG_2( OP_POW,      *y = CATCH_MATHS_ERROR( pow( x1, x2 ) ) )
            ARG_2( OP_RAND,     ran = Rand( rcontext, status );
                                *y = x1 * ran + x2 * ( 1.0 - ran ); )
            ARG_2( OP_SIGN,     *y = ( ( x1 >= 0.0 ) == ( x2 >= 0.0 ) ) ?
                                     x1 : -x1 )

/* Functions with three arguments. */
/* ------------------------------- */
/* These evaluate a function of the top three entries on the stack. */
            ARG_3B( OP_QIF,     *y = ( ( x1 ) ? ( x2 ) : ( x3 ) ) )


/* Functions with variable numbers of arguments. */
/* --------------------------------------------- */
/* These operations take a variable number of arguments, the actual
   number being determined by consuming a constant. We then loop to
   perform a 2-argument operation on the stack (as above) the required
   number of times. */
            case OP_MAX:
               narg = (int) ( con[ icon++ ] + 0.5 );
               for ( iarg = 0; iarg < ( narg - 1 ); iarg++ ) {
                  DO_ARG_2( *y = ( x1 >= x2 ) ? x1 : x2 )
               }
               break;
            case OP_MIN:
               narg = (int) ( con[ icon++ ] + 0.5 );
               for ( iarg = 0; iarg < ( narg - 1 ); iarg++ ) {
                  DO_ARG_2( *y = ( x1 <= x2 ) ? x1 : x2 )
               }
               break;

/* Unary arithmetic operators. */
/* --------------------------- */
            ARG_1( OP_NEG,      *y = -x )

/* Unary boolean operators. */
/* ------------------------ */
            ARG_1( OP_NOT,      *y = ( x == 0.0 ) )

/* Binary arithmetic operators. */
/* ---------------------------- */
            ARG_2( OP_ADD,      *y = SAFE_ADD( x1, x2 ) )
            ARG_2( OP_SUB,      *y = SAFE_SUB( x1, x2 ) )
            ARG_2( OP_MUL,      *y = SAFE_MUL( x1, x2 ) )
            ARG_2( OP_DIV ,     *y = SAFE_DIV( x1, x2 ) )

/* Bit-shift operators. */
/* -------------------- */
            ARG_2( OP_SHFTL,    *y = SHIFT_BITS( x1, x2 ) )
            ARG_2( OP_SHFTR,    *y = SHIFT_BITS( x1, -x2 ) )

/* Relational operators. */
/* --------------------- */
            ARG_2( OP_EQ,       *y = ( x1 == x2 ) )
            ARG_2( OP_GE,       *y = ( x1 >= x2 ) )
            ARG_2( OP_GT,       *y = ( x1 > x2 ) )
            ARG_2( OP_LE,       *y = ( x1 <= x2 ) )
            ARG_2( OP_LT,       *y = ( x1 < x2 ) )
            ARG_2( OP_NE,       *y = ( x1 != x2 ) )

/* Bit-wise operators. */
/* ------------------- */
            ARG_2( OP_BITOR,    BIT_OPER( |, x1, x2 ); *y = result )
            ARG_2( OP_BITXOR,   BIT_OPER( ^, x1, x2 ); *y = result )
            ARG_2( OP_BITAND,   BIT_OPER( &, x1, x2 ); *y = result )

/* Binary boolean operators. */
/* ------------------------- */
            ARG_2B( OP_AND,     *y = TRISTATE_AND( x1, x2 ) )
            ARG_2( OP_EQV,      *y = ( ( x1 != 0.0 ) == ( x2 != 0.0 ) ) )
            ARG_2B( OP_OR,      *y = TRISTATE_OR( x1, x2 ) )
            ARG_2( OP_XOR,      *y = ( ( x1 != 0.0 ) != ( x2 != 0.0 ) ) )
         }
      }
   }

/* When all opcodes have been processed, the result of the function
   evaluation will reside in the lowest stack entry - i.e. the output
   array. */

/* Free the workspace arrays. */
   work = astFree( work );
   stack = astFree( stack );

/* Undefine macros local to this function. */
#undef ARG_0
#undef ARG_1
#undef ARG_1B
#undef DO_ARG_2
#undef ARG_2
#undef ARG_2B
#undef ABS
#undef INT
#undef CATCH_MATHS_OVERFLOW
#undef CATCH_MATHS_ERROR
#undef TRISTATE_OR
#undef TRISTATE_AND
#undef SAFE_ADD
#undef SAFE_SUB
#undef SAFE_MUL
#undef SAFE_DIV
#undef SHIFT_BITS
#undef BIT_OPER
#undef GAUSS
}

static void EvaluationSort( const double con[], int nsym, int symlist[],
                            int **code, int *stacksize, int *status ) {
/*
*  Name:
*     EvaluationSort

*  Purpose:
*     Perform an evaluation-order sort on parsed expression symbols.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     void EvaluationSort( const double con[], int nsym, int symlist[],
*                          int **code, int *stacksize, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     This function sorts a sequence of numbers representing symbols
*     identified in an expression. The symbols (i.e. the expression syntax)
*     must have been fully validated beforehand, as no validation is
*     performed here.
*
*     The symbols are sorted into the order in which corresponding
*     operations must be performed on a push-down arithmetic stack in order
*     to evaluate the expression. Operation codes (opcodes), as defined in
*     the "Oper" enum, are then substituted for the symbol numbers.

*  Parameters:
*     con
*        Pointer to an array of double containing the set of constants
*        generated while parsing the expression (these are required in order
*        to determine the number of arguments associated with functions which
*        take a variable number of arguments).
*     nsym
*        The number of symbols identified while parsing the expression.
*     symlist
*        Pointer to an array of int, with "nsym" elements. On entry, this
*        should contain the indices in the static "symbol" array of the
*        symbols identified while parsing the expression. On exit, the
*        contents are undefined.
*     code
*        Address of a pointer which will be set to point at a dynamically
*        allocated array of int containing the set of opcodes (cast to int)
*        produced by this function. The first element of this array will
*        contain a count of the number of opcodes which follow.
*
*        The allocated space must be freed by the caller (using astFree) when
*        no longer required.
*     stacksize
*        Pointer to an int in which to return the size of the push-down stack
*        required to evaluate the expression using the returned opcodes.
*     status
*        Pointer to the inherited status variable.

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

/* Local Variables: */
   int flush;                    /* Flush parenthesised symbol sequence? */
   int icon;                     /* Input constant counter */
   int isym;                     /* Input symbol counter */
   int ncode;                    /* Number of opcodes generated */
   int nstack;                   /* Evaluation stack size */
   int push;                     /* Push a new symbol on to stack? */
   int sym;                      /* Variable for symbol number */
   int tos;                      /* Top of sort stack index */

/* Initialise */
   *code = NULL;
   *stacksize = 0;

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

/* Further initialisation. */
   flush = 0;
   icon = 0;
   isym = 0;
   ncode = 0;
   nstack = 0;
   tos = -1;

/* Loop to generate output opcodes until the sort stack is empty and
   there are no further symbols to process, or an error is detected.  */
   while ( astOK && ( ( tos > -1 ) || ( isym < nsym ) ) ) {

/* Decide whether to push a symbol on to the sort stack (which
   "diverts" it so that higher-priority symbols can be output), or to pop
   the top symbol off the sort stack and send it to the output
   stream... */

/* We must push a symbol on to the sort stack if the stack is
   currently empty. */
      if ( tos == -1 ) {
         push = 1;

/* We must pop the top symbol off the sort stack if there are no more
   input symbols to process. */
      } else if ( isym >= nsym ) {
         push = 0;

/* If the sort stack is being flushed to complete the evaluation of a
   parenthesised expression, then the top symbol (which will be the
   opening parenthesis or function call) must be popped. This is only
   done once, so reset the "flush" flag before the next loop. */
      } else if ( flush ) {
         push = 0;
         flush = 0;

/* In all other circumstances, we must push a symbol on to the sort
   stack if its evaluation priority (seen from the left) is higher than
   that of the current top of stack symbol (seen from the right). This
   means it will eventually be sent to the output stream ahead of the
   current top of stack symbol. */
      } else {
         push = ( symbol[ symlist[ isym ] ].leftpriority >
                  symbol[ symlist[ tos ] ].rightpriority );
      }

/* If a symbol is being pushed on to the sort stack, then get the next
   input symbol which is to be used. */
      if ( push ) {
         sym = symlist[ isym++ ];

/* If the symbol decreases the parenthesis level (a closing
   parenthesis), then all the sort stack entries down to the symbol which
   opened the current level of parenthesis (the matching opening
   parenthesis or function call) will already have been sent to the
   output stream as a consequence of the evaluation priority defined for
   a closing parenthesis in the symbol data. The opening parenthesis (or
   function call) must next be flushed from the sort stack, so set the
   "flush" flag which is interpreted on the next loop. Ignore the current
   symbol, which cancels with the opening parenthesis on the stack. */
         if ( symbol[ sym ].parincrement < 0 ) {
            flush = 1;

/* All other symbols are pushed on to the sort stack. The stack
   occupies that region of the "symlist" array from which the input
   symbol numbers have already been extracted. */
         } else {
            symlist[ ++tos ] = sym;
         }

/* If a symbol is being popped from the top of the sort stack, then
   the top of stack entry is transferred to the output stream. Obtain the
   symbol number from the stack. Increment the local constant counter if
   the associated operation will use a constant. */
      } else {
         sym = symlist[ tos-- ];
         icon += ( ( sym == symbol_ldvar ) || ( sym == symbol_ldcon ) );

/* If the output symbol does not represent a "null" operation,
   increase the size of the output opcode array to accommodate it,
   checking for errors. Note that we allocate one extra array element
   (the first) which will eventually hold a count of all the opcodes
   generated. */
         if ( symbol[ sym ].opcode != OP_NULL ) {
            *code = astGrow( *code, ncode + 2, sizeof( int ) );
            if ( astOK ) {

/* Append the new opcode to the end of this array. */
               ( *code )[ ++ncode ] = (int) symbol[ sym ].opcode;

/* Increment/decrement the counter representing the stack size
   required for evaluation of the expression.  If the symbol is a
   function with a variable number of arguments (indicated by a negative
   "nargs" entry in the symbol data table), then the change in stack size
   must be determined from the argument number stored in the constant
   table. */
               if ( symbol[ sym ].nargs >= 0 ) {
                  nstack += symbol[ sym ].stackincrement;
               } else {
                  nstack -= (int) ( con[ icon++ ] + 0.5 ) - 1;
               }

/* Note the maximum size of the stack. */
               *stacksize = ( nstack > *stacksize ) ? nstack : *stacksize;
            }
         }
      }
   }

/* If no "*code" array has been allocated, then allocate one simply to
   store the number of opcodes generated, i.e. zero (this shouldn't
   normally happen as this represents an invalid expression). */
   if ( !*code ) *code = astMalloc( sizeof( int ) );

/* If no error has occurred, store the count of opcodes generated in
   the first element of the "*code" array and re-allocate the array to
   its final size (since astGrow may have over-allocated space). */
   if ( astOK ) {
      ( *code )[ 0 ] = ncode;
      *code = astRealloc( *code, sizeof( int ) * (size_t) ( ncode + 1 ) );
   }

/* If an error occurred, free any memory that was allocated and reset
   the output values. */
   if ( !astOK ) {
      *code = astFree( *code );
      *stacksize = 0;
   }
}

static void ExtractExpressions( const char *method, const char *class,
                                int nfun, const char *fun[], int forward,
                                char ***exprs, int *status ) {
/*
*  Name:
*     ExtractExpressions

*  Purpose:
*     Extract and validate expressions.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     void ExtractExpressions( const char *method, const char *class,
*                              int nfun, const char *fun[], int forward,
*                              char ***exprs, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     This function extracts expressions from the right hand sides of a set
*     of functions. These expressions are then validated to check that they
*     are either all present, or all absent (absence indicating an undefined
*     transformation). An error is reported if anything is found to be
*     wrong.
*
*     Note that the syntax of the expressions is not checked by this function
*     (i.e. they are not compiled).

*  Parameters:
*     method
*        Pointer to a constant null-terminated character string
*        containing the name of the method that invoked this function.
*        This method name is used solely for constructing error messages.
*     class
*        Pointer to a constant null-terminated character string containing the
*        class name of the Object being processed. This name is used solely
*        for constructing error messages.
*     nfun
*        The number of functions to be analysed.
*     fun
*        Pointer to an array, with "nfun" elements, of pointers to null
*        terminated strings which contain each of the functions. These
*        strings should contain no white space.
*     forward
*        A non-zero value indicates the the MathMap's forward transformation
*        functions are being processed, while a zero value indicates processing
*        of the inverse transformation functions. This value is used solely for
*        constructing error messages.
*     exprs
*        Address in which to return a pointer to an array (with "nfun"
*        elements) of pointers to null terminated strings containing the
*        extracted expressions (i.e. this returns an array of strings).
*
*        Both the returned array of pointers, and the strings to which they
*        point, will be stored in dynamically allocated memory and should
*        be freed by the caller (using astFree) when no longer required.
*
*        If the right hand sides (including the "=" sign) of all the supplied
*        functions are absent, then this indicates an undefined transformation
*        and the returned pointer value will be NULL. An error results if
*        an "=" sign is present but no expression follows it.
*     status
*        Pointer to the inherited status variable.

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

/* Local Variables: */
   char *ex;                     /* Pointer to start of expression string */
   int ifun;                     /* Loop counter for functions */
   int iud;                      /* Index of first undefined function */
   int nud;                      /* Number of undefined expressions */

/* Initialise. */
   *exprs = NULL;

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

/* Further initialisation. */
   nud = 0;
   iud = 0;

/* Allocate and initialise memory for the returned array of pointers. */
   MALLOC_POINTER_ARRAY( *exprs, char *, nfun )

/* Loop to inspect each function in turn. */
   if ( astOK ) {
      for ( ifun = 0; ifun < nfun; ifun++ ) {

/* Search for the first "=" sign. */
         if ( ( ex = strchr( fun[ ifun ], '=' ) ) ) {

/* If found, and there are more characters after the "=" sign, then
   find the length of the expression which follows. Allocate a string to
   hold this expression, storing its pointer in the array allocated
   above. Check for errors. */
            if ( *++ex ) {
               ( *exprs )[ ifun ] = astMalloc( strlen( ex ) + (size_t) 1 );
               if ( !astOK ) break;

/* If OK, extract the expression string. */
               (void) strcpy( ( *exprs )[ ifun ], ex );

/* If an "=" sign was found but there are no characters following it,
   then there is a missing right hand side to a function, so report an
   error and quit. */
            } else {
               astError( AST__NORHS,
                         "%s(%s): Missing right hand side in expression: "
                         "\"%s\".", status,
                         method, class, fun[ ifun ] );
               astError( astStatus,
                         "Error in %s transformation function %d.", status,
                         forward ? "forward" : "inverse", ifun + 1 );
               break;
            }

/* If no "=" sign was found, then the transformation may be undefined,
   in which case each function should only contain a variable name. Count
   the number of times this happens and record the index of the first
   instance. */
         } else {
            nud++;
            if ( nud == 1 ) iud = ifun;
         }
      }
   }

/* Either all functions should have an "=" sign (in which case the
   transformation is defined), or none of them should have (in which case
   it is undefined). If some do and some don't, then report an error,
   citing the first instance of a missing "=" sign. */
   if ( astOK && ( nud != 0 ) && ( nud != nfun ) ) {
      astError( AST__NORHS,
                "%s(%s): Missing right hand side in function: \"%s\".", status,
                method, class, fun[ iud ] );
      astError( astStatus,
                "Error in %s transformation function %d.", status,
                forward ? "forward" : "inverse", iud + 1 );
   }

/* If an error occurred, or all the expressions were absent, then free any
   allocated memory and reset the output value. */
   if ( !astOK || nud ) {
      FREE_POINTER_ARRAY( *exprs, nfun )
   }
}

static void ExtractVariables( const char *method, const char *class,
                              int nfun, const char *fun[],
                              int nin, int nout, int nfwd, int ninv,
                              int forward, char ***var, int *status ) {
/*
*  Name:
*     ExtractVariables

*  Purpose:
*     Extract and validate variable names.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     void ExtractVariables( const char *method, const char *class,
*                            int nfun, const char *fun[],
*                            int nin, int nout, int nfwd, int ninv,
*                            int forward, char ***var, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     This function extracts variable names from the left hand sides of a
*     set of transformation functions belonging to a MathMap. These variable
*     names are then validated to check for correct syntax and no
*     duplication. An error is reported if anything is wrong with the
*     variable names obtained.

*  Parameters:
*     method
*        Pointer to a constant null-terminated character string
*        containing the name of the method that invoked this function.
*        This method name is used solely for constructing error messages.
*     class
*        Pointer to a constant null-terminated character string containing the
*        class name of the Object being processed. This name is used solely
*        for constructing error messages.
*     nfun
*        The number of functions to be analysed.
*     fun
*        Pointer to an array, with "nfun" elements, of pointers to null
*        terminated strings which contain each of the functions. These strings
*        are case sensitive and should contain no white space.
*
*        The first elements of this array should point to the functions that
*        define the primary input/output variables (depending on direction).
*        These should be followed by any functions which define intermediate
*        variables (taken from the set of functions which transform in the
*        opposite direction to the first ones).
*     nin
*        Number of input variables for the MathMap.
*     nout
*        Number of output variables for the MathMap.
*     nfwd
*        Number of forward transformation functions for the MathMap.
*     ninv
*        Number of inverse transformation functions for the MathMap.
*     forward
*        A non-zero value indicates the the MathMap's forward transformation
*        functions are being processed, while a zero value indicates processing
*        of the inverse transformation functions. This value, together with
*        "nin", "nout", "nfwd" and "ninv" are used solely for constructing
*        error messages.
*     var
*        Address in which to return a pointer to an array (with "nfun"
*        elements) of pointers to null terminated strings containing the
*        extracted variable names (i.e. this returns an array of strings).
*
*        Both the returned array of pointers, and the strings to which they
*        point, will be stored in dynamically allocated memory and should
*        be freed by the caller (using astFree) when no longer required.
*     status
*        Pointer to the inherited status variable.

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

/* Local Variables: */
   char *duser1;                 /* Transformation direction for function */
   char *duser2;                 /* Transformation direction for function */
   char c;                       /* Extracted character */
   int i1;                       /* Loop counter for detecting duplicates */
   int i2;                       /* Loop counter for detecting duplicates */
   int i;                        /* Loop counter for characters */
   int iend;                     /* Last character index in parsed name */
   int ifun;                     /* Loop counter for functions */
   int iuser1;                   /* Function number as known to the user */
   int iuser2;                   /* Function number as known to the user */
   int nc;                       /* Character count */
   int nextra;                   /* Number of intermediate functions */
   int nprimary;                 /* Number of primary input/output variables */

/* Initialise. */
   *var = NULL;

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

/* Obtain the number of primary input/output variables, depending on
   the direction of the coordinate transformation. */
   nprimary = ( forward ? nin : nout );

/* Deterine the number of extra (intermediate) functions that come
   before these primary ones. These affect the numbering of
   transformation functions as known to the user, and must be accounted
   for when reporting error messages. */
   nextra = ( forward ? ninv - nin : nfwd - nout );

/* Allocate and initialise memory for the returned array of pointers. */
   MALLOC_POINTER_ARRAY( *var, char *, nfun )

/* Loop to process each function in turn. */
   if ( astOK ) {
      for ( ifun = 0; ifun < nfun; ifun++ ) {

/* Count the number of characters appearing before the "=" sign (or in
   the entire string if the "=" is absent). */
         for ( nc = 0; ( c = fun[ ifun ][ nc ] ); nc++ ) if ( c == '=' ) break;

/* If no characters were counted, then report an appropriate error
   message, depending on whether the function string was entirely
   blank. */
         if ( !nc ) {
            if ( c ) {
               astError( AST__MISVN,
                         "%s(%s): No left hand side in expression: \"%s\".", status,
                         method, class, fun[ ifun ] );
            } else {
               astError( AST__MISVN,
                         "%s: Transformation function contains no variable "
                         "name.", status,
                         method );
            }
            break;
         }

/* If OK, allocate memory to hold the output string and check for
   errors. */
         ( *var )[ ifun ] = astMalloc( sizeof( char ) * (size_t) ( nc + 1 ) ) ;
         if ( !astOK ) break;

/* If OK, copy the characters before the "=" sign to the new
   string. */
         nc = 0;
         for ( i = 0; ( c = fun[ ifun ][ i ] ); i++ ) {
            if ( c == '=' ) break;
            ( *var )[ ifun ][ nc++] = c;
         }

/* Null terminate the result. */
         ( *var )[ ifun ][ nc ] = '\0';

/* Try to parse the contents of the extracted string as a name. */
         ParseName( ( *var )[ ifun ], 0, &iend, status );

/* If unsuccessful, or if all the characters were not parsed, then we
   have an invalid variable name, so report an error and quit. */
         if ( ( iend < 0 ) || ( *var )[ ifun ][ iend + 1 ] ) {
            astError( AST__VARIN,
                      "%s(%s): Variable name is invalid: \"%s\".", status,
                      method, class, ( *var )[ ifun ] );
            break;
         }
      }

/* If an error occurred above, then determine the function number, and
   the direction of the transformation of which it forms part, as known
   to the user. */
      if ( !astOK ) {
         if ( ifun < nprimary ) {
            iuser1 = ifun + 1 + nextra;
            duser1 = ( forward ? "inverse" : "forward" );
         } else {
            iuser1 = ifun + 1 - nprimary;
            duser1 = ( forward ? "forward" : "inverse" );
         }

/* Report a contextual error message. */
         astError( astStatus,
                   "Error in %s transformation function %d.", status,
                   duser1, iuser1 );
      }
   }

/* If there has been no error, loop to compare all the variable names
   with each other to detect duplication. */
   if ( astOK ) {
      for ( i1 = 1; i1 < nfun; i1++ ) {
         for ( i2 = 0; i2 < i1; i2++ ) {

/* If a duplicate variable name is found, report an error. */
            if ( !strcmp( ( *var )[ i1 ], ( *var )[ i2 ] ) ) {
               astError( AST__DUVAR,
                         "%s(%s): Duplicate definition of variable name: "
                         "\"%s\".", status,
                         method, class, ( *var )[ i1 ] );

/* For each transformation function involved, determine the function
   number and the direction of the transformation of which it forms part,
   as known to the user. */
               if ( i1 < nprimary ) {
                  iuser1 = i1 + 1 + nextra;
                  duser1 = ( forward ? "inverse" : "forward" );
               } else {
                  iuser1 = i1 + 1 - nprimary;
                  duser1 = ( forward ? "forward" : "inverse" );
               }
               if ( i2 < nprimary ) {
                  iuser2 = i2 + 1 + nextra;
                  duser2 = ( forward ? "inverse" : "forward" );
               } else {
                  iuser2 = i2 + 1 - nprimary;
                  duser2 = ( forward ? "forward" : "inverse" );
               }

/* Report a contextual error message. */
               astError( astStatus,
                         "Conflict between %s function %d and %s function %d.", status,
                         duser1, iuser1, duser2, iuser2 );
               break;
            }
         }
         if ( !astOK ) break;
      }
   }

/* If an error occurred, free any allocated memory and reset the
   output value. */
   if ( !astOK ) {
      FREE_POINTER_ARRAY( *var, nfun )
   }
}

static double Gauss( Rcontext *context, int *status ) {
/*
*  Name:
*     Gauss

*  Purpose:
*     Produce a pseudo-random sample from a standard Gaussian distribution.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     double Gauss( Rcontext *context, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     On each invocation, this function returns a pseudo-random sample drawn
*     from a standard Gaussian distribution with mean zero and standard
*     deviation unity. The Box-Muller transformation method is used.

*  Parameters:
*     context
*        Pointer to an Rcontext structure which holds the random number
*        generator's context between invocations.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     A sample from a standard Gaussian distribution.

*  Notes:
*     - The sequence of numbers returned is determined by the "seed"
*     value in the Rcontext structure supplied.
*     - If the seed value is changed, the "active" flag must also be cleared
*     so that this function can re-initiallise the Rcontext structure before
*     generating the next pseudo-random number. The "active" flag should
*     also be clear to force initialisation the first time an Rcontext
*     structure is used.
*     - This function does not perform error checking and does not generate
*     errors. It will execute even if the global error status is set.
*/

/* Local Variables: */
   double rsq;                   /* Squared radius */
   double s;                     /* Scale factor */
   double x;                     /* First result value */
   static double y;              /* Second result value */
   static int ysaved = 0;        /* Previously-saved value available? */

   LOCK_MUTEX7

/* If the random number generator context is not active, then it will
   be (re)initialised on the first invocation of Rand (below). Ensure
   that any previously-saved value within this function is first
   discarded. */
   if ( !context->active ) ysaved = 0;

/* If there is a previously-saved value available, then use it and
   mark it as no longer available. */
   if ( ysaved ) {
      x = y;
      ysaved = 0;

/* Otherwise, loop until a suitable new pair of values has been
   obtained. */
   } else {
      while ( 1 ) {

/* Loop to obtain two random values uniformly distributed inside the
   unit circle, while avoiding the origin (which maps to an infinite
   result). */
         do {
            x = 2.0 * Rand( context, status ) - 1.0;
            y = 2.0 * Rand( context, status ) - 1.0;
            rsq = x * x + y * y;
         } while ( ( rsq >= 1.0 ) || ( rsq == 0.0 ) );

/* Perform the Box-Muller transformation, checking that this will not
   produce overflow (which is extremely unlikely). If overflow would
   occur, we simply repeat the above steps with a new pair of random
   numbers. */
         s = -2.0 * log( rsq );
         if ( ( DBL_MAX * rsq ) >= s ) {
            s = sqrt( s / rsq );

/* Scale the original random values to give a pair of results. One will be
   returned and the second kept until next time. */
            x *= s;
            y *= s;
            break;
         }
      }

/* Note that a saved value is available. */
      ysaved = 1;
   }

   UNLOCK_MUTEX7

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

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

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

*  Type:
*     Private function.

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

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

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

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

/* Initialise. */
   result = 0;

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

/* Obtain a pointers to the MathMap structure. */
   this = (AstMathMap *) this_object;

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

   SIZEOF_POINTER_ARRAY( this->fwdfun, this->nfwd )
   SIZEOF_POINTER_ARRAY( this->invfun, this->ninv )
   SIZEOF_POINTER_ARRAY( this->fwdcode, this->nfwd )
   SIZEOF_POINTER_ARRAY( this->invcode, this->ninv )
   SIZEOF_POINTER_ARRAY( this->fwdcon, this->nfwd )
   SIZEOF_POINTER_ARRAY( this->invcon, this->ninv )

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

*  Type:
*     Private function.

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

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

*  Parameters:
*     this
*        Pointer to the MathMap.
*     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 MathMap, 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 MathMap. 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 */
   AstMathMap *this;             /* Pointer to the MathMap structure */
   const char *result;           /* Pointer value to return */
   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 MathMap structure. */
   this = (AstMathMap *) 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. */

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

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

/* SimpIF. */
/* ------- */
   } else if ( !strcmp( attrib, "simpif" ) ) {
      ival = astGetSimpIF( this );
      if ( astOK ) {
         (void) sprintf( getattrib_buff, "%d", ival );
         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;

}

void astInitMathMapVtab_(  AstMathMapVtab *vtab, const char *name, int *status ) {
/*
*+
*  Name:
*     astInitMathMapVtab

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

*  Type:
*     Protected function.

*  Synopsis:
*     #include "mathmap.h"
*     void astInitMathMapVtab( AstMathMapVtab *vtab, const char *name )

*  Class Membership:
*     MathMap vtab initialiser.

*  Description:
*     This function initialises the component of a virtual function
*     table which is used by the MathMap 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 */
   AstMappingVtab *mapping;      /* Pointer to Mapping component of Vtab */
   AstObjectVtab *object;        /* Pointer to Object 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 astIsAMathMap) 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->ClearSeed = ClearSeed;
   vtab->ClearSimpFI = ClearSimpFI;
   vtab->ClearSimpIF = ClearSimpIF;
   vtab->GetSeed = GetSeed;
   vtab->GetSimpFI = GetSimpFI;
   vtab->GetSimpIF = GetSimpIF;
   vtab->SetSeed = SetSeed;
   vtab->SetSimpFI = SetSimpFI;
   vtab->SetSimpIF = SetSimpIF;
   vtab->TestSeed = TestSeed;
   vtab->TestSimpFI = TestSimpFI;
   vtab->TestSimpIF = TestSimpIF;

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

   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 copy constructor, destructor and class dump function. */
   astSetCopy( vtab, Copy );
   astSetDelete( vtab, Delete );
   astSetDump( vtab, Dump, "MathMap",
               "Transformation using mathematical functions" );

/* 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 double LogGamma( double x, int *status ) {
/*
*  Name:
*     LogGamma

*  Purpose:
*     Calculate the logarithm of the gamma function.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     double LogGamma( double x, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     This function returns the natural logarithm of the gamma function
*     for real arguments x>0. It uses the approximation of Lanczos, with
*     constants from Press et al. (Numerical Recipes), giving a maximum
*     fractional error (on the gamma function) of less than 2e-10.

*  Parameters:
*     x
*        The function argument, which must be greater than zero.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     The natural logarithm of the gamma function with "x" as argument,
*     or AST__BAD if "x" is not greater than zero.

*  Notes:
*     - This function does not generate errors and does not perform error
*     reporting. It will execute even if the global error status is set.
*/

/* Local Constants: */
   const double c0 = 1.000000000190015; /* Coefficients for series sum... */
   const double c1 = 76.18009172947146;
   const double c2 = -86.50532032941677;
   const double c3 = 24.01409824083091;
   const double c4 = -1.231739572450155;
   const double c5 = 0.1208650973866179e-2;
   const double c6 = -0.5395239384953e-5;
   const double g = 5.0;

/* Local Variables: */
   double result;                /* Result value to return */
   double sum;                   /* Series sum */
   double xx;                    /* Denominator for summing series */
   static double root_twopi;     /* sqrt( 2.0 * pi ) */
   static int init = 0;          /* Initialisation performed? */

/* If initialisation has not yet been performed, calculate the
   constant required below. */
   LOCK_MUTEX3
   if ( !init ) {
      root_twopi = sqrt( 2.0 * acos( -1.0 ) );

/* Note that initialisation has been performed. */
      init = 1;
   }
   UNLOCK_MUTEX3

/* Return a bad value if "x" is not greater than zero. */
   if ( x <= 0.0 ) {
      result = AST__BAD;

/* Otherwise, form the series sum. Since we only use 6 terms, the loop
   that would normally be used has been completely unrolled here. */
   } else {
      xx = x;
      sum = c0;
      sum += c1 / ++xx;
      sum += c2 / ++xx;
      sum += c3 / ++xx;
      sum += c4 / ++xx;
      sum += c5 / ++xx;
      sum += c6 / ++xx;

/* Calculate the result. */
      result = x + g + 0.5;
      result -= ( x + 0.5 ) * log( result );
      result = log( root_twopi * sum / x ) - result;
   }

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

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 MathMap.

*  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:
*     MathMap 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 MathMap 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 MathMap 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 MathMap which is to be merged with
*        its neighbours. This should be a cloned copy of the MathMap
*        pointer contained in the array element "(*map_list)[where]"
*        (see below). This pointer will not be annulled, and the
*        MathMap it identifies will not be modified by this function.
*     where
*        Index in the "*map_list" array (below) at which the pointer
*        to the nominated MathMap 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: */
   AstMapping *new;              /* Pointer to replacement Mapping */
   AstMathMap *mathmap1;         /* Pointer to first MathMap */
   AstMathMap *mathmap2;         /* Pointer to second MathMap */
   char **fwd1;                  /* Pointer to first forward function array */
   char **fwd2;                  /* Pointer to second forward function array */
   char **inv1;                  /* Pointer to first inverse function array */
   char **inv2;                  /* Pointer to second inverse function array */
   int ifun;                     /* Loop counter for functions */
   int imap1;                    /* Index of first Mapping */
   int imap2;                    /* Index of second Mapping */
   int imap;                     /* Loop counter for Mappings */
   int invert1;                  /* Invert flag for first MathMap */
   int invert2;                  /* Invert flag for second MathMap */
   int nfwd1;                    /* No. forward functions for first MathMap */
   int nfwd2;                    /* No. forward functions for second MathMap */
   int nin1;                     /* Number input coords for first MathMap */
   int ninv1;                    /* No. inverse functions for first MathMap */
   int ninv2;                    /* No. inverse functions for second MathMap */
   int nout2;                    /* Number output coords for second MathMap */
   int result;                   /* Result value to return */
   int simplify;                 /* Mappings may simplify? */

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

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

/* Initialise variables to avoid "used of uninitialised variable"
   messages from dumb compilers. */
   mathmap1 = NULL;
   mathmap2 = NULL;
   imap1 = 0;
   imap2 = 0;
   invert1 = 0;
   invert2 = 0;
   nfwd1 = 0;
   nin1 = 0;
   ninv1 = 0;

/* MathMaps are only worth simplifying if they occur in series. */
   simplify = series;

/* If simplification appears possible, then obtain the indices of the
   nominated mapping and of the one which follows it. Check that a
   mapping exists for the second index. */
   if ( simplify ) {
      imap1 = where;
      imap2 = imap1 + 1;
      simplify = ( imap2 < *nmap );
   }

/* If OK, check whether the class of both Mappings is "MathMap" (a
   MathMap can only combine with another MathMap). */
   if ( simplify ) {
      simplify = !strcmp( astGetClass( ( *map_list )[ imap1 ] ), "MathMap" );
   }
   if ( astOK && simplify ) {
      simplify = !strcmp( astGetClass( ( *map_list )[ imap2 ] ), "MathMap" );
   }

/* If still OK, obtain pointers to the two MathMaps and the associated
   invert flag values. */
   if ( astOK && simplify ) {
      mathmap1 = (AstMathMap *) ( *map_list )[ imap1 ];
      mathmap2 = (AstMathMap *) ( *map_list )[ imap2 ];
      invert1 = ( *invert_list )[ imap1 ];
      invert2 = ( *invert_list )[ imap2 ];

/* Depending on the invert flag values, obtain the SimpFI or SimpIF
   attribute value from each MathMap and check whether they are set so as
   to permit simplification. */
      simplify = ( ( invert1 ? astGetSimpIF( mathmap1 ) :
                               astGetSimpFI( mathmap1 ) ) &&
                   ( invert2 ? astGetSimpFI( mathmap2 ) :
                               astGetSimpIF( mathmap2 ) ) );
   }

/* If still OK, obtain the effective numbers of input coordinates for
   the first MathMap and output coordinates for the second. Take account
   of the associated invert flags and the way the Invert attribute of
   each MathMap is currently set. */
   if ( astOK && simplify ) {
      nin1 = ( invert1 == astGetInvert( mathmap1 ) ) ?
             astGetNin( mathmap1 ) : astGetNout( mathmap1 );
      nout2 = ( invert2 == astGetInvert( mathmap2 ) ) ?
              astGetNout( mathmap2 ) : astGetNin( mathmap2 );

/* Simplification is only possible if these two numbers are equal
   (otherwise the the two MathMaps cannot be identical). */
      simplify = ( nin1 == nout2 );
   }

/* If still OK, obtain the effective number of forward transformation
   functions for the first MathMap (allowing for the associated invert
   flag). Similarly, obtain the effective number of inverse
   transformation functions for the second MathMap. */
   if ( astOK && simplify ) {
      nfwd1 = !invert1 ? mathmap1->nfwd : mathmap1->ninv;
      ninv2 = !invert2 ? mathmap2->ninv : mathmap2->nfwd;

/* Check whether these values are equal. The MathMaps cannot be
   identical if they are not. */
      simplify = ( nfwd1 == ninv2 );
   }

/* As above, obtain pointers to the array of effective forward
   transformation functions for the first MathMap, and the effective
   inverse transformation functions for the second MathMap. */
   if ( astOK && simplify ) {
      fwd1 = !invert1 ? mathmap1->fwdfun : mathmap1->invfun;
      inv2 = !invert2 ? mathmap2->invfun : mathmap2->fwdfun;

/* Loop to check whether these two sets of functions are
   identical. The MathMaps cannot be merged unless they are. */
      for ( ifun = 0; ifun < nfwd1; ifun++ ) {
         simplify = !strcmp( fwd1[ ifun ], inv2[ ifun ] );
         if ( !simplify ) break;
      }
   }

/* If OK, repeat the above process to compare the effective inverse
   transformation functions of the first MathMap with the forward
   functions of the second one. */
   if ( astOK && simplify ) {
      ninv1 = !invert1 ? mathmap1->ninv : mathmap1->nfwd;
      nfwd2 = !invert2 ? mathmap2->nfwd : mathmap2->ninv;
      simplify = ( ninv1 == nfwd2 );
   }
   if ( astOK && simplify ) {
      inv1 = !invert1 ? mathmap1->invfun : mathmap1->fwdfun;
      fwd2 = !invert2 ? mathmap2->fwdfun : mathmap2->invfun;
      for ( ifun = 0; ifun < ninv1; ifun++ ) {
         simplify = !strcmp( inv1[ ifun ], fwd2[ ifun ] );
         if ( !simplify ) break;
      }
   }

/* If the two MathMaps can be merged, create a UnitMap as a
   replacement. */
   if ( astOK && simplify ) {
      new = (AstMapping *) astUnitMap( nin1, "", status );

/* If OK, annul the pointers to the original MathMaps. */
      if ( astOK ) {
         ( *map_list )[ imap1 ] = astAnnul( ( *map_list )[ imap1 ] );
         ( *map_list )[ imap2 ] = astAnnul( ( *map_list )[ imap2 ] );

/* Insert the pointer to the replacement UnitMap and store the
   associated invert flag. */
         ( *map_list )[ imap1 ] = new;
         ( *invert_list )[ imap1 ] = 0;

/* Loop to move the following Mapping pointers and invert flags down
   in their arrays to close the gap. */
         for ( imap = imap2 + 1; imap < *nmap; imap++ ) {
            ( *map_list )[ imap - 1 ] = ( *map_list )[ imap ];
            ( *invert_list )[ imap - 1 ] = ( *invert_list )[ imap ];
         }

/* Clear the final entry in each array. */
         ( *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 = imap1;
      }
   }

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

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

static void ParseConstant( const char *method, const char *class,
                           const char *exprs, int istart, int *iend,
                           double *con, int *status ) {
/*
*  Name:
*     ParseConstant

*  Purpose:
*     Parse a constant.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     void ParseConstant( const char *method, const char *class,
*                         const char *exprs, int istart, int *iend,
*                         double *con, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     This routine parses an expression, looking for a constant starting at
*     the character with index "istart" in the string "exprs". If it
*     identifies the constant successfully, "*con" it will return its value
*     and "*iend" will be set to the index of the final constant character
*     in "exprs".
*
*     If the characters encountered are clearly not part of a constant (it
*     does not begin with a numeral or decimal point) the function returns
*     with "*con" set to zero and "*iend" set to -1, but without reporting
*     an error. However, if the first character appears to be a constant but
*     its syntax proves to be invalid, then an error is reported.
*
*     The expression must be in lower case with no embedded white space.
*     The constant must not have a sign (+ or -) in front of it.

*  Parameters:
*     method
*        Pointer to a constant null-terminated character string
*        containing the name of the method that invoked this function.
*        This method name is used solely for constructing error messages.
*     class
*        Pointer to a constant null-terminated character string containing the
*        class name of the Object being processed. This name is used solely
*        for constructing error messages.
*     exprs
*        Pointer to a null-terminated string containing the expression
*        to be parsed.
*     istart
*        Index of the first character in "exprs" to be considered by this
*        function.
*     iend
*        Pointer to an int in which to return the index in "exprs" of the
*        final character which forms part of the constant. If no constant is
*        found, a value of -1 is returned.
*     con
*        Pointer to a double, in which the value of the constant, if found,
*        will be returned.
*     status
*        Pointer to the inherited status variable.
*/

/* Local Variables: */
   char *str;                    /* Pointer to temporary string */
   char c;                       /* Single character from the expression */
   int dpoint;                   /* Decimal point encountered? */
   int expon;                    /* Exponent character encountered? */
   int i;                        /* Loop counter for characters */
   int iscon;		         /* Character is part of the constant? */
   int n;                        /* Number of values read by astSscanf */
   int nc;                       /* Number of characters read by astSscanf */
   int numer;                    /* Numeral encountered in current field? */
   int sign;                     /* Sign encountered? */
   int valid;		         /* Constant syntax valid? */

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

/* Initialise. */
   *con = 0.0;
   *iend = -1;

/* Check if the expression starts with a numeral or a decimal point. */
   c = exprs[ istart ];
   numer = isdigit( c );
   dpoint = ( c == '.' );

/* If it begins with any of these, the expression is clearly intended
   to be a constant, so any failure beyond this point will result in an
   error. Otherwise, failure to find a constant is not an error. */
   if ( numer || dpoint ) {

/* Initialise remaining variables specifying the parser context. */
      expon = 0;
      sign = 0;
      valid = 1;

/* Loop to increment the last constant character position until the
   following character in the expression does not look like part of the
   constant. */
      *iend = istart;
      iscon = 1;
      while ( ( c = exprs[ *iend + 1 ] ) && iscon ) {
         iscon = 0;

/* It may be part of a numerical constant if it is a numeral, wherever
   it occurs. */
         if ( isdigit( c ) ) {
            numer = 1;
            iscon = 1;

/* Or a decimal point, so long as it is the first one and is not in
   the exponent field. Otherwise it is invalid. */
         } else if ( c == '.' ) {
            if ( !( dpoint || expon ) ) {
               dpoint = 1;
               iscon = 1;
            } else {
               valid = 0;
            }

/* Or if it is a 'd' or 'e' exponent character, so long as it is the
   first one and at least one numeral has been encountered first.
   Otherwise it is invalid. */
          } else if ( ( c == 'd' ) || ( c == 'e' ) ) {
             if ( !expon && numer ) {
                expon = 1;
                numer = 0;
                iscon = 1;
             } else {
                valid = 0;
             }

/* Or if it is a sign, so long as it is in the exponent field and is
   the first sign with no previous numerals in the same field. Otherwise
   it is invalid (unless numerals have been encountered, in which case it
   marks the end of the constant). */
          } else if ( ( c == '+' ) || ( c == '-' ) ) {
             if ( expon && !sign && !numer ) {
                sign = 1;
                iscon = 1;
             } else if ( !numer ) {
                valid = 0;
             }
          }

/* Increment the character count if the next character may be part of
   the constant, or if it was invalid (it will then form part of the
   error message). */
          if ( iscon || !valid ) ( *iend )++;
      }

/* Finally, check that the last field contained a numeral. */
      valid = ( valid && numer );

/* If the constant appears valid, allocate a temporary string to hold
   it. */
      if ( valid ) {
         str = astMalloc( (size_t) ( *iend - istart + 2 ) );
         if ( astOK ) {

/* Copy the constant's characters, changing 'd' to 'e' so that
   "astSscanf" will recognise it as an exponent character. */
            for ( i = istart; i <= *iend; i++ ) {
               str[ i - istart ] = ( exprs[ i ] == 'd' ) ? 'e' : exprs[ i ];
            }
            str[ *iend - istart + 1 ] = '\0';

/* Attempt to read the constant as a double, noting how many values
   are read and how many characters consumed. */
            n = astSscanf( str, "%lf%n", con, &nc );

/* Check that one value was read and all the characters consumed. If
   not, then the constant's syntax is invalid. */
            if ( ( n != 1 ) || ( nc < ( *iend - istart + 1 ) ) ) valid = 0;
         }

/* Free the temporary string. */
         str = astFree( str );
      }

/* If the constant syntax is invalid, and no other error has occurred,
   then report an error. */
      if ( astOK && !valid ) {
         astError( AST__CONIN,
                   "%s(%s): Invalid constant syntax in the expression "
                   "\"%.*s\".", status,
                   method, class, *iend + 1, exprs );
      }

/* If an error occurred, reset the output values. */
      if ( !astOK ) {
         *iend = -1;
         *con = 0.0;
      }
   }
}

static void ParseName( const char *exprs, int istart, int *iend, int *status ) {
/*
*  Name:
*     ParseName

*  Purpose:
*     Parse a name.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     void ParseName( const char *exprs, int istart, int *iend, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     This routine parses an expression, looking for a name starting at the
*     character with index "istart" in the string "exprs". If it identifies
*     a name successfully, "*iend" will return the index of the final name
*     character in "exprs". A name must begin with an alphabetic character
*     and subsequently contain only alphanumeric characters or underscores.
*
*     If the expression does not contain a name at the specified location,
*     "*iend" is set to -1. No error results.
*
*     The expression should not contain embedded white space.

*  Parameters:
*     exprs
*        Pointer to a null-terminated string containing the expression
*        to be parsed.
*     istart
*        Index of the first character in "exprs" to be considered by this
*        function.
*     iend
*        Pointer to an int in which to return the index in "exprs" of the
*        final character which forms part of the name. If no name is
*        found, a value of -1 is returned.
*     status
*        Pointer to the inherited status variable.
*/

/* Local Variables: */
   char c;                       /* Single character from expression */

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

/* Initialise. */
   *iend = -1;

/* Check the first character is valid for a name (alphabetic). */
   if ( isalpha( exprs[ istart ] ) ) {

/* If so, loop to inspect each subsequent character until one is found
   which is not part of a name (not alphanumeric or underscore). */
      for ( *iend = istart; ( c = exprs[ *iend + 1 ] ); ( *iend )++ ) {
         if ( !( isalnum( c ) || ( c == '_' ) ) ) break;
      }
   }
}

static void ParseVariable( const char *method, const char *class,
                           const char *exprs, int istart, int nvar,
                           const char *var[], int *ivar, int *iend, int *status ) {
/*
*  Name:
*     ParseVariable

*  Purpose:
*     Parse a variable name.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     void ParseVariable( const char *method, const char *class,
*                         const char *exprs, int istart, int nvar,
*                         const char *var[], int *ivar, int *iend, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     This routine parses an expression, looking for a recognised variable
*     name starting at the character with index "istart" in the string
*     "exprs". If it identifies a variable name successfully, "*ivar" will
*     return a value identifying it and "*iend" will return the index of the
*     final variable name character in "exprs". To be recognised, a name
*     must begin with an alphabetic character and subsequently contain only
*     alphanumeric characters or underscores. It must also appear in the
*     list of defined variable names supplied to this function.
*
*     If the expression does not contain a name at the specified location,
*     "*ivar" and "*iend" are set to -1 and no error results. However, if
*     the expression contains a name but it is not in the list of defined
*     variable names supplied, then an error is reported.
*
*     This function is case sensitive. The expression should not contain
*     embedded white space.

*  Parameters:
*     method
*        Pointer to a constant null-terminated character string
*        containing the name of the method that invoked this function.
*        This method name is used solely for constructing error messages.
*     class
*        Pointer to a constant null-terminated character string containing the
*        class name of the Object being processed. This name is used solely
*        for constructing error messages.
*     exprs
*        Pointer to a null-terminated string containing the expression
*        to be parsed.
*     istart
*        Index of the first character in "exprs" to be considered by this
*        function.
*     nvar
*        The number of defined variable names.
*     var
*        An array of pointers (with "nvar" elements) to null-terminated
*        strings. Each of these should contain a variable name to be
*        recognised. These strings are case sensitive and should contain
*        no white space.
*     ivar
*        Pointer to an int in which to return the index in "vars" of the
*        variable name found. If no variable name is found, a value of -1
*        is returned.
*     iend
*        Pointer to an int in which to return the index in "exprs" of the
*        final character which forms part of the variable name. If no variable
*        name is found, a value of -1 is returned.
*     status
*        Pointer to the inherited status variable.
*/

/* Local Variables: */
   int found;                    /* Variable name recognised? */
   int nc;                       /* Number of characters in variable name */

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

/* Initialise. */
   *ivar = -1;
   *iend = -1;

/* Determine if the characters in the expression starting at index
   "istart" constitute a valid name. */
   ParseName( exprs, istart, iend, status );

/* If so, calculate the length of the name. */
   if ( *iend >= istart ) {
      nc = *iend - istart + 1;

/* Loop to compare the name with the list of variable names
   supplied. */
      found = 0;
      for ( *ivar = 0; *ivar < nvar; ( *ivar )++ ) {
         found = ( nc == (int) strlen( var[ *ivar ] ) ) &&
                 !strncmp( exprs + istart, var[ *ivar ], (size_t) nc );

/* Break if the name is recognised. */
         if ( found ) break;
      }

/* If it was not recognised, then report an error and reset the output
   values. */
      if ( !found ) {
         astError( AST__UDVOF,
                   "%s(%s): Undefined variable or function in the expression "
                   "\"%.*s\".", status,
                   method, class, *iend + 1, exprs );
         *ivar = -1;
         *iend = -1;
      }
   }
}

static double Poisson( Rcontext *context, double mean, int *status ) {
/*
*  Name:
*     Poisson

*  Purpose:
*     Produce a pseudo-random sample from a Poisson distribution.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     double Poisson( Rcontext *context, double mean, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     On each invocation, this function returns a pseudo-random sample drawn
*     from a Poisson distribution with a specified mean. A combination of
*     methods is used, depending on the value of the mean. The algorithm is
*     based on that given by Press et al. (Numerical Recipes), but
*     re-implemented and extended.

*  Parameters:
*     context
*        Pointer to an Rcontext structure which holds the random number
*        generator's context between invocations.
*     mean
*        The mean of the Poisson distribution, which should not be
*        negative.
*     status
*        Pointer to the inherited status variable.

*  Returned Value:
*     A sample (which will only take integer values) from the Poisson
*     distribution, or AST__BAD if the mean supplied is negative.

*  Notes:
*     - The sequence of numbers returned is determined by the "seed"
*     value in the Rcontext structure supplied.
*     - If the seed value is changed, the "active" flag must also be cleared
*     so that this function can re-initiallise the Rcontext structure before
*     generating the next pseudo-random number. The "active" flag should
*     also be clear to force initialisation the first time an Rcontext
*     structure is used.
*     - This function does not perform error checking and does not generate
*     errors. It will execute even if the global error status is set.
*/

/* Local Constants: */
   const double small = 9.3;     /* "Small" distribution mean value */

/* Local Variables: */
   double pfract;                /* Probability of accepting sample */
   double product;               /* Product of random samples */
   double ran;                   /* Sample from Lorentzian distribution */
   double result;                /* Result value to return */
   static double beta;           /* Constant for forming acceptance ratio */
   static double huge;           /* Large mean where std. dev. is negligible */
   static double last_mean;      /* Value of "mean" on last invocation */
   static double log_mean;       /* Logarithm of "mean" */
   static double pi;             /* Value of pi */
   static double ranmax;         /* Maximum safe value of "ran" */
   static double root_2mean;     /* sqrt( 2.0 * mean ) */
   static double sqrt_point9;    /* Square root of 0.9 */
   static double thresh;         /* Threshold for product of samples */
   static int init = 0;          /* Local initialisation performed? */

   LOCK_MUTEX6

/* If initialisation has not yet been performed, then perform it
   now. */
   if ( !init ) {

/* Initialise the mean value from the previous invocation. */
      last_mean = -1.0;

/* Calculate simple constants. */
      pi = acos( -1.0 );
      sqrt_point9 = sqrt( 0.9 );

/* Calculate the value of the distribution mean for which the smallest
   representable deviation from the mean permitted by the machine
   precision is one thousand standard deviations. */
      huge = pow( 1.0e3 / DBL_EPSILON, 2.0 );

/* Calculate the largest value such that
   (0.9+(sqrt_point9*ranmax)*(sqrt_point9*ranmax)) doesn't overflow,
   allowing a small margin for rounding error. */
      ranmax = ( sqrt( DBL_MAX - 0.9 ) / sqrt( 0.9 ) ) *
               ( 1.0 - 4.0 * DBL_EPSILON );

/* Note that initialisation has been performed. */
      init = 1;
   }

/* If the distribution mean is less than zero, then return a bad
   result. */
   if ( mean < 0.0 ) {
      result = AST__BAD;

/* If the mean is zero, then the result can only be zero. */
   } else if ( mean == 0.0 ) {
      result = 0.0;

/* Otherwise, if the mean is sufficiently small, we can use the direct
   method of summing a series of exponentially distributed random samples
   and counting the number which occur before the mean is exceeded. This
   is equivalent to multiplying a series of uniformly distributed
   samples and counting the number which occur before the product
   becomes less then an equivalent threshold. */
   } else if ( mean <= small ) {

/* If the mean has changed since the last invocation, store the new
   mean and calculate a new threshold. */
      if ( mean != last_mean ) {
         last_mean = mean;
         thresh = exp( -mean );
      }

/* Initialise the product and the result. */
      product = 1.0;
      result = -1.0;

/* Multiply the random samples, counting the number needed to reach
   the threshold. */
      do {
         product *= Rand( context, status );
         result += 1.0;
      } while ( product > thresh );

/* Otherwise, if the distribution mean is large (but not huge), we
   must use an indirect rejection method. */
   } else if ( mean <= huge ) {

/* If the mean has changed since the last invocation, then
   re-calculate the constants required below. Note that because of the
   restrictions we have placed on "mean", these calculations are safe
   against overflow. */
      if ( mean != last_mean ) {
         last_mean = mean;
         log_mean = log( mean );
         root_2mean = sqrt( 2.0 * mean );
         beta = mean * log_mean - LogGamma( mean + 1.0, status );
      }

/* Loop until a suitable random sample has been generated. */
      do {
         do {

/* First transform a sample from a uniform distribution to obtain a
   sample from a Lorentzian distribution. Check that the result is not so
   large as to cause overflow later. Also check for overflow in the maths
   library. If necessary, obtain a new sample. */
            do {
               errno = 0;
               ran = tan( pi * Rand( context, status ) );
            } while ( ( ran > ranmax ) ||
                      ( ( errno == ERANGE ) &&
                        ( ( ( ran >= 0.0 ) ? ran : -ran ) == HUGE_VAL ) ) );

/* If OK, scale the sample and add a constant so that the sample's
   distribution approximates the Poisson distribution we
   require. Overflow is prevented by the check on "ran" above, together
   with the restricted value of "mean". */
            result = ran * root_2mean + mean;

/* If the result is less than zero (where the Poisson distribution has
   value zero), then obtain a new sample. */
         } while ( result < 0.0 );

/* Round down to an integer, so that the sample is valid for a Poisson
   distribution. */
         result = floor( result );

/* Calculate the ratio between the required Poisson distribution and
   the Lorentzian from which we have sampled (the factor of 0.9 prevents
   this exceeding 1.0, and overflow is again prevented by the checks
   performed above). */
         ran *= sqrt_point9;
         pfract = ( 0.9 + ran * ran ) *
                  exp( result * log_mean - LogGamma( result + 1.0, status ) - beta );

/* Accept the sample with this fractional probability, otherwise
   obtain a new sample. */
      } while ( Rand( context, status ) > pfract );

/* If the mean is huge, the relative standard deviation will be
   negligible compared to the machine precision. In such cases, the
   probability of getting a result that differs from the mean is
   effectively zero, so we can simply return the mean. */
   } else {
      result = mean;
   }

   UNLOCK_MUTEX6

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

static double Rand( Rcontext *context, int *status ) {
/*
*  Name:
*     Rand

*  Purpose:
*     Produce a uniformly distributed pseudo-random number.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     double Rand( Rcontext *context, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     On each invocation, this function returns a pseudo-random number
*     uniformly distributed in the range 0.0 to 1.0 (inclusive). The
*     underlying algorithm is that used by the "ran2" function of Press et
*     al. (Numerical Recipes), which has a long period and good statistical
*     properties. This independent implementation returns double precision
*     values.

*  Parameters:
*     context
*        Pointer to an Rcontext structure which holds the random number
*        generator's context between invocations.
*     status
*        Pointer to the inherited status variable.

*  Notes:
*     - The sequence of numbers returned is determined by the "seed"
*     value in the Rcontext structure supplied.
*     - If the seed value is changed, the "active" flag must also be cleared
*     so that this function can re-initiallise the Rcontext structure before
*     generating the next pseudo-random number. The "active" flag should
*     also be clear to force initialisation the first time an Rcontext
*     structure is used.
*     - This function does not perform error checking and does not generate
*     errors. It will execute even if the global error status is set.
*/

/* Local Constants: */
   const long int a1 = 40014L;   /* Random number generator constants... */
   const long int a2 = 40692L;
   const long int m1 = 2147483563L;
   const long int m2 = 2147483399L;
   const long int q1 = 53668L;
   const long int q2 = 52774L;
   const long int r1 = 12211L;
   const long int r2 = 3791L;
   const int ntab =              /* Size of shuffle table */
      AST_MATHMAP_RAND_CONTEXT_NTAB_;
   const int nwarm = 8;          /* Number of warm-up iterations */

/* Local Variables: */
   double result;                /* Result value to return */
   double scale;                 /* Scale factor for random integers */
   double sum;                   /* Sum for forming normalisation constant */
   int dbits;                    /* Approximate bits in double mantissa */
   int irand;                    /* Loop counter for random integers */
   int itab;                     /* Loop counter for shuffle table */
   int lbits;                    /* Approximate bits used by generators */
   long int seed;                /* Random number seed */
   long int tmp;                 /* Temporary variable */
   static double norm;           /* Normalisation constant */
   static double scale0;         /* Scale decrement for successive integers */
   static int init = 0;          /* Local initialisation performed? */
   static int nrand;             /* Number of random integers to use */

/* If the random number generator context is not active, then
   initialise it. */
   if ( !context->active ) {

/* First, perform local initialisation for this function, if not
   already done. */
      LOCK_MUTEX4
      if ( !init ) {

/* Obtain the approximate number of bits used by the random integer
   generator from the value "m1". */
         (void) frexp( (double) m1, &lbits );

/* Obtain the approximate number of bits used by the mantissa of the
   double value we want to produce, allowing for the (unlikely)
   possibility that the mantissa's radix isn't 2. */
         dbits = (int) ceil( (double) DBL_MANT_DIG *
                             log( (double) FLT_RADIX ) / log( 2.0 ) );

/* Hence determine how many random integers we need to combine to
   produce each double value, so that all the mantissa's bits will be
   used. */
         nrand = ( dbits + lbits - 1 ) / lbits;

/* Calculate the scale factor by which each successive random
   integer's contribution to the result is reduced so as to generate
   progressively less significant bits. */
         scale0 = 1.0 / (double) ( m1 - 1L );

/* Loop to sum the maximum contributions from each random integer
   (assuming that each takes the largest possible value, of "m1-1",
   from which we will later subtract 1). This produces the normalisation
   factor by which the result must be scaled so as to lie between 0.0 and
   1.0 (inclusive). */
         sum = 0.0;
         scale = 1.0;
         for ( irand = 0; irand < nrand; irand++ ) {
            scale *= scale0;
            sum += scale;
         }
         norm = 1.0 / ( sum * (double) ( m1 - 2L ) );

/* Note that local initialisation has been done. */
         init = 1;
      }
      UNLOCK_MUTEX4

/* Obtain the seed value, enforcing positivity. */
      seed = (long int) context->seed;
      if ( seed < 1 ) seed = seed + LONG_MAX;
      if ( seed < 1 ) seed = LONG_MAX;

/* Initialise the random number generators with this seed. */
      context->rand1 = context->rand2 = seed;

/* Now loop to initialise the shuffle table with an initial set of
   random values. We generate more values than required in order to "warm
   up" the generator before recording values in the table. */
      for ( itab = ntab + nwarm - 1; itab >= 0; itab-- ) {

/* Repeatedly update "rand1" from the expression "(rand1*a1)%m1" while
   avoiding overflow. */
         tmp = context->rand1 / q1;
         context->rand1 = a1 * ( context->rand1 - tmp * q1 ) - tmp * r1;
         if ( context->rand1 < 0L ) context->rand1 += m1;

/* After warming up, start recording values in the table. */
         if ( itab < ntab ) context->table[ itab ] = context->rand1;
      }

/* Record the last entry in the table as the "previous" random
   integer. */
      context->random_int = context->table[ 0 ];

/* Note the random number generator context is active. */
      context->active = 1;
   }

/* Generate a random value. */
/* ------------------------ */
/* Initialise. */
   result = 0.0;

/* Loop to generate sufficient random integers to combine into a
   double value. */
   scale = norm;
   for ( irand = 0; irand < nrand; irand++ ) {

/* Update the first generator "rand1" from the expression
   "(a1*rand1)%m1" while avoiding overflow. */
      tmp = context->rand1 / q1;
      context->rand1 = a1 * ( context->rand1 - tmp * q1 ) - tmp * r1;
      if ( context->rand1 < 0L ) context->rand1 += m1;

/* Similarly, update the second generator "rand2" from the expression
   "(a2*rand2)%m2". */
      tmp = context->rand2 / q2;
      context->rand2 = a2 * ( context->rand2 - tmp * q2 ) - tmp * r2;
      if ( context->rand2 < 0L ) context->rand2 += m2;

/* Use the previous random integer to generate an index into the
   shuffle table. */
      itab = (int) ( context->random_int /
                     ( 1L + ( m1 - 1L ) / (long int) ntab ) );

/* The algorithm left by RFWS seems to have a bug that "itab" can
   sometimes be outside the range of [0.,ntab-1] causing the context->table
   array to be addressed out of bounds. To avoid this, use the
   following sticking plaster, since I'm not sure what the correct fix is. */
      if( itab < 0 ) itab = -itab;
      itab = itab % ntab;

/* Extract the table entry and replace it with a new random value from
   the first generator "rand1". This is the Bays-Durham shuffle. */
      context->random_int = context->table[ itab ];
      context->table[ itab ] = context->rand1;

/* Combine the extracted value with the latest value from the second
   generator "rand2". */
      context->random_int -= context->rand2;
      if ( context->random_int < 1L ) context->random_int += m1 - 1L;

/* Update the scale factor to apply to the resulting random integer
   and accumulate its contribution to the result. */
      scale *= scale0;
      result += scale * (double) ( context->random_int - 1L );
   }

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

*  Type:
*     Private function.

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

*  Class Membership:
*     MathMap member function (extends the astSetAttrib method inherited from
*     the Mapping class).

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

*  Returned Value:
*     void
*/

/* Local Vaiables: */
   AstMathMap *this;             /* Pointer to the MathMap structure */
   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 MathMap structure. */
   this = (AstMathMap *) this_object;

/* Obtain the length of the setting string. */
   len = 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. */

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

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

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

/* Pass any unrecognised setting 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 MathMap.

*  Type:
*     Private function.

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

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

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

/* Initialise. */
   result = 0;

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

/* Obtain a pointer to the MathMap structure. */
   this = (AstMathMap *) this_object;

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

/* Seed. */
/* ----- */
   if ( !strcmp( attrib, "seed" ) ) {
      result = astTestSeed( this );

/* SimpFI. */
/* ------- */
   } else if ( !strcmp( attrib, "simpfi" ) ) {
      result = astTestSimpFI( this );

/* SimpIF. */
/* ------- */
   } else if ( !strcmp( attrib, "simpif" ) ) {
      result = astTestSimpIF( this );

/* If the attribute is 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 *map, AstPointSet *in,
                               int forward, AstPointSet *out, int *status ) {
/*
*  Name:
*     Transform

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

*  Type:
*     Private function.

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

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

*  Description:
*     This function takes a MathMap and a set of points encapsulated in a
*     PointSet and transforms the points so as to apply the required coordinate
*     transformation.

*  Parameters:
*     map
*        Pointer to the MathMap.
*     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 coordinates for the MathMap 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: */
   AstMathMap *this;             /* Pointer to MathMap to be applied */
   AstPointSet *result;          /* Pointer to output PointSet */
   double **data_ptr;            /* Array of pointers to coordinate data */
   double **ptr_in;              /* Pointer to input coordinate data */
   double **ptr_out;             /* Pointer to output coordinate data */
   double *work;                 /* Workspace for intermediate results */
   int idata;                    /* Loop counter for data pointer elements */
   int ifun;                     /* Loop counter for functions */
   int ncoord_in;                /* Number of coordinates per input point */
   int ncoord_out;               /* Number of coordinates per output point */
   int ndata;                    /* Number of data pointer elements filled */
   int nfun;                     /* Number of functions to evaluate */
   int npoint;                   /* Number of points */

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

/* Initialise variables to avoid "used of uninitialised variable"
   messages from dumb compilers. */
   work = NULL;

/* Obtain a pointer to the MathMap. */
   this = (AstMathMap *) map;

/* 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)( map, in, forward, out, status );

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

/* 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. */
   ncoord_in = astGetNcoord( in );
   ncoord_out = astGetNcoord( result );
   npoint = astGetNpoint( in );
   ptr_in = astGetPoints( in );
   ptr_out = astGetPoints( result );

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

/* Obtain the number of transformation functions that must be
   evaluated to perform the transformation. This will include any that
   produce intermediate results from which the final results are
   calculated. */
   nfun = forward ? this->nfwd : this->ninv;

/* If intermediate results are to be calculated, then allocate
   workspace to hold them (each intermediate result being a vector of
   "npoint" double values). */
   if ( nfun > ncoord_out ) {
      work = astMalloc( sizeof( double) *
                        (size_t) ( npoint * ( nfun - ncoord_out ) ) );
   }

/* Also allocate space for an array to hold pointers to the input
   data, intermediate results and output data. */
   data_ptr = astMalloc( sizeof( double * ) * (size_t) ( ncoord_in + nfun ) );

/* We now set up the "data_ptr" array to locate the data to be
   processed. */
   if ( astOK ) {

/* The first elements of this array point at the input data
   vectors. */
      ndata = 0;
      for ( idata = 0; idata < ncoord_in; idata++ ) {
         data_ptr[ ndata++ ] = ptr_in[ idata ];
      }

/* The following elements point at successive vectors within the
   workspace array (if allocated). These vectors will act first as output
   arrays for intermediate results, and then as input arrays for
   subsequent calculations which use these results. */
      for ( idata = 0; idata < ( nfun - ncoord_out ); idata++ ) {
         data_ptr[ ndata++ ] = work + ( idata * npoint );
      }

/* The final elements point at the output coordinate data arrays into
   which the final results will be written. */
      for ( idata = 0; idata < ncoord_out; idata++ ) {
         data_ptr[ ndata++ ] = ptr_out[ idata ];
      }

/* Perform coordinate transformation. */
/* ---------------------------------- */
/* Loop to evaluate each transformation function in turn. */
      for ( ifun = 0; ifun < nfun; ifun++ ) {

/* Invoke the function that evaluates compiled expressions. Pass the
   appropriate code and constants arrays, depending on the direction of
   coordinate transformation, together with the required stack size. The
   output array is the vector located by successive elements of the
   "data_ptr" array (skipping the input data elements), while the
   function has access to all previous elements of the "data_ptr" array
   to locate the required input data. */
         EvaluateFunction( &this->rcontext, npoint, (const double **) data_ptr,
                           forward ? this->fwdcode[ ifun ] :
                                     this->invcode[ ifun ],
                           forward ? this->fwdcon[ ifun ] :
                                     this->invcon[ ifun ],
                           forward ? this->fwdstack : this->invstack,
                           data_ptr[ ifun + ncoord_in ], status );
      }
   }

/* Free the array of data pointers and any workspace allocated for
   intermediate results. */
   data_ptr = astFree( data_ptr );
   if ( nfun > ncoord_out ) work = astFree( work );

/* If an error occurred, then return a NULL pointer. If no output
   PointSet was supplied, also delete any new one that may have been
   created. */
   if ( !astOK ) {
      result = ( result == out ) ? NULL : astDelete( result );
   }

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

static void ValidateSymbol( const char *method, const char *class,
                            const char *exprs, int iend, int sym,
                            int *lpar, int **argcount, int **opensym,
                            int *ncon, double **con, int *status ) {
/*
*  Name:
*     ValidateSymbol

*  Purpose:
*     Validate a symbol in an expression.

*  Type:
*     Private function.

*  Synopsis:
*     #include "mathmap.h"
*     void ValidateSymbol( const char *method, const char *class,
*                          const char *exprs, int iend, int sym, int *lpar,
*                          int **argcount, int **opensym, int *ncon,
*                          double **con, int *status )

*  Class Membership:
*     MathMap member function.

*  Description:
*     This function validates an identified standard symbol during
*     compilation of an expression. Its main task is to keep track of the
*     level of parenthesis in the expression and to count the number of
*     arguments supplied to functions at each level of parenthesis (for
*     nested function calls). On this basis it is able to interpret and
*     accept or reject symbols which represent function calls, parentheses
*     and delimiters. Other symbols are accepted automatically.

*  Parameters:
*     method
*        Pointer to a constant null-terminated character string
*        containing the name of the method that invoked this function.
*        This method name is used solely for constructing error messages.
*     class
*        Pointer to a constant null-terminated character string containing the
*        class name of the Object being processed. This name is used solely
*        for constructing error messages.
*     exprs
*        Pointer to a null-terminated string containing the expression
*        being parsed. This is only used for constructing error messages.
*     iend
*        Index in "exprs" of the last character belonging to the most
*        recently identified symbol. This is only used for constructing error
*        messages.
*     sym
*        Index in the static "symbol" array of the most recently identified
*        symbol in the expression. This is the symbol to be verified.
*     lpar
*        Pointer to an int which holds the current level of parenthesis. On
*        the first invocation, this should be zero. The returned value should
*        be passed to subsequent invocations.
*     argcount
*        Address of a pointer to a dynamically allocated array of int in
*        which argument count information is maintained for each level of
*        parenthesis (e.g. for nested function calls). On the first invocation,
*        "*argcount" should be NULL. This function will allocate the required
*        space as needed and update this pointer. The returned pointer value
*        should be passed to subsequent invocations.
*
*        The allocated space must be freed by the caller (using astFree) when
*        no longer required.
*     opensym
*        Address of a pointer to a dynamically allocated array of int, in which
*        information is maintained about the functions associated with each
*        level of parenthesis (e.g. for nested function calls). On the first
*        invocation, "*opensym" should be NULL. This function will allocate the
*        required space as needed and update this pointer. The returned pointer
*        value should be passed to subsequent invocations.
*
*        The allocated space must be freed by the caller (using astFree) when
*        no longer required.
*     ncon
*        Pointer to an int which holds a count of the constants associated
*        with the expression (and determines the size of the "*con" array).
*        This function will update the count to reflect any new constants
*        appended to the "*con" array and the returned value should be passed
*        to subsequent invocations.
*     con
*        Address of a pointer to a dynamically allocated array of double, in
*        which the constants associated with the expression being parsed are
*        accumulated. On entry, "*con" should point at a dynamic array with
*        at least "*ncon" elements containing existing constants (or may be
*        NULL if no constants have yet been stored). This function will
*        allocate the required space as needed and update this pointer (and
*        "*ncon") appropriately. The returned pointer value should be passed
*        to subsequent invocations.
*
*        The allocated space must be freed by the caller (using astFree) when
*        no longer required.
*     status
*        Pointer to the inherited status variable.

*  Notes:
*     - The dynamically allocated arrays normally returned by this function
*     will be freed and NULL pointers 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, but do not return at this point
   because dynamic arrays may require freeing. */
   if ( astOK ) {

/* Check if the symbol is a comma. */
      if ( ( symbol[ sym ].text[ 0 ] == ',' ) &&
           ( symbol[ sym ].text[ 1 ] == '\0' ) ) {

/* A comma is only used to delimit function arguments. If the current
   level of parenthesis is zero, or the symbol which opened the current
   level of parenthesis was not a function call (indicated by an argument
   count of zero at the current level of parenthesis), then report an
   error. */
         if ( ( *lpar <= 0 ) || ( ( *argcount )[ *lpar - 1 ] == 0 ) ) {
            astError( AST__COMIN,
                      "%s(%s): Spurious comma encountered in the expression "
                      "\"%.*s\".", status,
                      method, class, iend + 1, exprs );

/* If a comma is valid, then increment the argument count at the
   current level of parenthesis. */
         } else {
            ( *argcount )[ *lpar - 1 ]++;
         }

/* If the symbol is not a comma, check if it increases the current
   level of parenthesis. */
      } else if ( symbol[ sym ].parincrement > 0 ) {

/* Increase the size of the arrays which hold parenthesis level
   information and check for errors. */
         *argcount = astGrow( *argcount, *lpar + 1, sizeof( int ) );
         *opensym = astGrow( *opensym, *lpar + 1, sizeof( int ) );
         if ( astOK ) {

/* Increment the level of parenthesis and initialise the argument
   count at the new level. This count is set to zero if the symbol which
   opens the parenthesis level is not a function call (indicated by a
   zero "nargs" entry in the symbol data), and it subsequently remains at
   zero. If the symbol is a function call, the argument count is
   initially set to 1 and increments whenever a comma is encountered at
   this parenthesis level. */
            ( *argcount )[ ++( *lpar ) - 1 ] = ( symbol[ sym ].nargs != 0 );

/* Remember the symbol which opened this parenthesis level. */
            ( *opensym )[ *lpar - 1 ] = sym;
         }

/* Check if the symbol decreases the current parenthesis level. */
      } else if ( symbol[ sym ].parincrement < 0 ) {

/* Ensure that the parenthesis level is not already at zero. If it is,
   then there is a missing left parenthesis in the expression being
   compiled, so report an error. */
         if ( *lpar == 0 ) {
            astError( AST__MLPAR,
                      "%s(%s): Missing left parenthesis in the expression "
                      "\"%.*s\".", status,
                      method, class, iend + 1, exprs );

/* If the parenthesis level is valid and the symbol which opened this
   level of parenthesis was a function call with a fixed number of
   arguments (indicated by a positive "nargs" entry in the symbol data),
   then we must check the number of function arguments which have been
   encountered. */
         } else if ( symbol[ ( *opensym )[ *lpar - 1 ] ].nargs > 0 ) {

/* Report an error if the number of arguments is wrong. */
            if ( ( *argcount )[ *lpar - 1 ] !=
                 symbol[ ( *opensym )[ *lpar - 1 ] ].nargs ) {
               astError( AST__WRNFA,
                         "%s(%s): Wrong number of function arguments in the "
                         "expression \"%.*s\".", status,
                         method, class, iend + 1, exprs );

/* If the number of arguments is valid, decrement the parenthesis
   level. */
            } else {
               ( *lpar )--;
            }

/* If the symbol which opened this level of parenthesis was a function
   call with a variable number of arguments (indicated by a negative
   "nargs" entry in the symbol data), then we must check and process the
   number of function arguments. */
         } else if ( symbol[ ( *opensym )[ *lpar - 1 ] ].nargs < 0 ) {

/* Check that the minimum required number of arguments have been
   supplied. Report an error if they have not. */
            if ( ( *argcount )[ *lpar - 1 ] <
                 ( -symbol[ ( *opensym )[ *lpar - 1 ] ].nargs ) ) {
               astError( AST__WRNFA,
                         "%s(%s): Insufficient function arguments in the "
                         "expression \"%.*s\".", status,
                         method, class, iend + 1, exprs );

/* If the number of arguments is valid, increase the size of the
   constants array and check for errors. */
            } else {
               *con = astGrow( *con, *ncon + 1, sizeof( double ) );
               if ( astOK ) {

/* Append the argument count to the end of the array of constants and
   decrement the parenthesis level. */
                  ( *con )[ ( *ncon )++ ] =
                     (double) ( *argcount )[ --( *lpar ) ];
               }
            }

/* Finally, if the symbol which opened this level of parenthesis was
   not a function call ("nargs" entry in the symbol data is zero), then
   decrement the parenthesis level. In this case there is no need to
   check the argument count, because it will not have been
   incremented. */
         } else {
            ( *lpar )--;
         }
      }
   }

/* If an error occurred (or the global error status was set on entry),
   then reset the parenthesis level and free any memory which may have
   been allocated. */
   if ( !astOK ) {
      *lpar = 0;
      if ( *argcount ) *argcount = astFree( *argcount );
      if ( *opensym ) *opensym = astFree( *opensym );
      if ( *con ) *con = astFree( *con );
   }
}

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

*  Purpose:
*     Random number seed for a MathMap.

*  Type:
*     Public attribute.

*  Synopsis:
*     Integer.

*  Description:
*     This attribute, which may take any integer value, determines the
*     sequence of random numbers produced by the random number functions in
*     MathMap expressions. It is set to an unpredictable default value when
*     a MathMap is created, so that by default each MathMap uses a different
*     set of random numbers.
*
*     If required, you may set this Seed attribute to a value of your
*     choosing in order to produce repeatable behaviour from the random
*     number functions. You may also enquire the Seed value (e.g. if an
*     initially unpredictable value has been used) and then use it to
*     reproduce the resulting sequence of random numbers, either from the
*     same MathMap or from another one.
*
*     Clearing the Seed attribute gives it a new unpredictable default
*     value.

*  Applicability:
*     MathMap
*        All MathMaps have this attribute.
*att--
*/
/* Clear the Seed value by setting it to a new unpredictable value
   produced by DefaultSeed and clearing the "seed_set" flag in the
   MathMap's random number generator context. Also clear the "active"
   flag, so that the generator will be re-initialised to use this seed
   when it is next invoked. */
astMAKE_CLEAR(MathMap,Seed,rcontext.seed,( this->rcontext.seed_set = 0,
                                           this->rcontext.active = 0,
                                           DefaultSeed( &this->rcontext, status ) ))

/* Return the "seed" value from the random number generator
   context. */
astMAKE_GET(MathMap,Seed,int,0,this->rcontext.seed)

/* Store the new seed value in the MathMap's random number generator
   context and set the context's "seed_set" flag. Also clear the "active"
   flag, so that the generator will be re-initialised to use this seed
   when it is next invoked. */
astMAKE_SET(MathMap,Seed,int,rcontext.seed,( this->rcontext.seed_set = 1,
                                             this->rcontext.active = 0,
                                             value ))

/* Test the "seed_set" flag in the random number generator context. */
astMAKE_TEST(MathMap,Seed,( this->rcontext.seed_set ))

/*
*att++
*  Name:
*     SimpFI

*  Purpose:
*     Forward-inverse MathMap pairs simplify?

*  Type:
*     Public attribute.

*  Synopsis:
*     Integer (boolean).

*  Description:
c     This attribute should be set to a non-zero value if applying a
c     MathMap's forward transformation, followed immediately by the matching
c     inverse transformation will always restore the original set of
c     coordinates. It indicates that AST may replace such a sequence of
c     operations by an identity Mapping (a UnitMap) if it is encountered
c     while simplifying a compound Mapping (e.g. using astSimplify).
f     This attribute should be set to a non-zero value if applying a
f     MathMap's forward transformation, followed immediately by the matching
f     inverse transformation will always restore the original set of
f     coordinates. It indicates that AST may replace such a sequence of
f     operations by an identity Mapping (a UnitMap) if it is encountered
f     while simplifying a compound Mapping (e.g. using AST_SIMPLIFY).
*
*     By default, the SimpFI attribute is zero, so that AST will not perform
*     this simplification unless you have set SimpFI to indicate that it is
*     safe to do so.

*  Applicability:
*     MathMap
*        All MathMaps have this attribute.

*  Notes:
*     - For simplification to occur, the two MathMaps must be in series and
*     be identical (with textually identical transformation
*     functions). Functional equivalence is not sufficient.
*     - The consent of both MathMaps is required before simplification can
*     take place. If either has a SimpFI value of zero, then simplification
*     will not occur.
*     - The SimpFI attribute controls simplification only in the case where
*     a MathMap's forward transformation is followed by the matching inverse
*     transformation. It does not apply if an inverse transformation is
*     followed by a forward transformation. This latter case is controlled
*     by the SimpIF attribute.
c     - The "forward" and "inverse" transformations referred to are those
c     defined when the MathMap is created (corresponding to the "fwd" and
c     "inv" parameters of its constructor function). If the MathMap is
c     inverted (i.e. its Invert attribute is non-zero), then the role of the
c     SimpFI and SimpIF attributes will be interchanged.
f     - The "forward" and "inverse" transformations referred to are those
f     defined when the MathMap is created (corresponding to the FWD and
f     INV arguments of its constructor function). If the MathMap is
f     inverted (i.e. its Invert attribute is non-zero), then the role of the
f     SimpFI and SimpIF attributes will be interchanged.
*att--
*/
/* Clear the SimpFI value by setting it to -INT_MAX. */
astMAKE_CLEAR(MathMap,SimpFI,simp_fi,-INT_MAX)

/* Supply a default of 0 if no SimpFI value has been set. */
astMAKE_GET(MathMap,SimpFI,int,0,( ( this->simp_fi != -INT_MAX ) ?
                                   this->simp_fi : 0 ))

/* Set a SimpFI value of 1 if any non-zero value is supplied. */
astMAKE_SET(MathMap,SimpFI,int,simp_fi,( value != 0 ))

/* The SimpFI value is set if it is not -INT_MAX. */
astMAKE_TEST(MathMap,SimpFI,( this->simp_fi != -INT_MAX ))

/*
*att++
*  Name:
*     SimpIF

*  Purpose:
*     Inverse-forward MathMap pairs simplify?

*  Type:
*     Public attribute.

*  Synopsis:
*     Integer (boolean).

*  Description:
c     This attribute should be set to a non-zero value if applying a
c     MathMap's inverse transformation, followed immediately by the matching
c     forward transformation will always restore the original set of
c     coordinates. It indicates that AST may replace such a sequence of
c     operations by an identity Mapping (a UnitMap) if it is encountered
c     while simplifying a compound Mapping (e.g. using astSimplify).
f     This attribute should be set to a non-zero value if applying a
f     MathMap's inverse transformation, followed immediately by the matching
f     forward transformation will always restore the original set of
f     coordinates. It indicates that AST may replace such a sequence of
f     operations by an identity Mapping (a UnitMap) if it is encountered
f     while simplifying a compound Mapping (e.g. using AST_SIMPLIFY).
*
*     By default, the SimpIF attribute is zero, so that AST will not perform
*     this simplification unless you have set SimpIF to indicate that it is
*     safe to do so.

*  Applicability:
*     MathMap
*        All MathMaps have this attribute.

*  Notes:
*     - For simplification to occur, the two MathMaps must be in series and
*     be identical (with textually identical transformation
*     functions). Functional equivalence is not sufficient.
*     - The consent of both MathMaps is required before simplification can
*     take place. If either has a SimpIF value of zero, then simplification
*     will not occur.
*     - The SimpIF attribute controls simplification only in the case where
*     a MathMap's inverse transformation is followed by the matching forward
*     transformation. It does not apply if a forward transformation is
*     followed by an inverse transformation. This latter case is controlled
*     by the SimpFI attribute.
c     - The "forward" and "inverse" transformations referred to are those
c     defined when the MathMap is created (corresponding to the "fwd" and
c     "inv" parameters of its constructor function). If the MathMap is
c     inverted (i.e. its Invert attribute is non-zero), then the role of the
c     SimpFI and SimpIF attributes will be interchanged.
f     - The "forward" and "inverse" transformations referred to are those
f     defined when the MathMap is created (corresponding to the FWD and
f     INV arguments of its constructor function). If the MathMap is
f     inverted (i.e. its Invert attribute is non-zero), then the role of the
f     SimpFI and SimpIF attributes will be interchanged.
*att--
*/
/* Clear the SimpIF value by setting it to -INT_MAX. */
astMAKE_CLEAR(MathMap,SimpIF,simp_if,-INT_MAX)

/* Supply a default of 0 if no SimpIF value has been set. */
astMAKE_GET(MathMap,SimpIF,int,0,( ( this->simp_if != -INT_MAX ) ?
                                   this->simp_if : 0 ))

/* Set a SimpIF value of 1 if any non-zero value is supplied. */
astMAKE_SET(MathMap,SimpIF,int,simp_if,( value != 0 ))

/* The SimpIF value is set if it is not -INT_MAX. */
astMAKE_TEST(MathMap,SimpIF,( this->simp_if != -INT_MAX ))

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

*  Purpose:
*     Copy constructor for MathMap objects.

*  Type:
*     Private function.

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

*  Description:
*     This function implements the copy constructor for MathMap 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.
*/

/* Local Variables: */
   AstMathMap *in;               /* Pointer to input MathMap */
   AstMathMap *out;              /* Pointer to output MathMap */
   int ifun;                     /* Loop counter for functions */

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

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

/* For safety, first clear any references to the input memory from
   the output MathMap. */
   out->fwdfun = NULL;
   out->invfun = NULL;
   out->fwdcode = NULL;
   out->invcode = NULL;
   out->fwdcon = NULL;
   out->invcon = NULL;

/* Now allocate and initialise each of the output pointer arrays
   required. */
   if ( in->fwdfun ) {
      MALLOC_POINTER_ARRAY( out->fwdfun, char *, out->nfwd )
   }
   if ( in->invfun ) {
      MALLOC_POINTER_ARRAY( out->invfun, char *, out->ninv )
   }
   if ( in->fwdcode ) {
      MALLOC_POINTER_ARRAY( out->fwdcode, int *, out->nfwd )
   }
   if ( in->invcode ) {
      MALLOC_POINTER_ARRAY( out->invcode, int *, out->ninv )
   }
   if ( in->fwdcon ) {
      MALLOC_POINTER_ARRAY( out->fwdcon, double *, out->nfwd )
   }
   if ( in->invcon ) {
      MALLOC_POINTER_ARRAY( out->invcon, double *, out->ninv )
   }

/* If OK, loop to make copies of the data (where available) associated
   with each forward transformation function, storing pointers to the
   copy in the output pointer arrays allocated above. */
   if ( astOK ) {
      for ( ifun = 0; ifun < out->nfwd; ifun++ ) {
         if ( in->fwdfun && in->fwdfun[ ifun ] ) {
            out->fwdfun[ ifun ] = astStore( NULL, in->fwdfun[ ifun ],
                                            astSizeOf( in->fwdfun[ ifun ] ) );
         }
         if ( in->fwdcode && in->fwdcode[ ifun ] ) {
            out->fwdcode[ ifun ] = astStore( NULL, in->fwdcode[ ifun ],
                                            astSizeOf( in->fwdcode[ ifun ] ) );
         }
         if ( in->fwdcon && in->fwdcon[ ifun ] ) {
            out->fwdcon[ ifun ] = astStore( NULL, in->fwdcon[ ifun ],
                                            astSizeOf( in->fwdcon[ ifun ] ) );
         }
         if ( !astOK ) break;
      }
   }

/* Repeat this process for the inverse transformation functions. */
   if ( astOK ) {
      for ( ifun = 0; ifun < out->ninv; ifun++ ) {
         if ( in->invfun && in->invfun[ ifun ] ) {
            out->invfun[ ifun ] = astStore( NULL, in->invfun[ ifun ],
                                            astSizeOf( in->invfun[ ifun ] ) );
         }
         if ( in->invcode && in->invcode[ ifun ] ) {
            out->invcode[ ifun ] = astStore( NULL, in->invcode[ ifun ],
                                            astSizeOf( in->invcode[ ifun ] ) );
         }
         if ( in->invcon && in->invcon[ ifun ] ) {
            out->invcon[ ifun ] = astStore( NULL, in->invcon[ ifun ],
                                            astSizeOf( in->invcon[ ifun ] ) );
         }
         if ( !astOK ) break;
      }
   }

/* If an error occurred, clean up by freeing all output memory
   allocated above. */
   if ( !astOK ) {
      FREE_POINTER_ARRAY( out->fwdfun, out->nfwd )
      FREE_POINTER_ARRAY( out->invfun, out->ninv )
      FREE_POINTER_ARRAY( out->fwdcode, out->nfwd )
      FREE_POINTER_ARRAY( out->invcode, out->ninv )
      FREE_POINTER_ARRAY( out->fwdcon, out->nfwd )
      FREE_POINTER_ARRAY( out->invcon, out->ninv )
   }
}

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

*  Purpose:
*     Destructor for MathMap objects.

*  Type:
*     Private function.

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

*  Description:
*     This function implements the destructor for MathMap 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: */
   AstMathMap *this;             /* Pointer to MathMap */

/* Obtain a pointer to the MathMap structure. */
   this = (AstMathMap *) obj;

/* Free all memory allocated by the MathMap. */
   FREE_POINTER_ARRAY( this->fwdfun, this->nfwd )
   FREE_POINTER_ARRAY( this->invfun, this->ninv )
   FREE_POINTER_ARRAY( this->fwdcode, this->nfwd )
   FREE_POINTER_ARRAY( this->invcode, this->ninv )
   FREE_POINTER_ARRAY( this->fwdcon, this->nfwd )
   FREE_POINTER_ARRAY( this->invcon, this->ninv )
}

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

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

*  Parameters:
*     this
*        Pointer to the MathMap 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 COMMENT_LEN 150          /* Maximum length of a comment string */
#define KEY_LEN 50               /* Maximum length of a keyword */

/* Local Variables: */
   AstMathMap *this;             /* Pointer to the MathMap structure */
   char comment[ COMMENT_LEN + 1 ]; /* Buffer for comment strings */
   char key[ KEY_LEN + 1 ];      /* Buffer for keyword strings */
   int ifun;                     /* Loop counter for functions */
   int invert;                   /* MathMap inverted? */
   int ival;                     /* Integer attribute value */
   int nin;                      /* True number of input coordinates */
   int nout;                     /* True number of output coordinates */
   int set;                      /* Attribute value set? */

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

/* Obtain a pointer to the MathMap structure. */
   this = (AstMathMap *) this_object;

/* Determine if the MathMap is inverted and obtain the "true" number
   of input and output coordinates by un-doing the effects of any
   inversion. */
   invert = astGetInvert( this );
   nin = !invert ? astGetNin( this ) : astGetNout( this );
   nout = !invert ? astGetNout( this ) : astGetNin( this );

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

/* In the case of attributes, we first use the appropriate (private)
   Test...  member function to see if they are set. If so, we then use
   the (private) Get... function to obtain the value to be written
   out.

   For attributes which are not set, we use the astGet... method to
   obtain the value instead. This will supply a default value
   (possibly provided by a derived class which over-rides this method)
   which is more useful to a human reader as it corresponds to the
   actual default attribute value.  Since "set" will be zero, these
   values are for information only and will not be read back. */

/* Number of forward transformation functions. */
/* ------------------------------------------- */
/* We regard this value as set if it differs from the number of output
   coordinates for the MathMap. */
   set = ( this->nfwd != nout );
   astWriteInt( channel, "Nfwd", set, 0, this->nfwd,
                "Number of forward transformation functions" );

/* Forward transformation functions. */
/* --------------------------------- */
/* Loop to write out each forward transformation function, generating
   a suitable keyword and comment for each one. */
   for ( ifun = 0; ifun < this->nfwd; ifun++ ) {
      (void) sprintf( key, "Fwd%d", ifun + 1 );
      (void) sprintf( comment, "Forward function %d", ifun + 1 );
      astWriteString( channel, key, 1, 1, this->fwdfun[ ifun ], comment );
   }

/* Number of inverse transformation functions. */
/* ------------------------------------------- */
/* We regard this value as set if it differs from the number of input
   coordinates for the MathMap. */
   set = ( this->ninv != nin );
   astWriteInt( channel, "Ninv", set, 0, this->ninv,
                "Number of inverse transformation functions" );

/* Inverse transformation functions. */
/* --------------------------------- */
/* Similarly, loop to write out each inverse transformation
   function. */
   for ( ifun = 0; ifun < this->ninv; ifun++ ) {
      (void) sprintf( key, "Inv%d", ifun + 1 );
      (void) sprintf( comment, "Inverse function %d", ifun + 1 );
      astWriteString( channel, key, 1, 1, this->invfun[ ifun ], comment );
   }

/* SimpFI. */
/* ------- */
/* Write out the forward-inverse simplification flag. */
   set = TestSimpFI( this, status );
   ival = set ? GetSimpFI( this, status ) : astGetSimpFI( this );
   astWriteInt( channel, "SimpFI", set, 0, ival,
                ival ? "Forward-inverse pairs may simplify" :
                       "Forward-inverse pairs do not simplify" );

/* SimpIF. */
/* ------- */
/* Write out the inverse-forward simplification flag. */
   set = TestSimpIF( this, status );
   ival = set ? GetSimpIF( this, status ) : astGetSimpIF( this );
   astWriteInt( channel, "SimpIF", set, 0, ival,
                ival ? "Inverse-forward pairs may simplify" :
                       "Inverse-forward pairs do not simplify" );

/* Seed. */
/* ----- */
/* Write out any random number seed value which is set. Prefix this with
   a separate flag which indicates if the seed has been set. */
   set = TestSeed( this, status );
   ival = set ? GetSeed( this, status ) : astGetSeed( this );
   astWriteInt( channel, "Seeded", set, 0, set,
                set? "Explicit random number seed set" :
                     "No random number seed set" );
   astWriteInt( channel, "Seed", set, 0, ival,
                set ? "Random number seed value" :
                      "Default random number seed used" );

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

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

AstMathMap *astMathMap_( int nin, int nout,
                         int nfwd, const char *fwd[],
                         int ninv, const char *inv[],
                         const char *options, int *status, ...) {
/*
*+
*  Name:
*     astMathMap

*  Purpose:
*     Create a MathMap.

*  Type:
*     Protected function.

*  Synopsis:
*     #include "mathmap.h"
*     AstMathMap *astMathMap( int nin, int nout,
*                             int nfwd, const char *fwd[],
*                             int ninv, const char *inv[],
*                             const char *options, ..., int *status )

*  Class Membership:
*     MathMap constructor.

*  Description:
*     This function creates a new MathMap and optionally initialises its
*     attributes.

*  Parameters:
*     nin
*        Number of input variables for the MathMap.
*     nout
*        Number of output variables for the MathMap.
*     nfwd
*        The number of forward transformation functions being supplied.
*        This must be at least equal to "nout".
*     fwd
*        Pointer to an array, with "nfwd" elements, of pointers to null
*        terminated strings which contain each of the forward transformation
*        functions.
*     ninv
*        The number of inverse transformation functions being supplied.
*        This must be at least equal to "nin".
*     inv
*        Pointer to an array, with "ninv" elements, of pointers to null
*        terminated strings which contain each of the inverse transformation
*        functions.
*     options
*        Pointer to a null terminated string containing an optional
*        comma-separated list of attribute assignments to be used for
*        initialising the new MathMap. The syntax used is the same as
*        for the astSet method and may include "printf" format
*        specifiers identified by "%" symbols in the normal way.
*     status
*        Pointer to the inherited status variable.
*     ...
*        If the "options" string contains "%" format specifiers, then
*        an optional list of arguments may follow it in order to
*        supply values to be substituted for these specifiers. The
*        rules for supplying these are identical to those for the
*        astSet method (and for the C "printf" function).

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

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

*  Implementation Notes:
*     - This function implements the basic MathMap constructor which is
*     available via the protected interface to the MathMap class.  A
*     public interface is provided by the astMathMapId_ function.
*/

/* Local Variables: */
   astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
   AstMathMap *new;              /* Pointer to new MathMap */
   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 MathMap, allocating memory and initialising the
   virtual function table as well if necessary. */
   new = astInitMathMap( NULL, sizeof( AstMathMap ), !class_init, &class_vtab,
                         "MathMap", nin, nout, nfwd, fwd, ninv, inv );

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

AstMathMap *astMathMapId_( int nin, int nout,
                           int nfwd, const char *fwd[],
                           int ninv, const char *inv[],
                           const char *options, ... ) {
/*
*++
*  Name:
c     astMathMap
f     AST_MATHMAP

*  Purpose:
*     Create a MathMap.

*  Type:
*     Public function.

*  Synopsis:
c     #include "mathmap.h"
c     AstMathMap *astMathMap( int nin, int nout,
c                             int nfwd, const char *fwd[],
c                             int ninv, const char *inv[],
c                             const char *options, ... )
f     RESULT = AST_MATHMAP( NIN, NOUT, NFWD, FWD, NINV, INV, OPTIONS, STATUS )

*  Class Membership:
*     MathMap constructor.

*  Description:
*     This function creates a new MathMap and optionally initialises its
*     attributes.
*
c     A MathMap is a Mapping which allows you to specify a set of forward
c     and/or inverse transformation functions using arithmetic operations
c     and mathematical functions similar to those available in C. The
c     MathMap interprets these functions at run-time, whenever its forward
c     or inverse transformation is required. Because the functions are not
c     compiled in the normal sense (unlike an IntraMap), they may be used to
c     describe coordinate transformations in a transportable manner. A
c     MathMap therefore provides a flexible way of defining new types of
c     Mapping whose descriptions may be stored as part of a dataset and
c     interpreted by other programs.
f     A MathMap is a Mapping which allows you to specify a set of forward
f     and/or inverse transformation functions using arithmetic operations
f     and mathematical functions similar to those available in Fortran. The
f     MathMap interprets these functions at run-time, whenever its forward
f     or inverse transformation is required. Because the functions are not
f     compiled in the normal sense (unlike an IntraMap), they may be used to
f     describe coordinate transformations in a transportable manner. A
f     MathMap therefore provides a flexible way of defining new types of
f     Mapping whose descriptions may be stored as part of a dataset and
f     interpreted by other programs.

*  Parameters:
c     nin
f     NIN = INTEGER
*        Number of input variables for the MathMap. This determines the
*        value of its Nin attribute.
c     nout
f     NOUT = INTEGER
*        Number of output variables for the MathMap. This determines the
*        value of its Nout attribute.
c     nfwd
f     NFWD = INTEGER
*        The number of forward transformation functions being supplied.
c        This must be at least equal to "nout", but may be increased to
f        This must be at least equal to NOUT, but may be increased to
*        accommodate any additional expressions which define intermediate
*        variables for the forward transformation (see the "Calculating
*        Intermediate Values" section below).
c     fwd
f     FWD = CHARACTER * ( * )( NFWD )
c        An array (with "nfwd" elements) of pointers to null terminated strings
c        which contain the expressions defining the forward transformation.
f        An array which contains the expressions defining the forward
f        transformation.
*        The syntax of these expressions is described below.
c     ninv
f     NINV = INTEGER
*        The number of inverse transformation functions being supplied.
c        This must be at least equal to "nin", but may be increased to
f        This must be at least equal to NIN, but may be increased to
*        accommodate any additional expressions which define intermediate
*        variables for the inverse transformation (see the "Calculating
*        Intermediate Values" section below).
c     inv
f     INV = CHARACTER * ( * )( NINV )
c        An array (with "ninv" elements) of pointers to null terminated strings
c        which contain the expressions defining the inverse transformation.
f        An array which contains the expressions defining the inverse
f        transformation.
*        The syntax of these expressions is described below.
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 MathMap. 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.
c        If no initialisation is required, a zero-length string may be
c        supplied.
f        A character string containing an optional comma-separated
f        list of attribute assignments to be used for initialising the
f        new MathMap. The syntax used is identical to that for the
f        AST_SET routine. If no initialisation is required, a blank
f        value may be supplied.
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     astMathMap()
f     AST_MATHMAP = INTEGER
*        A pointer to the new MathMap.

*  Defining Transformation Functions:
c     A MathMap's transformation functions are supplied as a set of
c     expressions in an array of character strings. Normally you would
c     supply the same number of expressions for the forward transformation,
c     via the "fwd" parameter, as there are output variables (given by the
c     MathMap's Nout attribute). For instance, if Nout is 2 you might use:
c     - "r = sqrt( x * x + y * y )"
c     - "theta = atan2( y, x )"
c
c     which defines a transformation from Cartesian to polar
c     coordinates. Here, the variables that appear on the left of each
c     expression ("r" and "theta") provide names for the output variables
c     and those that appear on the right ("x" and "y") are references to
c     input variables.
f     A MathMap's transformation functions are supplied as a set of
f     expressions in an array of character strings. Normally you would
f     supply the same number of expressions for the forward transformation,
f     via the FWD argument, as there are output variables (given by the
f     MathMap's Nout attribute). For instance, if Nout is 2 you might use:
f     - 'R = SQRT( X * X + Y * Y )'
f     - 'THETA = ATAN2( Y, X )'
f
f     which defines a transformation from Cartesian to polar
f     coordinates. Here, the variables that appear on the left of each
f     expression (R and THETA) provide names for the output variables and
f     those that appear on the right (X and Y) are references to input
f     variables.
*
c     To complement this, you must also supply expressions for the inverse
c     transformation via the "inv" parameter.  In this case, the number of
c     expressions given would normally match the number of MathMap input
c     coordinates (given by the Nin attribute).  If Nin is 2, you might use:
c     - "x = r * cos( theta )"
c     - "y = r * sin( theta )"
c
c     which expresses the transformation from polar to Cartesian
c     coordinates. Note that here the input variables ("x" and "y") are
c     named on the left of each expression, and the output variables ("r"
c     and "theta") are referenced on the right.
f     To complement this, you must also supply expressions for the inverse
f     transformation via the INV argument.  In this case, the number of
f     expressions given would normally match the number of MathMap input
f     coordinates (given by the Nin attribute).  If Nin is 2, you might use:
f     - 'X = R * COS( THETA )'
f     - 'Y = R * SIN( THETA )'
f
f     which expresses the transformation from polar to Cartesian
f     coordinates. Note that here the input variables (X and Y) are named on
f     the left of each expression, and the output variables (R and THETA)
f     are referenced on the right.
*
*     Normally, you cannot refer to a variable on the right of an expression
*     unless it is named on the left of an expression in the complementary
*     set of functions. Therefore both sets of functions (forward and
*     inverse) must be formulated using the same consistent set of variable
*     names. This means that if you wish to leave one of the transformations
*     undefined, you must supply dummy expressions which simply name each of
*     the output (or input) variables.  For example, you might use:
c     - "x"
c     - "y"
f     - 'X'
f     - 'Y'
*
*     for the inverse transformation above, which serves to name the input
*     variables but without defining an inverse transformation.

*  Calculating Intermediate Values:
c     It is sometimes useful to calculate intermediate values and then to
c     use these in the final expressions for the output (or input)
c     variables. This may be done by supplying additional expressions for
c     the forward (or inverse) transformation functions. For instance, the
c     following array of five expressions describes 2-dimensional pin-cushion
c     distortion:
c     - "r = sqrt( xin * xin + yin * yin )"
c     - "rout = r * ( 1 + 0.1 * r * r )"
c     - "theta = atan2( yin, xin )"
c     - "xout = rout * cos( theta )"
c     - "yout = rout * sin( theta )"
f     It is sometimes useful to calculate intermediate values and then to
f     use these in the final expressions for the output (or input)
f     variables. This may be done by supplying additional expressions for
f     the forward (or inverse) transformation functions. For instance, the
f     following array of five expressions describes 2-dimensional pin-cushion
f     distortion:
f     - 'R = SQRT( XIN * XIN + YIN * YIN )'
f     - 'ROUT = R * ( 1 + 0.1 * R * R )'
f     - 'THETA = ATAN2( YIN, XIN )',
f     - 'XOUT = ROUT * COS( THETA )'
f     - 'YOUT = ROUT * SIN( THETA )'
*
c     Here, we first calculate three intermediate results ("r", "rout"
c     and "theta") and then use these to calculate the final results ("xout"
c     and "yout"). The MathMap knows that only the final two results
c     constitute values for the output variables because its Nout attribute
c     is set to 2. You may define as many intermediate variables in this
c     way as you choose. Having defined a variable, you may then refer to it
c     on the right of any subsequent expressions.
f     Here, we first calculate three intermediate results (R, ROUT
f     and THETA) and then use these to calculate the final results (XOUT
f     and YOUT). The MathMap knows that only the final two results
f     constitute values for the output variables because its Nout attribute
f     is set to 2. You may define as many intermediate variables in this
f     way as you choose. Having defined a variable, you may then refer to it
f     on the right of any subsequent expressions.
*
c     Note that when defining the inverse transformation you may only refer
c     to the output variables "xout" and "yout".  The intermediate variables
c     "r", "rout" and "theta" (above) are private to the forward
c     transformation and may not be referenced by the inverse
c     transformation. The inverse transformation may, however, define its
c     own private intermediate variables.
f     Note that when defining the inverse transformation you may only refer
f     to the output variables XOUT and YOUT. The intermediate variables R,
f     ROUT and THETA (above) are private to the forward transformation and
f     may not be referenced by the inverse transformation. The inverse
f     transformation may, however, define its own private intermediate
f     variables.

*  Expression Syntax:
c     The expressions given for the forward and inverse transformations
c     closely follow the syntax of the C programming language (with some
c     extensions for compatibility with Fortran). They may contain
c     references to variables and literal constants, together with
c     arithmetic, boolean, relational and bitwise operators, and function
c     invocations. A set of symbolic constants is also available. Each of
c     these is described in detail below. Parentheses may be used to
c     over-ride the normal order of evaluation. There is no built-in limit
c     to the length of expressions and they are insensitive to case or the
c     presence of additional white space.
f     The expressions given for the forward and inverse transformations
f     closely follow the syntax of Fortran (with some extensions for
f     compatibility with the C language). They may contain references to
f     variables and literal constants, together with arithmetic, logical,
f     relational and bitwise operators, and function invocations. A set of
f     symbolic constants is also available. Each of these is described in
f     detail below. Parentheses may be used to over-ride the normal order of
f     evaluation. There is no built-in limit to the length of expressions
f     and they are insensitive to case or the presence of additional white
f     space.

*  Variables:
*     Variable names must begin with an alphabetic character and may contain
*     only alphabetic characters, digits, and the underscore character
*     "_". There is no built-in limit to the length of variable names.

*  Literal Constants:
c     Literal constants, such as "0", "1", "0.007" or "2.505e-16" may appear
c     in expressions, with the decimal point and exponent being optional (a
c     "D" may also be used as an exponent character for compatibility with
c     Fortran). A unary minus "-" may be used as a prefix.
f     Literal constants, such as "0", "1", "0.007" or "2.505E-16" may appear
f     in expressions, with the decimal point and exponent being optional (a
f     "D" may also be used as an exponent character). A unary minus "-" may
f     be used as a prefix.

*  Arithmetic Precision:
*     All arithmetic is floating point, performed in double precision.

*  Propagation of Missing Data:
*     Unless indicated otherwise, if any argument of a function or operator
*     has the value AST__BAD (indicating missing data), then the result of
*     that function or operation is also AST__BAD, so that such values are
*     propagated automatically through all operations performed by MathMap
*     transformations.  The special value AST__BAD can be represented in
*     expressions by the symbolic constant "<bad>".
*
*     A <bad> result (i.e. equal to AST__BAD) is also produced in response
*     to any numerical error (such as division by zero or numerical
*     overflow), or if an invalid argument value is provided to a function
*     or operator.

*  Arithmetic Operators:
*     The following arithmetic operators are available:
c     - x1 + x2: Sum of "x1" and "x2".
f     - X1 + X2: Sum of X1 and X2.
c     - x1 - x2: Difference of "x1" and "x2".
f     - X1 - X2: Difference of X1 and X2.
c     - x1 * x2: Product of "x1" and "x1".
f     - X1 * X2: Product of X1 and X2.
c     - x1 / x2: Ratio of "x1" and "x2".
f     - X1 / X2: Ratio of X1 and X2.
c     - x1 ** x2: "x1" raised to the power of "x2".
f     - X1 ** X2: X1 raised to the power of X2.
c     - + x: Unary plus, has no effect on its argument.
f     - + X: Unary plus, has no effect on its argument.
c     - - x: Unary minus, negates its argument.
f     - - X: Unary minus, negates its argument.

c  Boolean Operators:
f  Logical Operators:
c     Boolean values are represented using zero to indicate false and
c     non-zero to indicate true. In addition, the value AST__BAD is taken to
c     mean "unknown". The values returned by boolean operators may therefore
c     be 0, 1 or AST__BAD. Where appropriate, "tri-state" logic is
c     implemented. For example, "a||b" may evaluate to 1 if "a" is non-zero,
c     even if "b" has the value AST__BAD. This is because the result of the
c     operation would not be affected by the value of "b", so long as "a" is
c     non-zero.
f     Logical values are represented using zero to indicate .FALSE. and
f     non-zero to indicate .TRUE.. In addition, the value AST__BAD is taken to
f     mean "unknown". The values returned by logical operators may therefore
f     be 0, 1 or AST__BAD. Where appropriate, "tri-state" logic is
f     implemented. For example, A.OR.B may evaluate to 1 if A is non-zero,
f     even if B has the value AST__BAD. This is because the result of the
f     operation would not be affected by the value of B, so long as A is
f     non-zero.
*
c     The following boolean operators are available:
f     The following logical operators are available:
c     - x1 && x2: Boolean AND between "x1" and "x2", returning 1 if both "x1"
c     and "x2" are non-zero, and 0 otherwise. This operator implements
c     tri-state logic. (The synonym ".and." is also provided for compatibility
c     with Fortran.)
f     - X1 .AND. X2: Logical AND between X1 and X2, returning 1 if both X1
f     and X2 are non-zero, and 0 otherwise. This operator implements
f     tri-state logic. (The synonym "&&" is also provided for compatibility
f     with C.)
c     - x1 || x2: Boolean OR between "x1" and "x2", returning 1 if either "x1"
c     or "x2" are non-zero, and 0 otherwise. This operator implements
c     tri-state logic. (The synonym ".or." is also provided for compatibility
c     with Fortran.)
f     - X1 .OR. X2: Logical OR between X1 and X2, returning 1 if either X1
f     or X2 are non-zero, and 0 otherwise. This operator implements
f     tri-state logic. (The synonym "||" is also provided for compatibility
f     with C.)
c     - x1 ^^ x2: Boolean exclusive OR (XOR) between "x1" and "x2", returning
c     1 if exactly one of "x1" and "x2" is non-zero, and 0 otherwise. Tri-state
c     logic is not used with this operator. (The synonyms ".neqv." and ".xor."
c     are also provided for compatibility with Fortran, although the second
c     of these is not standard.)
f     - X1 .NEQV. X2: Logical exclusive OR (XOR) between X1 and X2,
f     returning 1 if exactly one of X1 and X2 is non-zero, and 0
f     otherwise. Tri-state logic is not used with this operator. (The
f     synonym ".XOR." is also provided, although this is not standard
f     Fortran. In addition, the C-like synonym "^^" may be used, although
f     this is also not standard.)
c     - x1 .eqv. x2: This is provided only for compatibility with Fortran
c     and tests whether the boolean states of "x1" and "x2" (i.e. true/false)
c     are equal. It is the negative of the exclusive OR (XOR) function.
c     Tri-state logic is not used with this operator.
f     - X1 .EQV. X2: Tests whether the logical states of X1 and X2
f     (i.e. .TRUE./.FALSE.) are equal. It is the negative of the exclusive OR
f     (XOR) function.  Tri-state logic is not used with this operator.
c     - ! x: Boolean unary NOT operation, returning 1 if "x" is zero, and
c     0 otherwise. (The synonym ".not." is also provided for compatibility
c     with Fortran.)
f     - .NOT. X: Logical unary NOT operation, returning 1 if X is zero, and
f     0 otherwise. (The synonym "!" is also provided for compatibility with
f     C.)

*  Relational Operators:
c     Relational operators return the boolean result (0 or 1) of comparing
c     the values of two floating point values for equality or inequality. The
c     value AST__BAD may also be returned if either argument is <bad>.
f     Relational operators return the logical result (0 or 1) of comparing
f     the values of two floating point values for equality or inequality. The
f     value AST__BAD may also be returned if either argument is <bad>.
*
*     The following relational operators are available:
c     - x1 == x2: Tests whether "x1" equals "x1". (The synonym ".eq." is
c     also provided for compatibility with Fortran.)
f     - X1 .EQ. X2: Tests whether X1 equals X2. (The synonym "==" is also
f     provided for compatibility with C.)
c     - x1 != x2: Tests whether "x1" is unequal to "x2". (The synonym ".ne."
c     is also provided for compatibility with Fortran.)
f     - X1 .NE. X2: Tests whether X1 is unequal to X2. (The synonym "!=" is
f     also provided for compatibility with C.)
c     - x1 > x2: Tests whether "x1" is greater than "x2". (The synonym
c     ".gt." is also provided for compatibility with Fortran.)
f     - X1 .GT. X2: Tests whether X1 is greater than X2. (The synonym ">" is
f     also provided for compatibility with C.)
c     - x1 >= x2: Tests whether "x1" is greater than or equal to "x2". (The
c     synonym ".ge."  is also provided for compatibility with Fortran.)
f     - X1 .GE. X2: Tests whether X1 is greater than or equal to X2. (The
f     synonym ">=" is also provided for compatibility with C.)
c     - x1 < x2: Tests whether "x1" is less than "x2". (The synonym ".lt."
c     is also provided for compatibility with Fortran.)
f     - X1 .LT. X2: Tests whether X1 is less than X2. (The synonym "<" is also
f     provided for compatibility with C.)
c     - x1 <= x2: Tests whether "x1" is less than or equal to "x2". (The
c     synonym ".le." is also provided for compatibility with Fortran.)
f     - X1 .LE. X2: Tests whether X1 is less than or equal to X2. (The synonym
f     "<=" is also provided for compatibility with C.)
*
c     Note that relational operators cannot usefully be used to compare
c     values with the <bad> value (representing missing data), because the
c     result is always <bad>. The isbad() function should be used instead.
f     Note that relational operators cannot usefully be used to compare
f     values with the <bad> value (representing missing data), because the
f     result is always <bad>. The ISBAD() function should be used instead.
f
f     Note, also, that because logical operators can operate on floating
f     point values, care must be taken to use parentheses in some cases
f     where they would not normally be required in Fortran. For example,
f     the expresssion:
f     - .NOT. A .EQ. B
f
f     must be written:
f     - .NOT. ( A .EQ. B )
f
f     to prevent the .NOT. operator from associating with the variable A.

*  Bitwise Operators:
c     The bitwise operators provided by C are often useful when operating on
c     raw data (e.g. from instruments), so they are also provided for use in
c     MathMap expressions. In this case, however, the values on which they
c     operate are floating point values rather than pure integers. In order
c     to produce results which match the pure integer case, the operands are
c     regarded as fixed point binary numbers (i.e. with the binary
c     equivalent of a decimal point) with negative numbers represented using
c     twos-complement notation. For integer values, the resulting bit
c     pattern corresponds to that of the equivalent signed integer (digits
c     to the right of the point being zero). Operations on the bits
c     representing the fractional part are also possible, however.
f     Bitwise operators are often useful when operating on raw data
f     (e.g. from instruments), so they are provided for use in MathMap
f     expressions. In this case, however, the values on which they operate
f     are floating point values rather than the more usual pure integers. In
f     order to produce results which match the pure integer case, the
f     operands are regarded as fixed point binary numbers (i.e. with the
f     binary equivalent of a decimal point) with negative numbers
f     represented using twos-complement notation. For integer values, the
f     resulting bit pattern corresponds to that of the equivalent signed
f     integer (digits to the right of the point being zero). Operations on
f     the bits representing the fractional part are also possible, however.
*
*     The following bitwise operators are available:
c     - x1 >> x2: Rightward bit shift. The integer value of "x2" is taken
c     (rounding towards zero) and the bits representing "x1" are then
c     shifted this number of places to the right (or to the left if the
c     number of places is negative). This is equivalent to dividing "x1" by
c     the corresponding power of 2.
f     - X1 >> X2: Rightward bit shift. The integer value of X2 is taken
f     (rounding towards zero) and the bits representing X1 are then
f     shifted this number of places to the right (or to the left if the
f     number of places is negative). This is equivalent to dividing X1 by
f     the corresponding power of 2.
c     - x1 << x2: Leftward bit shift. The integer value of "x2" is taken
c     (rounding towards zero), and the bits representing "x1" are then
c     shifted this number of places to the left (or to the right if the
c     number of places is negative). This is equivalent to multiplying "x1"
c     by the corresponding power of 2.
f     - X1 << X2: Leftward bit shift. The integer value of X2 is taken
f     (rounding towards zero), and the bits representing X1 are then
f     shifted this number of places to the left (or to the right if the
f     number of places is negative). This is equivalent to multiplying X1
f     by the corresponding power of 2.
c     - x1 & x2: Bitwise AND between the bits of "x1" and those of "x2"
c     (equivalent to a boolean AND applied at each bit position in turn).
f     - X1 & X2: Bitwise AND between the bits of X1 and those of X2
f     (equivalent to a logical AND applied at each bit position in turn).
c     - x1 | x2: Bitwise OR between the bits of "x1" and those of "x2"
c     (equivalent to a boolean OR applied at each bit position in turn).
f     - X1 | X2: Bitwise OR between the bits of X1 and those of X2
f     (equivalent to a logical OR applied at each bit position in turn).
c     - x1 ^ x2: Bitwise exclusive OR (XOR) between the bits of "x1" and
c     those of "x2" (equivalent to a boolean XOR applied at each bit
c     position in turn).
f     - X1 ^ X2: Bitwise exclusive OR (XOR) between the bits of X1 and
f     those of X2 (equivalent to a logical XOR applied at each bit
f     position in turn).
*
c     Note that no bit inversion operator ("~" in C) is provided. This is
c     because inverting the bits of a twos-complement fixed point binary
c     number is equivalent to simply negating it. This differs from the
c     pure integer case because bits to the right of the binary point are
c     also inverted. To invert only those bits to the left of the binary
c     point, use a bitwise exclusive OR with the value -1 (i.e. "x^-1").
f     Note that no bit inversion operator is provided. This is
f     because inverting the bits of a twos-complement fixed point binary
f     number is equivalent to simply negating it. This differs from the
f     pure integer case because bits to the right of the binary point are
f     also inverted. To invert only those bits to the left of the binary
f     point, use a bitwise exclusive OR with the value -1 (i.e. X^-1).

*  Functions:
*     The following functions are available:
c     - abs(x): Absolute value of "x" (sign removal), same as fabs(x).
f     - ABS(X): Absolute value of X (sign removal), same as FABS(X).
c     - acos(x): Inverse cosine of "x", in radians.
f     - ACOS(X): Inverse cosine of X, in radians.
c     - acosd(x): Inverse cosine of "x", in degrees.
f     - ACOSD(X): Inverse cosine of X, in degrees.
c     - acosh(x): Inverse hyperbolic cosine of "x".
f     - ACOSH(X): Inverse hyperbolic cosine of X.
c     - acoth(x): Inverse hyperbolic cotangent of "x".
f     - ACOTH(X): Inverse hyperbolic cotangent of X.
c     - acsch(x): Inverse hyperbolic cosecant of "x".
f     - ACSCH(X): Inverse hyperbolic cosecant of X.
c     - aint(x): Integer part of "x" (round towards zero), same as int(x).
f     - AINT(X): Integer part of X (round towards zero), same as INT(X).
c     - asech(x): Inverse hyperbolic secant of "x".
f     - ASECH(X): Inverse hyperbolic secant of X.
c     - asin(x): Inverse sine of "x", in radians.
f     - ASIN(X): Inverse sine of X, in radians.
c     - asind(x): Inverse sine of "x", in degrees.
f     - ASIND(X): Inverse sine of X, in degrees.
c     - asinh(x): Inverse hyperbolic sine of "x".
f     - ASINH(X): Inverse hyperbolic sine of X.
c     - atan(x): Inverse tangent of "x", in radians.
f     - ATAN(X): Inverse tangent of X, in radians.
c     - atand(x): Inverse tangent of "x", in degrees.
f     - ATAND(X): Inverse tangent of X, in degrees.
c     - atanh(x): Inverse hyperbolic tangent of "x".
f     - ATANH(X): Inverse hyperbolic tangent of X.
c     - atan2(x1, x2): Inverse tangent of "x1/x2", in radians.
f     - ATAN2(X1, X2): Inverse tangent of X1/X2, in radians.
c     - atan2d(x1, x2): Inverse tangent of "x1/x2", in degrees.
f     - ATAN2D(X1, X2): Inverse tangent of X1/X2, in degrees.
c     - ceil(x): Smallest integer value not less then "x" (round towards
c       plus infinity).
f     - CEIL(X): Smallest integer value not less then X (round towards
f       plus infinity).
c     - cos(x): Cosine of "x" in radians.
f     - COS(X): Cosine of X in radians.
c     - cosd(x): Cosine of "x" in degrees.
f     - COSD(X): Cosine of X in degrees.
c     - cosh(x): Hyperbolic cosine of "x".
f     - COSH(X): Hyperbolic cosine of X.
c     - coth(x): Hyperbolic cotangent of "x".
f     - COTH(X): Hyperbolic cotangent of X.
c     - csch(x): Hyperbolic cosecant of "x".
f     - CSCH(X): Hyperbolic cosecant of X.
c     - dim(x1, x2): Returns "x1-x2" if "x1" is greater than "x2", otherwise 0.
f     - DIM(X1, X2): Returns X1-X2 if X1 is greater than X2, otherwise 0.
c     - exp(x): Exponential function of "x".
f     - EXP(X): Exponential function of X.
c     - fabs(x): Absolute value of "x" (sign removal), same as abs(x).
f     - FABS(X): Absolute value of X (sign removal), same as ABS(X).
c     - floor(x): Largest integer not greater than "x" (round towards
c       minus infinity).
f     - FLOOR(X): Largest integer not greater than X (round towards
f       minus infinity).
c     - fmod(x1, x2): Remainder when "x1" is divided by "x2", same as
c       mod(x1, x2).
f     - FMOD(X1, X2): Remainder when X1 is divided by X2, same as
f       MOD(X1, X2).
c     - gauss(x1, x2): Random sample from a Gaussian distribution with mean
c       "x1" and standard deviation "x2".
f     - GAUSS(X1, X2): Random sample from a Gaussian distribution with mean
f       X1 and standard deviation X2.
c     - int(x): Integer part of "x" (round towards zero), same as aint(x).
f     - INT(X): Integer part of X (round towards zero), same as AINT(X).
c     - isbad(x): Returns 1 if "x" has the <bad> value (AST__BAD), otherwise 0.
f     - ISBAD(X): Returns 1 if X has the <bad> value (AST__BAD), otherwise 0.
c     - log(x): Natural logarithm of "x".
f     - LOG(X): Natural logarithm of X.
c     - log10(x): Logarithm of "x" to base 10.
f     - LOG10(X): Logarithm of X to base 10.
c     - max(x1, x2, ...): Maximum of two or more values.
f     - MAX(X1, X2, ...): Maximum of two or more values.
c     - min(x1, x2, ...): Minimum of two or more values.
f     - MIN(X1, X2, ...): Minimum of two or more values.
c     - mod(x1, x2): Remainder when "x1" is divided by "x2", same as
c       fmod(x1, x2).
f     - MOD(X1, X2): Remainder when X1 is divided by X2, same as
f       FMOD(X1, X2).
c     - nint(x): Nearest integer to "x" (round to nearest).
f     - NINT(X): Nearest integer to X (round to nearest).
c     - poisson(x): Random integer-valued sample from a Poisson
c       distribution with mean "x".
f     - POISSON(X): Random integer-valued sample from a Poisson
f       distribution with mean X.
c     - pow(x1, x2): "x1" raised to the power of "x2".
f     - POW(X1, X2): X1 raised to the power of X2.
c     - qif(x1, x2, x3): Returns "x2" if "x1" is true, and "x3" otherwise.
f     - QIF(x1, x2, x3): Returns X2 if X1 is true, and X3 otherwise.
c     - rand(x1, x2): Random sample from a uniform distribution in the
c       range "x1" to "x2" inclusive.
f     - RAND(X1, X2): Random sample from a uniform distribution in the
f       range X1 to X2 inclusive.
c     - sech(x): Hyperbolic secant of "x".
f     - SECH(X): Hyperbolic secant of X.
c     - sign(x1, x2): Absolute value of "x1" with the sign of "x2"
c       (transfer of sign).
f     - SIGN(X1, X2): Absolute value of X1 with the sign of X2
f       (transfer of sign).
c     - sin(x): Sine of "x" in radians.
f     - SIN(X): Sine of X in radians.
c     - sinc(x): Sinc function of "x" [= "sin(x)/x"].
f     - SINC(X): Sinc function of X [= SIN(X)/X].
c     - sind(x): Sine of "x" in degrees.
f     - SIND(X): Sine of X in degrees.
c     - sinh(x): Hyperbolic sine of "x".
f     - SINH(X): Hyperbolic sine of X.
c     - sqr(x): Square of "x" (= "x*x").
f     - SQR(X): Square of X (= X*X).
c     - sqrt(x): Square root of "x".
f     - SQRT(X): Square root of X.
c     - tan(x): Tangent of "x" in radians.
f     - TAN(X): Tangent of X in radians.
c     - tand(x): Tangent of "x" in degrees.
f     - TAND(X): Tangent of X in degrees.
c     - tanh(x): Hyperbolic tangent of "x".
f     - TANH(X): Hyperbolic tangent of X.

*  Symbolic Constants:
*     The following symbolic constants are available (the enclosing "<>"
*     brackets must be included):
c     - <bad>: The "bad" value (AST__BAD) used to flag missing data. Note
c     that you cannot usefully compare values with this constant because the
c     result is always <bad>. The isbad() function should be used instead.
f     - <bad>: The "bad" value (AST__BAD) used to flag missing data. Note
f     that you cannot usefully compare values with this constant because the
f     result is always <bad>. The ISBAD() function should be used instead.
c     - <dig>: Number of decimal digits of precision available in a
c     floating point (double) value.
f     - <dig>: Number of decimal digits of precision available in a
f     floating point (double precision) value.
*     - <e>: Base of natural logarithms.
*     - <epsilon>: Smallest positive number such that 1.0+<epsilon> is
*     distinguishable from unity.
c     - <mant_dig>: The number of base <radix> digits stored in the
c     mantissa of a floating point (double) value.
f     - <mant_dig>: The number of base <radix> digits stored in the
f     mantissa of a floating point (double precision) value.
c     - <max>: Maximum representable floating point (double) value.
f     - <max>: Maximum representable floating point (double precision) value.
c     - <max_10_exp>: Maximum integer such that 10 raised to that power
c     can be represented as a floating point (double) value.
f     - <max_10_exp>: Maximum integer such that 10 raised to that power
f     can be represented as a floating point (double precision) value.
c     - <max_exp>: Maximum integer such that <radix> raised to that
c     power minus 1 can be represented as a floating point (double) value.
f     - <max_exp>: Maximum integer such that <radix> raised to that
f     power minus 1 can be represented as a floating point (double precision)
f     value.
c     - <min>: Smallest positive number which can be represented as a
c     normalised floating point (double) value.
f     - <min>: Smallest positive number which can be represented as a
f     normalised floating point (double precision) value.
c     - <min_10_exp>: Minimum negative integer such that 10 raised to that
c     power can be represented as a normalised floating point (double) value.
f     - <min_10_exp>: Minimum negative integer such that 10 raised to that
f     power can be represented as a normalised floating point (double
f     precision) value.
c     - <min_exp>: Minimum negative integer such that <radix> raised to
c     that power minus 1 can be represented as a normalised floating point
c     (double) value.
f     - <min_exp>: Minimum negative integer such that <radix> raised to
f     that power minus 1 can be represented as a normalised floating point
f     (double precision) value.
*     - <pi>: Ratio of the circumference of a circle to its diameter.
c     - <radix>: The radix (number base) used to represent the mantissa of
c     floating point (double) values.
f     - <radix>: The radix (number base) used to represent the mantissa of
f     floating point (double precision) values.
*     - <rounds>: The mode used for rounding floating point results after
*     addition. Possible values include: -1 (indeterminate), 0 (toward
*     zero), 1 (to nearest), 2 (toward plus infinity) and 3 (toward minus
*     infinity). Other values indicate machine-dependent behaviour.

*  Evaluation Precedence and Associativity:
*     Items appearing in expressions are evaluated in the following order
*     (highest precedence first):
*     - Constants and variables
*     - Function arguments and parenthesised expressions
*     - Function invocations
*     - Unary + - ! .not.
*     - **
*     - * /
*     - + -
*     - << >>
*     - < .lt. <= .le. > .gt. >= .ge.
*     - == .eq. != .ne.
*     - &
*     - ^
*     - |
*     - && .and.
*     - ^^
*     - || .or
*     - .eqv. .neqv. .xor.
*
*     All operators associate from left-to-right, except for unary +,
*     unary -, !, .not. and ** which associate from right-to-left.

*  Notes:
*     - The sequence of numbers produced by the random number functions
*     available within a MathMap is normally unpredictable and different for
*     each MathMap. However, this behaviour may be controlled by means of
*     the MathMap's Seed attribute.
c     - Normally, compound Mappings (CmpMaps) which involve MathMaps will
c     not be subject to simplification (e.g. using astSimplify) because AST
c     cannot know how different MathMaps will interact. However, in the
c     special case where a MathMap occurs in series with its own inverse,
c     then simplification may be possible. Whether simplification does, in
c     fact, occur under these circumstances is controlled by the MathMap's
c     SimpFI and SimpIF attributes.
f     - Normally, compound Mappings (CmpMaps) which involve MathMaps will
f     not be subject to simplification (e.g. using AST_SIMPLIFY) because AST
f     cannot know how different MathMaps will interact. However, in the
f     special case where a MathMap occurs in series with its own inverse,
f     then simplification may be possible. Whether simplification does, in
f     fact, occur under these circumstances is controlled by the MathMap's
f     SimpFI and SimpIF attributes.
*     - 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.
*--

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

/* Local Variables: */
   astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
   AstMathMap *new;              /* Pointer to new MathMap */
   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 error status. */
   if ( !astOK ) return NULL;

/* Initialise the MathMap, allocating memory and initialising the virtual
   function table as well if necessary. */
   new = astInitMathMap( NULL, sizeof( AstMathMap ), !class_init, &class_vtab,
                         "MathMap", nin, nout, nfwd, fwd, ninv, inv );

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

AstMathMap *astInitMathMap_( void *mem, size_t size, int init,
                             AstMathMapVtab *vtab, const char *name,
                             int nin, int nout,
                             int nfwd, const char *fwd[],
                             int ninv, const char *inv[], int *status ) {
/*
*+
*  Name:
*     astInitMathMap

*  Purpose:
*     Initialise a MathMap.

*  Type:
*     Protected function.

*  Synopsis:
*     #include "mathmap.h"
*     AstMathMap *astInitMathMap_( void *mem, size_t size, int init,
*                                  AstMathMapVtab *vtab, const char *name,
*                                  int nin, int nout,
*                                  int nfwd, const char *fwd[],
*                                  int ninv, const char *inv[] )

*  Class Membership:
*     MathMap initialiser.

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

*  Parameters:
*     mem
*        A pointer to the memory in which the MathMap is to be initialised.
*        This must be of sufficient size to accommodate the MathMap data
*        (sizeof(MathMap)) 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 MathMap (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 MathMap
*        structure, so a valid value must be supplied even if not required for
*        allocating memory.
*     init
*        A logical flag indicating if the MathMap'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 MathMap.
*     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 Object
*        astClass function).
*     nin
*        Number of input variables for the MathMap.
*     nout
*        Number of output variables for the MathMap.
*     nfwd
*        The number of forward transformation functions being supplied.
*        This must be at least equal to "nout".
*     fwd
*        Pointer to an array, with "nfwd" elements, of pointers to null
*        terminated strings which contain each of the forward transformation
*        functions.
*     ninv
*        The number of inverse transformation functions being supplied.
*        This must be at least equal to "nin".
*     inv
*        Pointer to an array, with "ninv" elements, of pointers to null
*        terminated strings which contain each of the inverse transformation
*        functions.

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

*  Notes:
*     -  This function does not attempt to ensure that the forward and inverse
*     transformations performed by the resulting MathMap are consistent in any
*     way.
*     -  This function makes a copy of the contents of the strings supplied.
*     -  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: */
   AstMathMap *new;              /* Pointer to new MathMap */
   char **fwdfun;                /* Array of cleaned forward functions */
   char **invfun;                /* Array of cleaned inverse functions */
   double **fwdcon;              /* Constants for forward functions */
   double **invcon;              /* Constants for inverse functions */
   int **fwdcode;                /* Code for forward functions */
   int **invcode;                /* Code for inverse functions */
   int fwdstack;                 /* Stack size for forward functions */
   int invstack;                 /* Stack size for inverse functions */

/* Initialise. */
   new = NULL;

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

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

/* Check the numbers of input and output variables for validity,
   reporting an error if necessary. */
   if ( nin < 1 ) {
      astError( AST__BADNI,
                "astInitMathMap(%s): Bad number of input coordinates (%d).", status,
                name, nin );
      astError( AST__BADNI,
                "This number should be one or more." , status);
   } else if ( nout < 1 ) {
      astError( AST__BADNO,
                "astInitMathMap(%s): Bad number of output coordinates (%d).", status,
                name, nout );
      astError( AST__BADNI,
                "This number should be one or more." , status);

/* Check that sufficient number of forward and inverse transformation
   functions have been supplied and report an error if necessary. */
   } else if ( nfwd < nout ) {
      astError( AST__INNTF,
                "astInitMathMap(%s): Too few forward transformation functions "
                "given (%d).", status,
                name, nfwd );
      astError( astStatus,
                "At least %d forward transformation functions must be "
                "supplied. ", status,
                nout );
   } else if ( ninv < nin ) {
      astError( AST__INNTF,
                "astInitMathMap(%s): Too few inverse transformation functions "
                "given (%d).", status,
                name, ninv );
      astError( astStatus,
                "At least %d inverse transformation functions must be "
                "supplied. ", status,
                nin );

/* Of OK, clean the forward and inverse functions provided. This makes
   a lower-case copy with white space removed. */
   } else {
      CleanFunctions( nfwd, fwd, &fwdfun, status );
      CleanFunctions( ninv, inv, &invfun, status );

/* Compile the cleaned functions. From the returned pointers (if
   successful), we can now tell which transformations (forward and/or
   inverse) are defined. */
      CompileMapping( "astInitMathMap", name, nin, nout,
                      nfwd, (const char **) fwdfun,
                      ninv, (const char **) invfun,
                      &fwdcode, &invcode, &fwdcon, &invcon,
                      &fwdstack, &invstack, status );

/* Initialise a Mapping structure (the parent class) as the first
   component within the MathMap structure, allocating memory if
   necessary. Specify that the Mapping should be defined in the required
   directions. */
      new = (AstMathMap *) astInitMapping( mem, size, 0,
                                           (AstMappingVtab *) vtab, name,
                                           nin, nout,
                                           ( fwdcode != NULL ),
                                           ( invcode != NULL ) );


/* If an error has occurred, free all the memory which may have been
   allocated by the cleaning and compilation steps above. */
      if ( !astOK ) {
         FREE_POINTER_ARRAY( fwdfun, nfwd )
         FREE_POINTER_ARRAY( invfun, ninv )
         FREE_POINTER_ARRAY( fwdcode, nfwd )
         FREE_POINTER_ARRAY( invcode, ninv )
         FREE_POINTER_ARRAY( fwdcon, nfwd )
         FREE_POINTER_ARRAY( invcon, ninv )
      }

/* Initialise the MathMap data. */
/* ---------------------------- */
/* Store pointers to the compiled function information, together with
   other MathMap data. */
      if ( new ) {
         new->fwdfun = fwdfun;
         new->invfun = invfun;
         new->fwdcode = fwdcode;
         new->invcode = invcode;
         new->fwdcon = fwdcon;
         new->invcon = invcon;
         new->fwdstack = fwdstack;
         new->invstack = invstack;
         new->nfwd = nfwd;
         new->ninv = ninv;
         new->simp_fi = -INT_MAX;
         new->simp_if = -INT_MAX;

/* Initialise the random number generator context associated with the
   MathMap, using an unpredictable default seed value. */
         new->rcontext.active = 0;
         new->rcontext.random_int = 0;
         new->rcontext.seed_set = 0;
         new->rcontext.seed = DefaultSeed( &new->rcontext, status );

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

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

AstMathMap *astLoadMathMap_( void *mem, size_t size,
                             AstMathMapVtab *vtab, const char *name,
                             AstChannel *channel, int *status ) {
/*
*+
*  Name:
*     astLoadMathMap

*  Purpose:
*     Load a MathMap.

*  Type:
*     Protected function.

*  Synopsis:
*     #include "mathmap.h"
*     AstMathMap *astLoadMathMap( void *mem, size_t size,
*                                 AstMathMapVtab *vtab, const char *name,
*                                 AstChannel *channel )

*  Class Membership:
*     MathMap loader.

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


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

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

*  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: */
   AstMathMap *new;              /* Pointer to the new MathMap */
   char key[ KEY_LEN + 1 ];      /* Buffer for keyword strings */
   int ifun;                     /* Loop counter for functions */
   int invert;                   /* Invert attribute value */
   int nin;                      /* True number of input coordinates */
   int nout;                     /* True number of output coordinates */

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

/* If required, initialise the virtual function table for this class. */
      if ( !class_init ) {
         astInitMathMapVtab( 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 MathMap. */
   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, "MathMap" );

/* Determine if the MathMap is inverted and obtain the "true" number
   of input and output coordinates by un-doing the effects of any
   inversion. */
      invert = astGetInvert( new );
      nin = invert ? astGetNout( new ) : astGetNin( new );
      nout = invert ? astGetNin( new ) : astGetNout( new );

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

/* In the case of attributes, we first read the "raw" input value,
   supplying the "unset" value as the default. If a "set" value is
   obtained, we then use the appropriate (private) Set... member
   function to validate and set the value properly. */

/* Numbers of transformation functions. */
/* ------------------------------------ */
/* Read the numbers of forward and inverse transformation functions,
   supplying appropriate defaults. */
      new->nfwd = astReadInt( channel, "nfwd", nout );
      new->ninv = astReadInt( channel, "ninv", nin );
      if ( astOK ) {

/* Allocate memory for the MathMap's transformation function arrays. */
         MALLOC_POINTER_ARRAY( new->fwdfun, char *, new->nfwd )
         MALLOC_POINTER_ARRAY( new->invfun, char *, new->ninv )
         if ( astOK ) {

/* Forward transformation functions. */
/* --------------------------------- */
/* Create a keyword for each forward transformation function and read
   the function's value as a string. */
            for ( ifun = 0; ifun < new->nfwd; ifun++ ) {
               (void) sprintf( key, "fwd%d", ifun + 1 );
               new->fwdfun[ ifun ] = astReadString( channel, key, "" );
            }

/* Inverse transformation functions. */
/* --------------------------------- */
/* Repeat this process for the inverse transformation functions. */
            for ( ifun = 0; ifun < new->ninv; ifun++ ) {
               (void) sprintf( key, "inv%d", ifun + 1 );
               new->invfun[ ifun ] = astReadString( channel, key, "" );
            }

/* Forward-inverse simplification flag. */
/* ------------------------------------ */
            new->simp_fi = astReadInt( channel, "simpfi", -INT_MAX );
            if ( TestSimpFI( new, status ) ) SetSimpFI( new, new->simp_fi, status );

/* Inverse-forward simplification flag. */
/* ------------------------------------ */
            new->simp_if = astReadInt( channel, "simpif", -INT_MAX );
            if ( TestSimpIF( new, status ) ) SetSimpIF( new, new->simp_if, status );

/* Random number context. */
/* ---------------------- */
/* Initialise the random number generator context. */
            new->rcontext.active = 0;
            new->rcontext.random_int = 0;

/* Read the flag that determines if the Seed value is set, and the
   Seed value itself. */
            new->rcontext.seed_set = astReadInt( channel, "seeded", 0 );
            if ( TestSeed( new, status ) ) {
               new->rcontext.seed = astReadInt( channel, "seed", 0 );
               SetSeed( new, new->rcontext.seed, status );

/* Supply an unpredictable default Seed value if necessary. */
            } else {
               new->rcontext.seed = DefaultSeed( &new->rcontext, status );
            }

/* Compile the MathMap's transformation functions. */
            CompileMapping( "astLoadMathMap", name, nin, nout,
                            new->nfwd, (const char **) new->fwdfun,
                            new->ninv, (const char **) new->invfun,
                            &new->fwdcode, &new->invcode,
                            &new->fwdcon, &new->invcon,
                            &new->fwdstack, &new->invstack, status );
         }

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

/* Return the new MathMap pointer. */
   return new;

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