update ast

1 files changed, 0 insertions, 7421 deletions

diff --git a/ast/mathmap.c b/ast/mathmap.c
deleted file mode 100644
index cfcd033..0000000
--- a/ast/mathmap.c
+++ /dev/null
@@ -1,7421 +0,0 @@
-/*
-*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
-}
-
-
-
-
-