upgrade ast 8.7.1

1 files changed, 6218 insertions, 0 deletions

diff --git a/ast/unit.c b/ast/unit.c
new file mode 100644
index 0000000..0ce3184
--- /dev/null
+++ b/ast/unit.c
@@ -0,0 +1,6218 @@
+/*
+*  Name:
+*     unit.c
+
+*  Purpose:
+*     Implement unit conversion functions.
+
+*  Description:
+*     This file implements the Unit module which is used for identifying
+*     units and transforming between them. It follows the recommendations
+*     for unit handling contained within FITS WCS paper I (Greisen &
+*     Calabretta). All methods have protected access.
+
+*  Methods:
+*     astUnitMapper: Create a Mapping between two systems of units.
+*     astUnitLabel: Returns a label for a given unit symbol.
+
+*  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:
+*     DSB: D.S. Berry (Starlink)
+
+*  History:
+*     10-DEC-2002 (DSB):
+*        Original version.
+*     10-FEB-2004 (DSB):
+*        Added debug conditional code to keep track of memory leaks.
+*     15-JUL-2004 (DSB):
+*        In astUnitMapper: if no Mapping can be found from input to
+*        output units (e.g. because fo the less than perfect simplication
+*        algorithm in SimplifyTree), try finding a Mapping from output to
+*        input units and inverting the result.
+*     14-DEC-2004 (DSB):
+*        In CreateTree, move the invocation of FixConstants from after
+*        InvertConstants to before InvertConstants. This is because
+*        InvertConstants ignores nodes which contain all constant
+*        arguments. This results in constants not being inverted in
+*        expressions such as "1/60 deg" (because all arguments are
+*        constant in the the "1/60" node).
+*     18-JAN-2005 (DSB):
+*        Fix memory leaks.
+*     2-FEB-2005 (DSB):
+*        - Avoid using astStore to allocate more storage than is supplied
+*        in the "data" pointer. This can cause access violations since
+*        astStore will then read beyond the end of the "data" area.
+*     15-FEB-2005 (DSB):
+*        - Modified CleanExp to fix up some common units mistakes.
+*     21-FEB-2005 (DSB):
+*        - Modified CleanExp to accept <word><digit> as equivalent to
+*        <word>^<digit>.
+*        - Modified MakeTree to do case insensitive checking if case
+*        sensitive checking failsto produce a match to a multiplier/unit
+*        symbol combination. If this still produces no match, do a case
+*        insensitive match for multiplier/unit label.
+*     1-MAR-2006 (DSB):
+*        Replace astSetPermMap within DEBUG blocks by astBeginPM/astEndPM.
+*     6-APR-2006 (DSB):
+*        Modify CleanExp to convert "MJY/STER" to standard form ("MJy/sr").
+*     7-JUL-2006 (DSB):
+*        Correct initialisation of "word" flag in CleanExp.
+*     17-MAY-2007 (DSB):
+*        Simplify the units string returned by astUnitNormaliser.
+*        - Fix indexing bug in CombineFactors.
+*     26-MAY-2007 (DSB):
+*        - Correct error reporting in astUNitNormaliser.
+*        - Fix bug in CleanExp that caused (for example) "-2" to be
+*        converted to "^-2".
+*     2-DEC-2008 (DSB):
+*        Correct memory allocation bug in CleanExp.
+*     6-MAY-2011 (DSB):
+*        Include "adu" as basic unit.
+*     9-MAY-2011 (DSB):
+*        Change "A" to be Ampere (as defined by FITS-WCS paper 1) rather
+*        than "Angstrom".
+*/
+
+/* Module Macros. */
+/* ============== */
+/* Define the astCLASS macro (even although this is not a class
+   implementation) to obtain access to the protected error and memory
+   handling functions. */
+#define astCLASS
+#define PI 3.141592653589793
+#define E 2.718281828459045
+
+/* Macro which returns the nearest integer to a given floating point
+   value. */
+#define NINT(x) (int)((x)+(((x)>0.0)?0.5:-0.5))
+
+/* Macro identifying a character as lower or upper case letter, digit or
++ or -. */
+#define ISWORD(c) (isalnum(c)||((c)=='+')||((c)=='-'))
+
+/* The number of basic dimension quantities used for dimensional analysis.
+   In addition to the usual M, L and T, this includes pseudo-dimensions
+   describing strictly dimensionless quantities such as plane angle,
+   magnitude, etc. */
+#define NQUANT 10
+
+/* Include files. */
+/* ============== */
+/* Interface definitions. */
+/* ---------------------- */
+#include "error.h"
+#include "memory.h"
+#include "pointset.h"
+#include "mapping.h"
+#include "unitmap.h"
+#include "zoommap.h"
+#include "mathmap.h"
+#include "unit.h"
+
+/* Error code definitions. */
+/* ----------------------- */
+#include "ast_err.h"             /* AST error codes */
+
+/* C header files. */
+/* --------------- */
+#include <string.h>
+#include <stdio.h>
+#include <ctype.h>
+#include <limits.h>
+#include <math.h>
+
+#ifdef THREAD_SAFE
+#include <pthread.h>
+#endif
+
+/* Module Type Definitions. */
+/* ======================== */
+
+/* This declaration enumerates the codes for the operations which may
+   legally be included in a units expression. */
+typedef enum {
+   OP_LDCON,                     /* Load constant */
+   OP_LDVAR,                     /* Load variable */
+   OP_LOG,                       /* Common logarithm */
+   OP_LN,                        /* Natural logarithm */
+   OP_EXP,                       /* Natural exponential */
+   OP_SQRT,                      /* Square root */
+   OP_POW,                       /* Exponentiate */
+   OP_DIV,                       /* Division */
+   OP_MULT,                      /* Multiplication */
+   OP_LDPI,                      /* Load constant PI */
+   OP_LDE,                       /* Load constant E */
+   OP_NULL                       /* Null operation */
+} Oper;
+
+/* A structure describing a standard multiplier prefix. */
+typedef struct Multiplier {
+   const char *label;       /* Multipler label string (null terminated) */
+   const char *sym;         /* Multipler symbol string (null terminated) */
+   int symlen;              /* Length of symbol (without trailing null ) */
+   int lablen;              /* Length of label (without trailing null ) */
+   double scale;            /* The scale factor associated with the prefix */
+   struct Multiplier *next; /* Next Multiplier in linked list */
+} Multiplier;
+
+/* A structure describing a single node in a tree representation of a
+   single units expression. */
+typedef struct UnitNode {
+   Oper opcode;             /* Code for operation performed by this node */
+   int narg;                /* No. of arguments used by the operation */
+   struct UnitNode **arg;   /* Array of pointers to argument nodes */
+   double con;              /* Constant to be loaded by OP_LDCON operations */
+   struct KnownUnit *unit;  /* Known unit referred to by OP_LDVAR nodes */
+   Multiplier *mult;        /* Multiplier used by OP_LDVAR nodes */
+   const char *name;        /* User-defined unit referred to by OP_LDVAR
+                               nodes (no multiplier prefix included) */
+} UnitNode;
+
+/* A structure describing a known unit. */
+typedef struct KnownUnit {
+   const char *sym;         /* Unit symbol string (null terminated) */
+   const char *label;       /* Unit label string (null terminated) */
+   int symlen;              /* Length of symbol (without trailing null ) */
+   int lablen;              /* Length of label (without trailing null ) */
+   struct UnitNode *head;   /* Head of definition tree (NULL for basic units) */
+   struct KnownUnit *next;  /* Next KnownUnit in linked list */
+   struct KnownUnit *use;   /* KnownUnit to be used in place of this one */
+} KnownUnit;
+
+/* Module Variables. */
+/* ================= */
+
+/* A pointer to the KnownUnit structure at the head of a linked list of
+   such structures containing definitions of all known units. */
+static KnownUnit *known_units = NULL;
+
+/* An array of pointers to KnownUnits which list the basic quantities
+used in dimensional analysis. */
+static KnownUnit *quant_units[ NQUANT ];
+
+/* A pointer to the Multiplier structure at the head of a linked list of
+   such structures containing definitions of all known multipliers. */
+static Multiplier *multipliers = NULL;
+
+/* Set up mutexes */
+#ifdef THREAD_SAFE
+
+static pthread_mutex_t mutex1 = PTHREAD_MUTEX_INITIALIZER;
+#define LOCK_MUTEX1 pthread_mutex_lock( &mutex1 );
+#define UNLOCK_MUTEX1 pthread_mutex_unlock( &mutex1 );
+
+static pthread_mutex_t mutex2 = PTHREAD_MUTEX_INITIALIZER;
+#define LOCK_MUTEX2 pthread_mutex_lock( &mutex2 );
+#define UNLOCK_MUTEX2 pthread_mutex_unlock( &mutex2 );
+
+#else
+
+#define LOCK_MUTEX1
+#define UNLOCK_MUTEX1
+
+#define LOCK_MUTEX2
+#define UNLOCK_MUTEX2
+
+#endif
+
+/* Prototypes for Private Functions. */
+/* ================================= */
+static AstMapping *MakeMapping( UnitNode *, int * );
+static KnownUnit *GetKnownUnits( int, int * );
+static Multiplier *GetMultipliers( int * );
+static UnitNode *ConcatTree( UnitNode *, UnitNode *, int * );
+static UnitNode *CopyTree( UnitNode *, int * );
+static UnitNode *CreateTree( const char *, int, int, int * );
+static UnitNode *FixUnits( UnitNode *, UnitNode *, int * );
+static UnitNode *FreeTree( UnitNode *, int * );
+static UnitNode *MakeTree( const char *, int, int, int * );
+static UnitNode *MakeLabelTree( const char *, int, int * );
+static UnitNode *NewNode( UnitNode *, Oper, int * );
+static UnitNode *CombineFactors( UnitNode **, double *, int, double, int * );
+static const char *CleanExp( const char *, int * );
+static int EndsWith( const char *, int, const char *, int * );
+static int CmpTree( UnitNode *, UnitNode *, int, int * );
+static void FixConstants( UnitNode **, int, int * );
+static void InvertConstants( UnitNode **, int * );
+static UnitNode *InvertTree( UnitNode *, UnitNode *, int * );
+static void LocateUnits( UnitNode *, UnitNode ***, int *, int * );
+static void MakeKnownUnit( const char *, const char *, const char *, int * );
+static void MakeUnitAlias( const char *, const char *, int * );
+static void RemakeTree( UnitNode **, int * );
+static int SimplifyTree( UnitNode **, int, int * );
+static int ComplicateTree( UnitNode **, int * );
+static int ReplaceNode( UnitNode *, UnitNode *, UnitNode *, int * );
+static void FindFactors( UnitNode *, UnitNode ***, double **, int *, double *, int * );
+static const char *MakeExp( UnitNode *, int, int, int * );
+static int DimAnal( UnitNode *, double[NQUANT], double *, int * );
+static int Ustrcmp( const char *, const char *, int * );
+static int Ustrncmp( const char *, const char *, size_t, int * );
+static int SplitUnit( const char *, int, const char *, int, Multiplier **, int *, int * );
+static UnitNode *ModifyPrefix( UnitNode *, int * );
+static int ConStart( const char *, double *, int *, int * );
+
+/*  Debug functions...
+static const char *DisplayTree( UnitNode *, int );
+static void OpSym( UnitNode *, char * );
+static const char *OpName( Oper );
+static const char *TreeExp( UnitNode * );
+*/
+
+/* Function implementations. */
+/* ========================= */
+static const char *CleanExp( const char *exp, int *status ) {
+/*
+*  Name:
+*     CleanExp
+
+*  Purpose:
+*     Produce a clean copy of a units expression.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     const char *CleanExp( const char *exp )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function returns a pointer to dynamic memory containing a
+*     cleaned copy of the supplied units expression. Cleaning consists of
+*     the following operations:
+*        - removal of leading and trailing white space.
+*        - replacement of multiple adjacent spaces by a single space
+*        - removal of spaces adjacent to a parenthesis
+*        - removal of spaces adjacent to a binary operator
+*        - translates various common non-standard units into equivalent
+*          standard units.
+*
+*     Such carefull handling of spaces is necessary since a space is
+*     recognised by the MakeTree function as a multiplication operator.
+
+*  Parameters:
+*     exp
+*        A pointer to the expression to be cleaned.
+
+*  Returned Value:
+*     A pointer to the cleaned expression, which should be freed using
+*     astFree when no longer needed.
+
+*  Notes:
+*     - This function returns NULL if it is invoked with the global error
+*     status set, or if it should fail for any reason.
+*/
+
+/* Local Variables: */
+   char **tok;
+   char *p;
+   char *r;
+   char *result;
+   char *s;
+   char *t;
+   char *w;
+   const char *start;
+   int i;
+   int l;
+   int len;
+   char *tt;
+   int ntok;
+   int po;
+   int ps;
+   int word;
+
+/* Initialise */
+   result = NULL;
+
+/* Check inherited status */
+   if( !astOK ) return result;
+
+/* Split the supplied string up into tokens. Each block of contiguous
+   alphanumeric characters is a token. Each contiguous block of
+   non-alphanumerical characters is also a token. The + and - signs are
+   counted as alphanumeric. */
+   start = exp;
+   p = (char *) exp - 1;
+   word = ISWORD( *( p + 1 ) );
+   ntok = 0;
+   tok = NULL;
+   while( *(++p) ){
+      if( word ) {
+         if( !ISWORD( *p ) ) {
+            l = p - start;
+            t = astStore( NULL, start, l + 1 );
+            if( t ) t[ l ] = 0;
+            tok = astGrow( tok, ntok + 1, sizeof( char * ) );
+            if( tok ) tok[ ntok++ ] = t;
+            start = p;
+            word = 0;
+         }
+      } else {
+         if( ISWORD( *p ) ) {
+            l = p - start;
+            t = astStore( NULL, start, l + 1 );
+            if( t ) t[ l ] = 0;
+            tok = astGrow( tok, ntok + 1, sizeof( char * ) );
+            if( tok ) tok[ ntok++ ] = t;
+            start = p;
+            word = 1;
+         }
+      }
+   }
+
+   l = p - start;
+   t = astStore( NULL, start, l + 1 );
+   if( t ) t[ l ] = 0;
+   tok = astGrow( tok, ntok + 1, sizeof( char * ) );
+   if( tok ) tok[ ntok++ ] = t;
+
+/* Check the tokens for known non-standard unit syntax, and replace with the
+   equivalent standard syntax. Starlink SPLAT has a class called UnitUtilities
+   which has more of these common units mistakes. AST has to be a bit
+   more conservative than SPLAT though because of its wider remit. */
+   len = 0;
+   tt = NULL;
+   for( i = 0; i < ntok; i++ ) {
+      t = tok[ i ];
+      l = strlen( t );
+      tt = astStore( tt, t, l + 1 );
+
+/* Any alphabetical word followed by a digit is taken as <word>^<digit>.
+   Any alphabetical word followed by a sign and a digit is taken as
+   <word>^<sign><digit>. */
+      if( l > 1 && *t != '-' && *t != '+' &&
+          strcspn( t, "0123456789" ) == l - 1 ) {
+         tok[ i ] = astMalloc( l + 2 );
+         if( tok[ i ] ) {
+            strcpy( tok[ i ], t );
+            w = t + l - 2;
+            if( *w != '+' && *w != '-' ) {
+               tok[ i ][ l - 1 ] = '^';
+               strcpy( tok[ i ] + l, t + l - 1 );
+            } else {
+               tok[ i ][ l - 2 ] = '^';
+               strcpy( tok[ i ] + l - 1, t + l - 2 );
+            }
+            t = astFree( t );
+         }
+         l++;
+
+/* If the word ends with "micron" change to "(<start>m*1.0E-6)". Should be OK
+   for things like "Kmicron". */
+      } else if( ( s = strstr( t, "micron" ) ) ) {
+         tok[ i ] = astMalloc( s - t + 11 );
+         if( tok[ i ] ) {
+            w = tok[ i ];
+            *(w++) = '(';
+            if( s > t ) {
+               strncpy( w, t, s - t );
+               w += s - t;
+            }
+            strcpy( w, "m*1.0E-6)" );
+            l = s - t + 11;
+            t = astFree( t );
+         }
+
+/* Convert "STER" to "sr". */
+      } else if( !Ustrcmp( t, "STER", status ) ) {
+         tok[ i ] = astStore( NULL, "sr", 3 );
+         l = 2;
+         t = astFree( t );
+
+/* If the word ends with "JY" and is preceded by a single character, change
+   to "<start>Jy". Should be OK for things like "MJY". */
+      } else if( l == 3 && !strcmp( t + 1, "JY" ) ) {
+         tok[ i ][ 2 ] = 'y';
+
+/* If the word begins with "nano" (case-insensitive) change "nano" to
+   "n". Such changes are usually handled by SplitUnit, but we need to
+   handle this as a special case here since scanf seems to read "nan" as
+   a string representation of NaN. */
+      } else if( !Ustrncmp( t, "nano", 4, status ) ) {
+         tok[ i ] = astStore( NULL, t + 3, l - 2 );
+         if( tok[ i ] ) {
+            *(tok[ i ]) = 'n';
+            t = astFree( t );
+         }
+         l -= 3;
+      }
+
+/* Update the total length of the string. */
+      len += l;
+   }
+   tt = astFree( tt );
+
+/* Concatentate the tokens into a single string, freeing the individual
+   strings. */
+   result = astMalloc( len + 1 );
+   if( result ) {
+      p = result;
+      for( i = 0; i < ntok; i++ ) {
+         len = strlen( tok[ i ] );
+         memcpy( p, tok[ i ], len );
+         p += len;
+         tok[ i ] = astFree( tok[ i ] );
+      }
+      *p = 0;
+      tok = astFree( tok );
+
+/* Now do other cleaning.
+   ---------------------- */
+
+/* Initialise a pointer to the previous character read from the string. */
+      r = result - 1;
+
+/* Initialise a pointer to the next character to be written to the string. */
+      w = result;
+
+/* Pretend the previous character written to the string was a space. */
+      ps = 1;
+
+/* Read all the supplied string, copying it to earlier parts of the
+   string discarding leading spaces and multiple adjacent embedded spaces in
+   the process. */
+      while( *(++r) ) {
+
+/* If the character read is a space, only write it to the string if the
+   previous character written was not a space (in which case set a flag
+   to indicate that the previous character written to the string is now a
+   space). */
+         if( isspace( *r ) ) {
+            if( !ps ) {
+               *(w++) = *r;
+               ps = 1;
+            }
+
+/* Write all non-space characters to the string, and clear the flag which
+   indicates if the previous character written to the string was a space. */
+         } else {
+            *(w++) = *r;
+            ps = 0;
+         }
+      }
+
+/* If the last character written to the string was a space, reduce the
+   length of the string by one since we do not want any trailing spaces. */
+      if( ps ) w--;
+
+/* Terminate the string. */
+      *w = 0;
+
+/* We now need to pass through the string again, this time removing any
+   spaces which are adjacent to a binary operator or a parenthesis. */
+      r = result - 1;
+      w = result;
+      ps = 0;
+      po = 0;
+      while( *(++r) ) {
+
+/* If the current character is a space, only write it if the previous
+   written character was not an operator or parenthesis. */
+         if( isspace( *r ) ) {
+            if( !po ) {
+               *(w++) = *r;
+               po = 1;
+               ps = 1;
+            }
+
+/* If the current character is an operator or parenthesis, back up one
+   character before writing it out if the previous written character was
+   a space. */
+         } else if( *r == '*' || *r == '/' || *r == '^' || *r == '.' ||
+                    *r == ')' || *r == '(' ) {
+            if( ps ) w--;
+            *(w++) = *r;
+            po = 1;
+            ps = 0;
+
+/* If the current character is not a space and not an operator symbol,
+   just write it out. */
+         } else {
+            *(w++) = *r;
+            po = 0;
+            ps = 0;
+         }
+      }
+
+/* Terminate the string. */
+      if( ps ) w--;
+      *w = 0;
+
+   }
+
+/* Return the result. */
+   return (const char *) result;
+}
+
+static int CmpTree( UnitNode *tree1, UnitNode *tree2, int exact, int *status ) {
+/*
+*  Name:
+*     CmpTree
+
+*  Purpose:
+*     Compares two trees of UnitNodes.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     int CmpTree( UnitNode *tree1, UnitNode *tree2, int exact, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function returns a zero value if the two trees are
+*     equivalent. This requires the trees to have identical structure
+*     except that, if "exact" is zero, arguments for OP_MULT nodes can
+*     be swapped.
+*
+*     If the trees are not equivalent then a value of +1 or -1 is returned
+*     depending on whether tree1 should be placed before or after tree2
+*     in a sorted list of trees.
+
+*  Parameters:
+*     tree1
+*        A pointer to the UnitNode at the head of the first tree.
+*     tree2
+*        A pointer to the UnitNode at the head of the second tree.
+*     exact
+*        If non-zero, then OP_MULT nodes must have their arguments the
+*        same way round in order for the OP_MULT nodes to match. Otherwise,
+*        OP_MULT nodes with equivalent arguments match even if the
+*        arguments are swapped.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     Zero if the two trees are equal. +1 if tree1 should be placed before
+*     tree2 in a sorted list of trees. -1 if tree1 should be placed after
+*     tree2 in a sorted list of trees.
+
+* Notes:
+*     - Zero is returned if an error has already occurred, or
+*     if this function fails for any reason.
+
+*/
+
+/* Local Variables: */
+   int result;
+   int i;
+   Oper op;
+
+/* Initialise. */
+   result = 0;
+
+/* Check inherited status. */
+   if( !astOK ) return result;
+
+/* If the op codes differ, compare them as integers. */
+   op = tree1->opcode;
+   if( op != tree2->opcode ) {
+      result = ( op > tree2->opcode ) ? 1: -1;
+
+/* If both supplied nodes are OP_LDVAR nodes, compare the associated names. */
+   } else if( op == OP_LDVAR ){
+      result = strcmp( tree1->name, tree2->name );
+
+/* If both supplied nodes are constant nodes, compare the constant values. */
+   } else if( tree1->con != AST__BAD ){
+      result = astEQUAL( tree1->con, tree2->con ) ? 0 : (
+                 ( tree1->con > tree2->con ) ? 1 : -1 );
+
+/* Otherwise, compare the arguments for the node. */
+   } else {
+      for( i = 0; i < tree1->narg; i++ ) {
+         result = CmpTree( tree1->arg[ i ], tree2->arg[ i ], exact, status );
+         if( result ) break;
+      }
+
+/* If the head nodes of the two trees are OP_MULT nodes, and the above
+   check determined they are different, this may be just because they
+   have their operands swapped. If "exact" si zero, this is considered an
+   insignificant difference between the two trees which we should ignore.
+   To check for this try comparing the arguments again, this time swapping
+   the arguments of tree2. */
+      if( result && op == OP_MULT && !exact ) {
+         for( i = 0; i < tree1->narg; i++ ) {
+            result = CmpTree( tree1->arg[ i ], tree2->arg[ 1 - i ], 0, status );
+            if( result ) break;
+         }
+      }
+   }
+
+/* If an error has occurred, return zero. */
+   if( !astOK ) result = 0;
+
+/* Return the answer. */
+   return result;
+}
+
+static UnitNode *CombineFactors( UnitNode **factors, double *powers,
+                                 int nfactor, double coeff, int *status ) {
+/*
+*  Name:
+*     CombineFactors
+
+*  Purpose:
+*     Create a tree which represents the product of the supplied factors.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     UnitNode *CombineFactors( UnitNode **factors, double *powers,
+*                               int nfactor, double coeff, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function createa a tree of UnitNodes which represents the
+*     product of the supplied factors, and the supplied coefficient.
+*     The factors are sorted before being combined, using the sort order
+*     implemented by the CmpTree function.
+
+*  Parameters:
+*     factors
+*        A pointer to an array with "nfactor" elements, each element being
+*        a pointer to a UnitNode which is a factor of the required tree.
+*        On exit, the array is sorted.
+*     powers
+*        A pointer to an array with "nfactor" elements, each element being a
+*        double holding the power of the associated factor in "factors".
+*        On exit, the array reflects the sorting applied to "factors".
+*     nfactor
+*        The number of elements in the "factors" and "powers" arrays.
+*     coeff
+*        The overall coefficient to be applied to the product of the factors.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     A pointer to a UnitNode which is at the head of the new tree.
+
+* Notes:
+*     - A NULL pointer is returned if an error has already occurred, or
+*     if this function fails for any reason.
+
+*/
+
+/* Local Variables: */
+   UnitNode *result;
+   int i;
+   int j;
+   int jp;
+   int done;
+   UnitNode *ftmp;
+   UnitNode *node1;
+   UnitNode *node2;
+   UnitNode *pnode;
+   double ptmp;
+
+/* Initialise. */
+   result = NULL;
+
+/* Check inherited status. */
+   if( !astOK ) return result;
+
+/* Sort the supplied list of factors, modifying the powers array
+   correspondingly. A simple bubblesort algorithm is used since there
+   will only be a handfull of factors. */
+   for( i = nfactor - 1; i > 0; i-- ) {
+      done = 1;
+      for( j = 0, jp = 1; j < i; j++, jp++ ) {
+         if( CmpTree( factors[ j ], factors[ jp ], 0, status ) > 0 ) {
+            ftmp = factors[ j ];
+            factors[ j ] = factors[ jp ];
+            factors[ jp ] = ftmp;
+
+            ptmp = powers[ j ];
+            powers[ j ] = powers[ jp ];
+            powers[ jp ] = ptmp;
+
+            done = 0;
+         }
+      }
+      if( done ) break;
+   }
+
+/* The first root term of the returned tree is the coefficient, unless the
+   coefficient is 1.0, in which case it will be the first factor. */
+   if( coeff != 1.0 ) {
+      node1 = NewNode( NULL, OP_LDCON, status );
+      if( astOK ) node1->con = coeff;
+   } else {
+      node1 = NULL;
+   }
+
+/* Loop through the factors. */
+   for( i = 0; i < nfactor; i++ ) {
+
+/* If the power of this factor is zero, we ignore the factor. */
+      if( powers[ i ] != 0.0 ) {
+
+/* If the power of this factor is one, we use the factor directly. */
+         if( astEQUAL( powers[ i ], 1.0 ) ) {
+            node2 = CopyTree( factors[ i ], status );
+
+/* Otherwise, for non-zero, non-unity powers, we create a POW node for
+   the factor. */
+         } else {
+            node2 = NewNode( NULL, OP_POW, status );
+            pnode = NewNode( NULL, OP_LDCON, status );
+            if( astOK ) {
+               pnode->con = powers[ i ];
+               node2->arg[ 0 ] = CopyTree( factors[ i ], status );
+               node2->arg[ 1 ] = pnode;
+            }
+         }
+
+/* We now combine node1 and node2 using an OP_MULT node, which becomes
+   the "node1" for the next pass. On the first pass we may have no node1 (if
+   the supplied coefficient was 1.0), in which case we reserve the current
+   node2 as the node1 for the next pass. */
+         if( node1 ) {
+            result = NewNode( NULL, OP_MULT, status );
+            if( astOK ) {
+               result->arg[ 0 ] = node1;
+               result->arg[ 1 ] = node2;
+               node1 = result;
+            }
+         } else {
+            node1 = node2;
+         }
+      }
+   }
+
+/* Ensure we have a node to return. */
+   if( astOK ) {
+      if( !result ) result = node1;
+      if( !result ) {
+         result = NewNode( NULL, OP_LDCON, status );
+         if( astOK ) result->con = 1.0;
+      }
+   }
+
+/* If an error has occurred, free any new tree. */
+   if( !astOK ) result = FreeTree( result, status );
+
+/* Return the answer. */
+   return result;
+}
+
+static int ComplicateTree( UnitNode **node, int *status ) {
+/*
+*  Name:
+*     ComplicateTree
+
+*  Purpose:
+*     Removes standardisations introduced by SimplifyTree.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     int ComplicateTree( UnitNode **node )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function modifies a tree of UnitNodes by removing standardisations
+*     introduced by SimplifyTree. The standardisations removed are ones
+*     which would make the corresponding algebraic expression (as produced
+*     by MakeExp) unnatural to a human reader.
+
+*  Parameters:
+*     node
+*        The address of a pointer to the UnitNode at the head of the tree
+*        which is to be complicated. On exit the supplied tree is freed and
+*        a pointer to a new tree is placed at the given address.
+
+*  Returned Value:
+*     Non-zero if any change was introduced into the tree.
+
+*/
+
+/* Local Variables: */
+   int i;
+   UnitNode *newnode;
+   UnitNode *node1;
+   UnitNode *node2;
+   UnitNode *node3;
+   Oper op;
+   double con;
+   double fk;
+   int k;
+   int result;
+   double kcon;
+
+/* Initialise */
+   result = 0;
+
+/* Check inherited status. */
+   if( !astOK ) return result;
+
+/* Initiallially, we have no replacement node. */
+   newnode = NULL;
+   node1 = NULL;
+   node3 = NULL;
+
+/* Complicate the sub-trees corresponding to the arguments of the node at
+   the head of the supplied tree. */
+   for( i = 0; i < (*node)->narg; i++ ) {
+      if( ComplicateTree( &( (*node)->arg[ i ] ), status ) ) result = 1;
+   }
+
+/* Now undo specific simplifications appropriate to the nature of the node at
+   the head of the tree. */
+   op = (*node)->opcode;
+
+/* If the head is an OP_MULT node with a constant first argument and
+   a "LN" second argument, rearrange the nodes to represent ln(x**k) instead
+   of k*ln(x). If k is an integer multiple of "0.1/ln(10)" convert the "ln"
+   function into a "log" (base 10) function. Check for "k==1" in which
+   case we do not need a POW node. */
+   if( (*node)->opcode == OP_MULT ) {
+
+      con = (*node)->arg[ 0 ]->con;
+      if( con != AST__BAD && (*node)->arg[ 1 ]->opcode == OP_LN ) {
+         fk = 10.0*con*log( 10.0 );
+         k = NINT(fk);
+         if( astEQUAL(fk,((double)k)) ) {
+            newnode = NewNode( NULL, OP_LOG, status );
+            con = k/10.0;
+         } else {
+            newnode = NewNode( NULL, OP_LN, status );
+         }
+
+         node2 = CopyTree( (*node)->arg[ 1 ]->arg[ 0 ], status );
+         if( !astEQUAL( con, 1.0 ) ){
+            node1 = CopyTree( (*node)->arg[ 0 ], status );
+            node3 = NewNode( NULL, OP_POW, status );
+         }
+
+         if( astOK ) {
+            if( !astEQUAL( con, 1.0 ) ){
+               node1->con = con;
+               node3->arg[ 0 ] = node2;
+               node3->arg[ 1 ] = node1;
+               newnode->arg[ 0 ] = node3;
+            } else {
+               newnode->arg[ 0 ] = node2;
+            }
+         }
+
+/* Replace "(A**-1)*B" with "B/A" */
+      } else if( (*node)->arg[ 0 ]->opcode == OP_POW &&
+                 astEQUAL( (*node)->arg[ 0 ]->arg[ 1 ]->con, -1.0 )) {
+         newnode = NewNode( NULL, OP_DIV, status );
+         if( astOK ) {
+            newnode->arg[ 0 ] = CopyTree( (*node)->arg[ 1 ], status );
+            newnode->arg[ 1 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
+         }
+
+/* Replace "B*(A**-1)" with "B/A" */
+      } else if( (*node)->arg[ 1 ]->opcode == OP_POW &&
+                 astEQUAL( (*node)->arg[ 1 ]->arg[ 1 ]->con, -1.0 )) {
+         newnode = NewNode( NULL, OP_DIV, status );
+         if( astOK ) {
+            newnode->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ], status );
+            newnode->arg[ 1 ] = CopyTree( (*node)->arg[ 1 ]->arg[ 0 ], status );
+         }
+
+/* Convert (x**k)*(y**k) to (x*y)**k. */
+      } else if( (*node)->arg[ 0 ]->opcode == OP_POW &&
+                 (*node)->arg[ 1 ]->opcode == OP_POW &&
+                 astEQUAL( (*node)->arg[ 0 ]->arg[ 1 ]->con,
+                        (*node)->arg[ 1 ]->arg[ 1 ]->con )) {
+         newnode = NewNode( NULL, OP_POW, status );
+         node1 = NewNode( NULL, OP_MULT, status );
+         if( astOK ) {
+            node1->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
+            node1->arg[ 1 ] = CopyTree( (*node)->arg[ 1 ]->arg[ 0 ], status );
+            newnode->arg[ 0 ] = node1;
+            newnode->arg[ 1 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 1 ], status );
+         }
+
+/* Convert c*sqrt(x) to sqrt((c**2)*x) (if c > 0). */
+      } else if( (kcon=(*node)->arg[ 0 ]->con) != AST__BAD &&
+                 kcon > 0.0 && (*node)->arg[ 1 ]->opcode == OP_SQRT ) {
+         newnode = NewNode( NULL, OP_SQRT, status );
+         node1 = NewNode( NULL, OP_MULT, status );
+         node2 = NewNode( NULL, OP_LDCON, status );
+         if( astOK ) {
+            node2->con = kcon*kcon;
+            node1->arg[ 0 ] = node2;
+            node1->arg[ 1 ] = CopyTree( (*node)->arg[ 1 ]->arg[ 0 ], status );
+            newnode->arg[ 0 ] = node1;
+         }
+      }
+
+/* If the head node is a POW node, replace "x**0.5" by sqrt(x) */
+   } else if( (*node)->opcode == OP_POW ) {
+      if( astEQUAL( (*node)->arg[ 1 ]->con, 0.5 ) ) {
+         newnode = NewNode( NULL, OP_SQRT, status );
+         if( astOK ) {
+            newnode->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ], status );
+         }
+      }
+   }
+
+/* If we have produced a new node which is identical to the old node,
+   free it. Otherwise, indicate we have made some changes. */
+   if( newnode ) {
+      if( !CmpTree( newnode, *node, 1, status ) ) {
+         newnode = FreeTree( newnode, status );
+      } else {
+         result = 1;
+      }
+   }
+
+/* If an error has occurred, free any new node. */
+   if( !astOK ) {
+      newnode = FreeTree( newnode, status );
+      result = 0;
+   }
+
+/* If we have a replacement node, free the supplied tree and return a
+   pointer to the new tree. */
+   if( newnode ) {
+      FreeTree( *node, status );
+      *node = newnode;
+   }
+
+/* If the above produced some change, try simplifying (without
+   re-introducing the standardisation we have just got rid of!) and
+   then re-complicating the tree. */
+   if( result ) {
+      SimplifyTree( node, 0, status );
+      ComplicateTree( node, status );
+   }
+
+/* Return the result. */
+   return result;
+}
+
+static UnitNode *ConcatTree( UnitNode *tree1, UnitNode *tree2, int *status ) {
+/*
+*  Name:
+*     ConcatTree
+
+*  Purpose:
+*     Concatenate two trees together.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     UnitNode *ConcatTree( UnitNode *tree1, UnitNode *tree2, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function a pointer to the head of a new tree of UnitNodes which
+*     is formed by feeding the output of "tree1" (i.e. the quantity
+*     represented by the node at the head of tree1) into the (single)
+*     input of "tree2" (i.e. the single OP_LDVAR Node containined within
+*     tree2).
+
+*  Parameters:
+*     tree1
+*        A pointer to the first tree.
+*     tree2
+*        A pointer to the second tree. This should have no more than one
+*        OP_LDVAR node.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     A pointer to a UnitNode which is at the head of the new tree.
+
+* Notes:
+*     - If "tree2" contains zero units, a NULL pointer is returned but no
+*     error is reported.
+*     - If "tree2" contains more than one unit, an error is reported
+*     error is reported.
+*     - A NULL pointer is returned if an error has already occurred, or
+*     if this function fails for any reason.
+
+*/
+
+/* Local Variables: */
+   UnitNode *result;
+   UnitNode **units;
+   int nunits;
+
+/* Initialise. */
+   result = NULL;
+
+/* Check inherited status. */
+   if( !astOK ) return result;
+
+/* Produce a copy of tree2. */
+   result = CopyTree( tree2, status );
+
+/* Locate the OP_LDVAR node in the copy of tree2. */
+   units = NULL;
+   nunits = 0;
+   LocateUnits( result, &units, &nunits, status );
+
+/* If no OP_LDVAR nodes were found in tree2, we cannot concatenate the
+   trees. */
+   if( nunits > 0 ) {
+
+/* Report an error if the number of pointers returned is larger than 1. */
+      if( nunits > 1 && astOK ) {
+         astError( AST__INTER, "ConcatTree(unit): tree2 uses %d units - "
+                   "should be 1 (internal AST programming error).", status, nunits );
+      }
+
+/* Replace the OP_LDVAR node in the copy of tree2 with a copy of tree1. */
+      if( astOK ) {
+
+/* If the node at the head of the supplied tree2 is the node to be
+   replaced, just free the tree created earlier and return a copy of
+   tree1. */
+         if( units[ 0 ] == result ) {
+            FreeTree( result, status );
+            result = CopyTree( tree1, status );
+
+/* Otherwise, search for the node to be replaced and do the substitution
+   within the tree created earlier. */
+         } else {
+            ReplaceNode( result, units[ 0 ], CopyTree( tree1, status ), status );
+         }
+      }
+   }
+
+/* Free resources. */
+   units = astFree( units );
+
+/* If an error has occurred, free any new tree. */
+   if( !astOK ) result = FreeTree( result, status );
+
+/* Return the answer. */
+   return result;
+}
+
+static int ConStart( const char *text, double *val, int *nc, int *status ) {
+/*
+*  Name:
+*     ConStart
+
+*  Purpose:
+*     See if the supplied string starts with a literal numeric constant.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     int ConStart( const char *text, double *val, int *nc, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function checks if the supplied string starts with a literal
+*     numeric constant and returns it if it does. It is a wrap-up for scanf
+*     since scanf has non-standard behaviour on some platforms (e.g. Cygwin
+*     scanf interprets the character "n" as a floating point number!).
+
+*  Parameters:
+*     text
+*        The text to check.
+*     val
+*        Address of a double to receive any numerical constant read
+*        from the start of the string. Unity is returned if the string
+*        does not start with a numerical constant.
+*     nc
+*        Address of an int to receive the number of characters used to
+*        create the value returned in "val". Zero is returned if the
+*        string does not start with a numerical constant.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     Non-zero if the text started with a numerical constant.
+
+*/
+
+/* Local Variables: */
+   int result;
+   const char *c;
+
+/* Initialise */
+   *nc = 0;
+   *val = 1.0;
+
+/* Return zero if no text was supplied */
+   if( !text ) return 0;
+
+/* Use sscanf to see if the string begin with a numerical constant */
+   result = astSscanf( text, "%lf%n", val, nc );
+
+/* If so, check that the first non-blank character in the string
+   is not "N" (interpreted by Cygwin as numerical zero!). */
+   if( result ) {
+      c = text;
+      while( isspace( *c ) ) c++;
+      if( *c == 'n' || *c == 'N' ) {
+         result = 0;
+         *nc = 0;
+         *val = 1.0;
+      }
+   }
+
+/* Return the result. */
+   return result;
+}
+
+static UnitNode *CopyTree( UnitNode *tree, int *status ) {
+/*
+*  Name:
+*     CopyTree
+
+*  Purpose:
+*     Create a new tree of UnitNodes containing a copy of a given tree.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     UnitNode *CopyTree( UnitNode *tree, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function creates a copy of the supplied tree of UnitNodes.
+
+*  Parameters:
+*     tree
+*        The UnitNode at the head of the tree to be copied.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     A pointer to the UnitNode at the head of the new tree.
+
+*  Notes:
+*     - A value of NULL will be returned if this function is invoked with
+*     the global error status set, or if it should fail for any reason.
+*/
+
+/* Local Variables: */
+   UnitNode **args;
+   UnitNode *result;
+   int i;
+   int narg;
+
+/* Initialise. */
+   result = NULL;
+
+/* Check the inherited status. */
+   if( !astOK || !tree ) return result;
+
+/* Create a new node to represent the head of the supplied tree. */
+   result = astMalloc( sizeof( UnitNode ) );
+
+/* Check pointers can be used safely. */
+   if( astOK ) {
+
+/* Copy the fields of the supplied node. */
+      narg = tree->narg;
+
+      result->arg = NULL;
+      result->unit = tree->unit;
+      result->mult = tree->mult;
+      result->opcode = tree->opcode;
+      result->narg = narg;
+      result->con = tree->con;
+      result->name = tree->name ? astStore( NULL, tree->name,
+                                            strlen( tree->name ) + 1 ) : NULL;
+
+/* Create an array of UnitNode pointers for the arguments. */
+      args = astMalloc( narg*sizeof( UnitNode * ) );
+      if( astOK ) {
+         result->arg = args;
+
+/* Copy the sub-trees headed by the argument nodes. */
+         for( i = 0; i < narg; i++ ) {
+            args[ i ] = CopyTree( tree->arg[ i ], status );
+         }
+      }
+   }
+
+/* Free any result if an error occurred. */
+   if( !astOK ) result = FreeTree( result, status );
+
+/* Return the answer. */
+   return result;
+}
+
+static UnitNode *CreateTree( const char *exp, int basic, int lock, int *status ){
+/*
+*  Name:
+*     CreateTree
+
+*  Purpose:
+*     Convert an algebraic units expression into a tree of UnitNodes.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     UnitNode *CreateTree( const char *exp, int basic, int lock, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function converts the supplied algebraic units expression into
+*     a tree of UnitNodes. The result tree can optionally be expanded to
+*     create a tree in which the "roots" (LDVAR nodes) all refer to
+*     basic units.
+
+*  Parameters:
+*     exp
+*        The units expression. This should not include any leading or
+*        trailing spaces.
+*     basic
+*        Should the tree created from parsing "exp" be expanded so that
+*        the leaf nodes of the tree are all basic units?
+*     lock
+*        Use a mutex to guard access to the KnownUnits list?
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     A pointer to a UnitNode which forms the head of a tree of UnitNodes
+*     representing the supplied unit expression.
+
+*  Notes:
+*     -  A NULL value is returned if this function is invoked with the
+*     global error status set or if it should fail for any reason.
+*/
+
+/* Local Variables: */
+   UnitNode *result;
+   const char *cleanex;
+
+/* Initialise */
+   result = NULL;
+
+/* Check the global error status, and that we have a string. */
+   if ( !astOK ) return result;
+
+/* Produce a clean copy of the supplied string. This has no leading
+   or trailing white space, and any spaces adjacent to operators within
+   the string are removed (this is needed because spaces are treated as
+   multiplication symbols). */
+   cleanex = CleanExp( exp, status );
+
+/* If the string is blank, return the NULL pointer. Otherwise, create a
+   tree of UnitNodes describing the units. The returned tree has LDVAR
+   nodes which refer to the unit symbols contained in the supplied string. */
+   if( cleanex && (*cleanex) ) {
+      result = MakeTree( cleanex, strlen( cleanex ), lock, status );
+
+/* Replace each subtree which simply combines constants (i.e. which has no
+   OP_LDVAR nodes) with a single OP_LDCON node. */
+      FixConstants( &result, 0, status );
+
+/* Invert literal constant unit multipliers. */
+      InvertConstants( &result, status );
+
+/* Now replace each LDVAR node which refers to a known derived unit with
+   a sub-tree which defines the derived unit in terms of known basic units.
+   The LDVAR nodes in the resulting tree all refer to basic units. */
+      if( basic ) RemakeTree( &result, status );
+   }
+
+/* Free resources. */
+   cleanex = astFree( (void *) cleanex );
+
+/* Free any returned tree if an error has occurred. */
+   if( !astOK ) result = FreeTree( result, status );
+
+/* Return the result. */
+   return result;
+}
+
+static int DimAnal( UnitNode *node, double powers[NQUANT], double *scale, int *status ) {
+/*
+*  Name:
+*     DimAnal
+
+*  Purpose:
+*     Perform a dimensional analysis of a unit tree.
+
+*  Type:
+*     Protected function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     int DimAnal( UnitNode *node, double powers[NQUANT], double *scale, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function returns a set of powers and a scaling factor which
+*     represent the units tree.
+
+*  Parameters:
+*     node
+*        Pointer to the UnitNode at the head of the unit tree.
+*     powers
+*        An array in which are returned the powers for each of the following
+*        basic units (in the order shown): kilogramme, metre, second, radian,
+*        Kelvin, count, adu, photon, magnitude, pixel. If the supplied unit
+*        does not depend on a given basic unit a value of 0.0 will be returned
+*        in the array. The returns values represent a system of units which is
+*        a scaled form of the supplied units, expressed in the basic units of
+*        m, kg, s, rad, K, count, adu, photon, mag and pixel. For instance, a
+*        returned array of [1,0,-2,0,0,0,0,0,0,0] would represent "m/s**2".
+*     scale
+*        Pointer to a location at which to return a scaling factor for the
+*        supplied units. The is the value, in the units represented by the
+*        returned powers, which corresponds to a value of 1.0 in the supplied
+*        units.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     Non-zero if the tree was analysed succesfully. Zero otherwise.
+
+*  Notes:
+*     -  Zero is returned if this function is invoked with the
+*     global error status set or if it should fail for any reason.
+*/
+
+/* Local Variables; */
+   Oper oper;
+   int result;
+   int i;
+   double p0[ NQUANT ];
+   double p1[ NQUANT ];
+   double s0;
+   double s1;
+
+/* Check inherited status */
+   if( !astOK ) return 0;
+
+/* Initialise the powers of all dimensions to zero, and set the scaling
+   factor to unity. */
+   result = 1;
+   *scale = 1.0;
+   for( i = 0; i < NQUANT; i++ ) powers[ i ] = 0.0;
+
+/* Load constant: constant is dimensionaless so leave powers unchanged,
+   and set the scaling factor. */
+   oper = node->opcode;
+   if( oper == OP_LDCON ) {
+      *scale = 1.0/node->con;
+
+/* Load variable: check it is one of the basic known dimensional
+   quantities. If so, set the power of the quantity to unity and store
+   the scale factor. If the unit is "g" modify the scale factor so that
+   the analysis quantity is "kg". */
+   } else if( oper == OP_LDVAR ) {
+      result = 0;
+      for( i = 0; i < NQUANT; i++ ) {
+         if( node->unit == quant_units[ i ] ) {
+            powers[ i ] = 1.0;
+            *scale = node->mult ? 1.0/node->mult->scale : 1.0;
+            if( !strcmp( node->unit->sym, "g" ) ) *scale *= 0.001;
+            result = 1;
+            break;
+         }
+      }
+
+/* How does dimensional analysis handle log or exp units?*/
+   } else if( oper == OP_LOG ) {
+      result= 0;
+
+   } else if( oper == OP_LN ) {
+      result= 0;
+
+   } else if( oper == OP_EXP ) {
+      result= 0;
+
+/* Get the powers for the child unit and then multiply each by 0.5 and
+   take the square root of the scale factor. */
+   } else if( oper == OP_SQRT ) {
+      result = DimAnal( node->arg[0], powers, scale, status );
+      if( result ) {
+         for( i = 0; i < NQUANT; i++ ) powers[ i ]*= 0.5;
+         *scale = sqrt( *scale );
+      }
+
+/* Similarly for pow nodes. */
+   } else if( oper == OP_POW ) {
+      result = DimAnal( node->arg[0], powers, scale, status );
+      if( result ) {
+         double power = node->arg[1]->con;
+         for( i = 0; i < NQUANT; i++ ) powers[ i ]*= power;
+         *scale = pow( *scale, power );
+      }
+
+/* Binary operators. Analyses the operands dimensions and combine. */
+   } else if( oper == OP_DIV ) {
+      if( DimAnal( node->arg[0], p0, &s0, status ) &&
+          DimAnal( node->arg[1], p1, &s1, status ) ) {
+         for( i = 0; i < NQUANT; i++ ) powers[ i ] = p0[ i ] - p1[ i ];
+         *scale = s0/s1;
+      } else {
+         result = 0;
+      }
+
+   } else if( oper == OP_MULT ) {
+      if( DimAnal( node->arg[0], p0, &s0, status ) &&
+          DimAnal( node->arg[1], p1, &s1, status ) ) {
+         for( i = 0; i < NQUANT; i++ ) powers[ i ] = p0[ i ] + p1[ i ];
+         *scale = s0*s1;
+      } else {
+         result = 0;
+      }
+
+/* Named constants are dimensionless */
+   } else if( oper == OP_LDPI ) {
+      *scale = 1.0/PI;
+
+   } else if( oper == OP_LDE ) {
+      *scale = 1.0/E;
+
+   }
+
+   return result;
+
+}
+
+static int EndsWith( const char *c, int nc, const char *test, int *status ){
+/*
+*  Name:
+*     EndsWith
+
+*  Purpose:
+*     See if a string ends with another string
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     int EndsWith( const char *c, int nc, const char *test, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function sees if the string given by "c" ends with the string
+*     given by "test". The comparison is case-insensitive.
+
+*  Parameters:
+*     c
+*        A pointer to the last character in the string to be tested.
+*     nc
+*        The number of characters in the string to be tested.
+*     test
+*        A pointer to the string to be tested for.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     Non-zero if the string "c" ends with the string "test".
+
+*/
+
+/* Local Variables: */
+   const char *start;
+   int i;
+   int result;
+   int tlen;
+
+/* initialise. */
+   result = 0;
+
+/* Check inherited status. */
+   if( !astOK ) return result;
+
+/* Check the string being tested for is not longer than the string being
+   tested. */
+   tlen = strlen( test );
+   if( tlen <= nc ){
+
+/* Get a pointer to where the matching string would start if the string "c"
+   ends with the required string "test". */
+      start = c - tlen + 1;
+
+/* Do the comparison. */
+      result = 1;
+      for( i = 0; i < tlen; i++ ) {
+         if( tolower( start[ i ] ) != tolower( test[ i ] ) ) {
+            result = 0;
+            break;
+         }
+      }
+   }
+
+/* Return the result. */
+   return result;
+
+}
+
+static void FindFactors( UnitNode *node, UnitNode ***factors, double **powers,
+                         int *nfactor, double *coeff, int *status ){
+/*
+*  Name:
+*     FindFactors
+
+*  Purpose:
+*     Find the factors within an expression given by a tree of UnitNodes.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     void FindFactors( UnitNode *node, UnitNode ***factors, double **powers,
+*                       int *nfactor, double *coeff, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function analyses the supplied tree of UnitNoes and returns
+*     an array of pointers to nodes within the supplied tree which form
+*     factors of the tree. The power associated with each factor is also
+*     returned, together with an overall coefficient for the tree. The
+*     expression represented by the tree is thus the product of the
+*     coefficient with each of the factors, each raised to the associated
+*     power.
+
+*  Parameters:
+*     node
+*        A pointer to the UnitNode at the head of the tree which is to be
+*        analysed.
+*     factors
+*        The address at which to return a pointer to an array with "*nfactor"
+*        elements, each element being a pointer to a UnitNode within the
+*        supplied tree which is a factor of the supplied tree.
+*     powers
+*        The address at which to return a pointer to an array with "*nfactor"
+*        elements, each element being a double holding the power of the
+*        associated factor in "*factors".
+*     nfactor
+*        The address of an int containing the number of elements in the
+*        returned "*factors" and "*powers" arrays.
+*     coeff
+*        The address of a double containing the overall coefficient to be
+*        applied to the product of the factors.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Notes:
+*     - If the supplied node is a constant node, then "*coeff" is
+*     returned holding the value of the constant and "*nfactor" is returned
+*     equal to zero ("*factors" and "*powers" are returned holding NULL).
+*     - If an error has already occurred, or if this function fails, then
+*     "*factors" and "*powers" are returned holding NULL, "*nfactor" is
+*     returned holding zero and "*coeff" is returned holding 1.0.
+
+*/
+
+/* Local Variables: */
+   int i;
+   int j;
+   int found;
+   UnitNode **fact1;
+   double *pow1;
+   double coeff1;
+   int nfac1;
+   double con;
+
+/* Initialise */
+   *factors = NULL;
+   *powers = NULL;
+   *nfactor = 0;
+   *coeff = 1.0;
+
+/* Check inherited status. */
+   if( !astOK ) return;
+
+/* If the node at the head of the supplied tree is an OP_MULT node... */
+   if( node->opcode == OP_MULT ) {
+
+/* Find the factors of the two arguments of the OP_MULT node. */
+      FindFactors( node->arg[ 0 ], factors, powers, nfactor, coeff, status );
+      FindFactors( node->arg[ 1 ], &fact1, &pow1, &nfac1, &coeff1, status );
+
+/* Combine the two lists. Loop round the factors of the seocnd argument. */
+      for( i = 0; i < nfac1; i++ ) {
+
+/* See if there is already an equivalent factor in the returned list of
+   factors. */
+         found = 0;
+         for( j = 0; j < *nfactor; j++ ) {
+            if( !CmpTree( (*factors)[ j ], fact1[ i ], 0, status ) ){
+               found = 1;
+               break;
+            }
+         }
+
+/* If so, increment the power of the factor. */
+         if( found ) {
+            (*powers)[ j ] += pow1[ i ];
+
+/* Otherwise, add the factor to the end of the returned list. */
+         } else {
+            *factors = astGrow( *factors, *nfactor + 1, sizeof( UnitNode *) );
+            *powers = astGrow( *powers, *nfactor + 1, sizeof( double ) );
+            if( astOK ) {
+               (*factors)[ *nfactor ] = fact1[ i ];
+               (*powers)[ (*nfactor)++ ] = pow1[ i ];
+            }
+         }
+      }
+
+/* Modify the overall coefficient. */
+      *coeff *= coeff1;
+
+/* Free resources */
+      fact1 = astFree( fact1 );
+      pow1 = astFree( pow1 );
+
+/* If the node at the head of the supplied tree is an OP_POW node, */
+   } else if( node->opcode == OP_POW ) {
+
+/* Find the factors of the first argument. */
+      FindFactors( node->arg[ 0 ], factors, powers, nfactor, coeff, status );
+
+/* Multiply all the factor powers by the constant exponent of the POW
+   node. */
+      con = node->arg[ 1 ]->con;
+      for( j = 0; j < *nfactor; j++ ) {
+         (*powers)[ j ] *= con;
+      }
+
+/* Exponentiate the coefficient. */
+      if( *coeff >= 0.0 || (int) con == con ) {
+         *coeff = pow( *coeff, con );
+      } else {
+         astError( AST__BADUN, "Simplifying a units expression requires a "
+                   "negative value to be raised to a non-intergal power." , status);
+      }
+
+/* If the node at the head of the supplied tree is an OP_DIV node, */
+   } else if( node->opcode == OP_DIV ) {
+
+/* Find the factors of the two arguments of the OP_DIV node. */
+      FindFactors( node->arg[ 0 ], factors, powers, nfactor, coeff, status );
+      FindFactors( node->arg[ 1 ], &fact1, &pow1, &nfac1, &coeff1, status );
+
+/* Combine the two lists. Loop round the factors of the second argument
+   (the denominator). */
+      for( i = 0; i < nfac1; i++ ) {
+
+/* See if there is already an equivalent factor in the returned list of
+   factors. */
+         found = 0;
+         for( j = 0; j < *nfactor; j++ ) {
+            if( !CmpTree( (*factors)[ j ], fact1[ i ], 0, status ) ){
+               found = 1;
+               break;
+            }
+         }
+
+/* If so, decrement the power of the factor. */
+         if( found ) {
+            (*powers)[ j ] -= pow1[ i ];
+
+/* Otherwise, add the factor to the end of the returned list, with a
+   negated power. */
+         } else {
+            *factors = astGrow( *factors, *nfactor + 1, sizeof( UnitNode *) );
+            *powers = astGrow( *powers, *nfactor + 1, sizeof( double ) );
+            if( astOK ) {
+               (*factors)[ *nfactor ] = fact1[ i ];
+               (*powers)[ (*nfactor)++ ] = -pow1[ i ];
+            }
+         }
+      }
+
+/* Modify the overall coefficient. */
+      if( coeff1 != 0.0 ) {
+         *coeff /= coeff1;
+      } else {
+         astError( AST__BADUN, "Simplifying a units expression"
+                   "requires a division by zero." , status);
+      }
+
+/* Free resources */
+      fact1 = astFree( fact1 );
+      pow1 = astFree( pow1 );
+
+/* If the node at the head of the supplied tree is an OP_SQRT node, */
+   } else if( node->opcode == OP_SQRT ) {
+
+/* Find the factors of the argument. */
+      FindFactors( node->arg[ 0 ], factors, powers, nfactor, coeff, status );
+
+/* Multiply all the factor powers by 0.5. */
+      for( j = 0; j < *nfactor; j++ ) {
+         (*powers)[ j ] *= 0.5;
+      }
+
+/* Square root the coefficient. */
+      if( *coeff >= 0.0 ) {
+         *coeff = sqrt( *coeff );
+      } else {
+         astError( AST__BADUN, "Simplifying a units expression requires "
+                   "the square root of a negative value to be taken." , status);
+      }
+
+/* If the node at the head of the supplied tree is constant we have no
+   factors but we have a coeffcient. */
+   } else if( node->con != AST__BAD ) {
+      *coeff = node->con;
+
+/* Other nodes have no factors other than themselves, so just return a
+   pointer to the supplied node. */
+   } else {
+      *factors = astMalloc( sizeof( UnitNode *) );
+      *powers = astMalloc( sizeof( double ) );
+      if( astOK ) {
+         *nfactor = 1;
+         (*factors)[ 0 ] = node;
+         (*powers)[ 0 ] = 1.0;
+         *coeff = 1.0;
+      }
+   }
+
+/* If an error has occurred, free any returned resources. */
+   if( !astOK ) {
+      *factors = astFree( *factors );
+      *powers = astFree( *powers );
+      *nfactor = 0;
+      *coeff = 1.0;
+   }
+}
+
+static void FixConstants( UnitNode **node, int unity, int *status ) {
+/*
+*  Name:
+*     FixConstants
+
+*  Purpose:
+*     Take the reciprocal of all constants in a tree of UnitNodes.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     void FixConstants( UnitNode **node, int unity, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function replaces sub-trees which have a constant value by
+*     a single OP_LDCON node which loads the appropriate constant.
+
+*  Parameters:
+*     node
+*        The address of a pointer to the UnitNode at the head of the tree
+*        which is to be fixed. On exit the supplied tree is freed and a
+*        pointer to a new tree is palced at he given address.
+*     unity
+*        If non-zero, then all multiplicative constants are set to 1.0, and
+*        their original values are forgotten, but only if the other
+*        argument of the OP_MULT node is an OP_LDVAR, OP_POW or OP_SQRT Node.
+*     status
+*        Pointer to the inherited status variable.
+
+*/
+
+/* Local Variables: */
+   int i;
+   UnitNode *newnode;
+   int allcon;
+   Oper op;
+   double newcon;
+
+/* Check inherited status and pointer. */
+   if( !astOK || !node || !(*node) ) return;
+
+/* Initiallially, we have no replacement node */
+   newnode = NULL;
+   newcon = AST__BAD;
+
+/* There is nothing to fix if the node has no arguments. */
+   if( (*node)->narg > 0 ) {
+
+/* Note the op code for the node. */
+      op = (*node)->opcode;
+
+/* Fix up the argument nodes. Also note if all the arguments are
+   constants. */
+      allcon = 1;
+      for( i = 0; i < (*node)->narg; i++ ) {
+         FixConstants( &( (*node)->arg[ i ] ), unity, status );
+         if( (*node)->arg[ i ]->con == AST__BAD ) allcon = 0;
+      }
+
+/* If an OP_MULT nodes within a simplified tree has a constant argument,
+   it will always be argument zero.  If this is an OP_MULT node and arg[0]
+   is constant and "unity" is non-zero and arg[1] is an OP_LDVAR, OP_POW
+   or OP_SQRT node, replace the constant value by 1.0. */
+      if( unity && op == OP_MULT &&
+          (*node)->arg[ 0 ]->con != AST__BAD &&
+          ( (*node)->arg[ 1 ]->opcode == OP_LDVAR ||
+            (*node)->arg[ 1 ]->opcode == OP_SQRT ||
+            (*node)->arg[ 1 ]->opcode == OP_POW ) ) {
+         (*node)->arg[ 0 ]->con = 1.0;
+      }
+
+/* If the arguments of this node are all constants, replace the node by
+   an OP_LDCON node which loads the resulting constant value. */
+      if( allcon ) {
+         if( (*node)->narg > 0 ) {
+            newnode = NewNode( NULL, OP_LDCON, status );
+            if( astOK ) {
+               if( op == OP_LOG ) {
+                  if( (*node)->arg[ 0 ]->con > 0.0 ) {
+                     newcon = log10( (*node)->arg[ 0 ]->con );
+                  } else {
+                     astError( AST__BADUN, "Illegal negative or zero constant "
+                               "value '%g' encountered.", status,
+                               (*node)->arg[ 0 ]->con );
+                  }
+               } else if( op == OP_LN ){
+                  if( (*node)->arg[ 0 ]->con > 0.0 ) {
+                     newcon = log( (*node)->arg[ 0 ]->con );
+                  } else {
+                     astError( AST__BADUN, "Illegal negative or zero constant value "
+                               "'%g' encountered.", status, (*node)->arg[ 0 ]->con );
+                  }
+               } else if( op == OP_EXP ){
+                  newcon = exp( (*node)->arg[ 0 ]->con );
+
+               } else if( op == OP_SQRT ){
+                  if( (*node)->arg[ 0 ]->con >= 0.0 ) {
+                     newcon = sqrt( (*node)->arg[ 0 ]->con );
+                  } else {
+                     astError( AST__BADUN, "Illegal negative constant value "
+                               "'%g' encountered.", status, (*node)->arg[ 0 ]->con );
+                  }
+
+               } else if( op == OP_POW ){
+                  if( (*node)->arg[ 0 ]->con >= 0.0 ||
+                      (int) (*node)->arg[ 1 ]->con == (*node)->arg[ 1 ]->con ) {
+                     newcon = pow( (*node)->arg[ 0 ]->con,
+                                   (*node)->arg[ 1 ]->con );
+                  } else {
+                     astError( AST__BADUN, "Illegal negative constant value "
+                               "'%g' encountered.", status, (*node)->arg[ 0 ]->con );
+                  }
+
+               } else if( op == OP_DIV ){
+                  if( (*node)->arg[ 1 ]->con != 0.0 ) {
+                     newcon = (*node)->arg[ 0 ]->con / (*node)->arg[ 1 ]->con;
+                  } else {
+                     astError( AST__BADUN, "Illegal zero constant value encountered." , status);
+                  }
+
+               } else if( op == OP_MULT ){
+                  newcon = (*node)->arg[ 0 ]->con * (*node)->arg[ 1 ]->con;
+
+               }
+
+
+               if( astOK ) newnode->con = newcon;
+            }
+         }
+      }
+   }
+
+/* If an error has occurred, free any new node. */
+   if( !astOK ) newnode = FreeTree( newnode, status );
+
+/* If we have a replacement node, free the supplied tree and return a
+   pointer to the new tree. */
+   if( newnode ) {
+      FreeTree( *node, status );
+      *node = newnode;
+   }
+
+}
+
+static UnitNode *FixUnits( UnitNode *node, UnitNode *test, int *status ) {
+/*
+*  Name:
+*     FixUnits
+
+*  Purpose:
+*     Assign a constant value to all units except for one.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     UnitNode *FixUnits( UnitNode *node, UnitNode *test, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function returns a copy of the supplied tree of UnitNodes. All
+*     OP_LDVAR nodes within the copy which refer to units which differ
+*     from those referred to by the supplied test node are replaced by
+*     OP_LDCON nodes which load the constant value 1.0.
+
+*  Parameters:
+*     node
+*        A pointer to the UnitNode at the head of the tree to be used.
+*     test
+*        A pointer to an OP_LDVAR node which defines the units which are
+*        *not* to be replaced by a constant value of 1.0.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     A pointer to a UnitNode which is at the head of a tree of UnitNodes
+*     which forms the required copy of th einput tree.
+
+* Notes:
+*     - A NULL pointer is returned if an error has already occurred, or
+*     if this function fails for any reason.
+
+*/
+
+/* Local Variables: */
+   int i;
+   UnitNode *result;
+
+/* Initialise. */
+   result = NULL;
+
+/* Check inherited status. */
+   if( !astOK ) return result;
+
+/* Create a complete copy of the supplied tree. */
+   result = CopyTree( node, status );
+
+/* Is the node at the head of the supplied tree an OP_LDVAR node? */
+   if( node->opcode == OP_LDVAR ) {
+
+/* Does it refer to a unit which differs from that of the test node? If so
+   annul the copy created above and return a new OP_LDCON node which loads
+   the constant value 1.0. */
+      if( strcmp( test->name, node->name ) ) {
+         FreeTree( result, status );
+         result = NewNode( NULL, OP_LDCON, status );
+         if( astOK ) result->con = 1.0;
+      }
+
+/* If the supplied node is not an OP_LDVAR node, check each argument of
+   the head node. */
+   } else {
+      for( i = 0; i < node->narg; i++ ) {
+
+/* Free the resources used to hold this argument in the tree copy created
+   above. */
+         FreeTree( result->arg[ i ], status );
+
+/* Create a new argument tree by calling this function recursively to
+   fix units in the argument sub-trees. */
+         result->arg[ i ] = FixUnits( node->arg[ i ], test, status );
+      }
+   }
+
+/* If an error has occurred, free any new tree. */
+   if( !astOK ) result = FreeTree( result, status );
+
+/* Return the answer. */
+   return result;
+}
+
+static UnitNode *FreeTree( UnitNode *node, int *status ) {
+/*
+*  Name:
+*     FreeTree
+
+*  Purpose:
+*     Free resources used by a tree of UnitNodes.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     UnitNode *FreeTree( UnitNode *node, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function frees the memory used to store a tree of UnitNodes.
+
+*  Parameters:
+*     node
+*        A pointer to the UnitNode at the head of the tree which is to be
+*        freed.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     A NULL pointer is returned.
+
+*  Notes:
+*     - This function attempts to execute even if it is invoked with
+*     the global error status set.
+*/
+
+/* Local Variables: */
+   int i;
+
+/* Check a node was supplied. */
+   if( node ) {
+
+/* Recursively free any argument nodes. */
+      if( node->arg ) {
+         for( i = 0; i < node->narg; i++ ) {
+            (node->arg)[ i ] = FreeTree( (node->arg)[ i ], status );
+         }
+         node->arg = astFree( node->arg );
+      }
+
+/* Nullify other pointers for safety. */
+      node->unit = NULL;
+      node->mult = NULL;
+
+/* Free the copy of the symbol string (if any). */
+      node->name = astFree( (char *) node->name );
+
+/* Free the memory holding the node. */
+      node = astFree( node );
+   }
+
+/* Return a null pointer. */
+   return NULL;
+}
+
+static KnownUnit *GetKnownUnits( int lock, int *status ) {
+/*
+*  Name:
+*     GetKnownUnits
+
+*  Purpose:
+*     Get a pointer to the head of a linked list of known unit definitions.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     KnownUnit *GetKnownUnits( int lock, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function returns a pointer to the head of a linked list of known
+*     unit definitions. The unit definitions are created as static module
+*     variables if they have not previously been created.
+
+*  Parameters:
+*     lock
+*        If non-zero, then lock a mutex prior to accessing the list of
+*        known units.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     A pointer to the first known unit definition.
+
+*  Notes:
+*     - A NULL pointer is returned if it is invoked with the global error
+*     status set, or if an error occurs.
+*/
+
+/* Local Variables: */
+   int iq;
+   KnownUnit *result;
+
+/* Initialise. */
+   result = NULL;
+
+/* Check inherited status. */
+   if( !astOK ) return result;
+
+/* Ensure the known units list is only initialised once. */
+   if( lock ) {
+      LOCK_MUTEX1
+   }
+
+/* If the linked list of KnownUnit structures describing the known units
+   has not yet been created, create it now. A pointer to the head of the
+   linked list is put into the static variable "known_units". */
+   if( !known_units ) {
+
+/* At the same time we store pointers to the units describing the basic
+   quantities used in dimensional analysis. Initialise th index of the
+   next such unit. */
+      iq = 0;
+
+/* Create definitions for the known units. First do all IAU basic units.
+   We include "g" instead of "kg" because otherwise we would have to
+   refer to a gramme as a milli-kilogramme. */
+      MakeKnownUnit( "g", "gram", NULL, status );
+      quant_units[ iq++ ] = known_units;
+      MakeKnownUnit( "m", "metre", NULL, status );
+      quant_units[ iq++ ] = known_units;
+      MakeKnownUnit( "s", "second", NULL, status );
+      quant_units[ iq++ ] = known_units;
+      MakeKnownUnit( "rad", "radian", NULL, status );
+      quant_units[ iq++ ] = known_units;
+      MakeKnownUnit( "K", "Kelvin", NULL, status );
+      quant_units[ iq++ ] = known_units;
+      MakeKnownUnit( "A", "Ampere", NULL, status );
+      MakeKnownUnit( "mol", "mole", NULL, status );
+      MakeKnownUnit( "cd", "candela", NULL, status );
+
+/* Now do all IAU derived units. Unit definitions may only refer to units
+   which have already been defined. */
+      MakeKnownUnit( "sr", "steradian", "rad rad", status );
+      MakeKnownUnit( "Hz", "Hertz", "1/s", status );
+      MakeKnownUnit( "N", "Newton", "kg m/s**2", status );
+      MakeKnownUnit( "J", "Joule", "N m", status );
+      MakeKnownUnit( "W", "Watt", "J/s", status );
+      MakeKnownUnit( "C", "Coulomb", "A s", status );
+      MakeKnownUnit( "V", "Volt", "J/C", status );
+      MakeKnownUnit( "Pa", "Pascal", "N/m**2", status );
+      MakeKnownUnit( "Ohm", "Ohm", "V/A", status );
+      MakeKnownUnit( "S", "Siemens", "A/V", status );
+      MakeKnownUnit( "F", "Farad", "C/V", status );
+      MakeKnownUnit( "Wb", "Weber", "V s", status );
+      MakeKnownUnit( "T", "Tesla", "Wb/m**2", status );
+      MakeKnownUnit( "H", "Henry", "Wb/A", status );
+      MakeKnownUnit( "lm", "lumen", "cd sr", status );
+      MakeKnownUnit( "lx", "lux", "lm/m**2", status );
+
+/* Now do additional derived and basic units listed in the FITS-WCS paper. */
+      MakeKnownUnit( "deg", "degree", "pi/180 rad", status );
+      MakeKnownUnit( "arcmin", "arc-minute", "1/60 deg", status );
+      MakeKnownUnit( "arcsec", "arc-second", "1/3600 deg", status );
+      MakeKnownUnit( "mas", "milli-arcsecond", "1/3600000 deg", status );
+      MakeKnownUnit( "min", "minute", "60 s", status );
+      MakeKnownUnit( "h", "hour", "3600 s", status );
+      MakeKnownUnit( "d", "day", "86400 s", status );
+      MakeKnownUnit( "yr", "year", "31557600 s", status );
+      MakeKnownUnit( "a", "year", "31557600 s", status );
+      MakeKnownUnit( "eV", "electron-Volt", "1.60217733E-19 J", status );
+      MakeKnownUnit( "erg", "erg", "1.0E-7 J", status );
+      MakeKnownUnit( "Ry", "Rydberg", "13.605692 eV", status );
+      MakeKnownUnit( "solMass", "solar mass", "1.9891E30 kg", status );
+      MakeKnownUnit( "u", "unified atomic mass unit", "1.6605387E-27 kg", status );
+      MakeKnownUnit( "solLum", "solar luminosity", "3.8268E26 W", status );
+      MakeKnownUnit( "Angstrom", "Angstrom", "1.0E-10 m", status );
+      MakeKnownUnit( "micron", "micron", "1.0E-6 m", status );
+      MakeKnownUnit( "solRad", "solar radius", "6.9599E8 m", status );
+      MakeKnownUnit( "AU", "astronomical unit", "1.49598E11 m", status );
+      MakeKnownUnit( "lyr", "light year", "9.460730E15 m", status );
+      MakeKnownUnit( "pc", "parsec", "3.0867E16 m", status );
+      MakeKnownUnit( "count", "count", NULL, status );
+      quant_units[ iq++ ] = known_units;
+      MakeKnownUnit( "adu", "analogue-to-digital unit", NULL, status );
+      quant_units[ iq++ ] = known_units;
+      MakeKnownUnit( "photon", "photon", NULL, status );
+      quant_units[ iq++ ] = known_units;
+      MakeKnownUnit( "Jy", "Jansky", "1.0E-26 W /m**2 /Hz", status );
+      MakeKnownUnit( "mag", "magnitude", NULL, status );
+      quant_units[ iq++ ] = known_units;
+      MakeKnownUnit( "G", "Gauss", "1.0E-4 T", status );
+      MakeKnownUnit( "pixel", "pixel", NULL, status );
+      quant_units[ iq++ ] = known_units;
+      MakeKnownUnit( "barn", "barn", "1.0E-28 m**2", status );
+      MakeKnownUnit( "D", "Debye", "(1.0E-29/3) C.m", status );
+
+      if( iq != NQUANT && astOK ) {
+         astError( AST__INTER, "unit(GetKnownUnits): %d basic quantities "
+                   "noted but this should be %d (internal AST programming "
+                   "error).", status, iq, NQUANT );
+      }
+
+/* Unit aliases... */
+      MakeUnitAlias( "Angstrom", "Ang", status );
+      MakeUnitAlias( "count", "ct", status );
+      MakeUnitAlias( "photon", "ph", status );
+      MakeUnitAlias( "Jy", "Jan", status );
+      MakeUnitAlias( "pixel", "pix", status );
+      MakeUnitAlias( "s", "sec", status );
+      MakeUnitAlias( "m", "meter", status );
+   }
+
+/* If succesful, return the pointer to the head of the list. */
+   if( astOK ) result = known_units;
+
+/* Allow the next thread to proceed. */
+   if( lock ) {
+      UNLOCK_MUTEX1
+   }
+
+/* Return the result. */
+   return result;
+}
+
+static Multiplier *GetMultipliers( int *status ) {
+/*
+*  Name:
+*     GetMultiplier
+
+*  Purpose:
+*     Get a pointer to the head of a linked list of multiplier definitions.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     Multiplier *Multipliers( void )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function returns a pointer to the head of a linked list of known
+*     multiplier definitions. The multiplier definitions are created as
+*     static module variables if they have not previously been created.
+
+*  Returned Value:
+*     A pointer to the first known multiplier definition.
+
+*  Notes:
+*     - A NULL pointer is returned if it is invoked with the global error
+*     status set, or if an error occurs.
+*/
+
+/* Local Variables: */
+   Multiplier *result;
+   Multiplier *mult;
+
+/* Initialise. */
+   result = NULL;
+
+/* Check inherited status. */
+   if( !astOK ) return result;
+
+/* Ensure the list is only initialised by one thread. */
+   LOCK_MUTEX2
+
+/* If the linked list of Multiplier structures describing the known
+   multipliers has not yet been created, create it now. A pointer to the
+   head of the linked list is put into the static variable "multipliers". */
+   if( !multipliers ) {
+
+/* Define a macro to create a multiplier struncture and add it to the
+   linked list of multiplier structures. */
+#define MAKEMULT(s,sl,sc,lab,ll) \
+      mult = astMalloc( sizeof( Multiplier ) ); \
+      if( astOK ) { \
+         mult->sym = s; \
+         mult->symlen = sl; \
+         mult->lablen = ll; \
+         mult->scale = sc; \
+         mult->label = lab; \
+         mult->next = multipliers; \
+         multipliers = mult; \
+      }
+
+/* Use the above macro to create all the standard multipliers listed in the
+   FITS WCS paper I. */
+      MAKEMULT("d",1,1.0E-1,"deci",4)
+      MAKEMULT("c",1,1.0E-2,"centi",5)
+      MAKEMULT("m",1,1.0E-3,"milli",5)
+      MAKEMULT("u",1,1.0E-6,"micro",5)
+      MAKEMULT("n",1,1.0E-9,"nano",4)
+      MAKEMULT("p",1,1.0E-12,"pico",4)
+      MAKEMULT("f",1,1.0E-15,"femto",5)
+      MAKEMULT("a",1,1.0E-18,"atto",4)
+      MAKEMULT("z",1,1.0E-21,"zepto",5)
+      MAKEMULT("y",1,1.0E-24,"yocto",5)
+      MAKEMULT("da",2,1.0E1,"deca",4)
+      MAKEMULT("h",1,1.0E2,"hecto",5)
+      MAKEMULT("k",1,1.0E3,"kilo",4)
+      MAKEMULT("M",1,1.0E6,"mega",4)
+      MAKEMULT("G",1,1.0E9,"giga",4)
+      MAKEMULT("T",1,1.0E12,"tera",4)
+      MAKEMULT("P",1,1.0E15,"peta",4)
+      MAKEMULT("E",1,1.0E18,"exa",3)
+      MAKEMULT("Z",1,1.0E21,"zetta",5)
+      MAKEMULT("Y",1,1.0E24,"yotta",5)
+
+/* Undefine the macro. */
+#undef MAKEMULT
+
+   }
+
+/* If succesful, return the pointer to the head of the list. */
+   if( astOK ) result = multipliers;
+
+/* Allow the next thread to proceed. */
+   UNLOCK_MUTEX2
+
+/* Return the result. */
+   return result;
+}
+
+static void InvertConstants( UnitNode **node, int *status ) {
+/*
+*  Name:
+*     InvertConstants
+
+*  Purpose:
+*     Take the reciprocal of all constants in a tree of UnitNodes.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     void InvertConstants( UnitNode **node, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function replaces constant unit coefficients by their reciprocal.
+*     This is because a string such as "0.01 m" will be interpreted as
+*     meaning "multiply a value in metres by 0.01 to get the value in the
+*     required units", whereas what is actually meant is "use units of
+*     0.01 of a metre" which requires us to divide the value in metres by
+*     0.01, not multiply it.
+
+*  Parameters:
+*     node
+*        The address of a pointer to the UnitNode at the head of the tree.
+*        On exit the supplied tree is freed and a pointer to a new tree is
+*        placed at the given address.
+*     status
+*        Pointer to the inherited status variable.
+
+*/
+
+/* Local Variables: */
+   int i;
+   UnitNode *newnode;
+   int allcon;
+   Oper op;
+
+/* Check inherited status and pointer. */
+   if( !astOK || !node || !(*node) ) return;
+
+/* Initiallially, we have no replacement node */
+   newnode = NULL;
+
+/* There is nothing to fix if the node has no arguments. */
+   if( (*node)->narg > 0 ) {
+
+/* Note the op code for the node. */
+      op = (*node)->opcode;
+
+/* Fix up the argument nodes. Also note if all the arguments are
+   constants. */
+      allcon = 1;
+      for( i = 0; i < (*node)->narg; i++ ) {
+         InvertConstants( &( (*node)->arg[ i ] ), status );
+         if( (*node)->arg[ i ]->con == AST__BAD ) allcon = 0;
+      }
+
+/* If all nodes are constant, there are no co-efficients to invert. */
+      if( !allcon ) {
+
+/* Iif this is a multiplication node, see if either of its arguments
+   is a constant. If so, invert the constant. This is because a string like
+   "0.01 m" means "each unit is 0.01 of a metre". Therefore, to transform
+   a value in metres into required units means multiplying the metres
+   value by 100.0 (i.e the reciprocal of 0.01), not 0.01. */
+         if( op == OP_MULT ) {
+            for( i = 0; i < 2; i++ ) {
+               if( (*node)->arg[ i ]->con != AST__BAD ) {
+                  if( (*node)->arg[ i ]->con != 0.0 ) {
+
+                     (*node)->arg[ i ]->con = 1.0/(*node)->arg[ i ]->con;
+                  } else {
+                     astError( AST__BADUN, "Illegal zero constant encountered." , status);
+                  }
+               }
+            }
+
+/* Likewise, check for division nodes in which the denominator is
+   constant. */
+         } else if( op == OP_DIV ) {
+            if( (*node)->arg[ 1 ]->con != AST__BAD ) {
+               if( (*node)->arg[ 1 ]->con != 0.0 ) {
+                  (*node)->arg[ 1 ]->con = 1.0/(*node)->arg[ 1 ]->con;
+               } else {
+                  astError( AST__BADUN, "Illegal zero constant encountered." , status);
+               }
+            }
+
+/* If this is a "pow" node check that the second argument is constant
+   (required by FITS WCS paper I). */
+         } else if( op == OP_POW ) {
+            if( (*node)->arg[ 1 ]->con == AST__BAD ) {
+               astError( AST__BADUN, "Illegal variable exponent." , status);
+            }
+         }
+      }
+   }
+
+/* If an error has occurred, free any new node. */
+   if( !astOK ) newnode = FreeTree( newnode, status );
+
+/* If we have a replacement node, free the supplied tree and return a
+   pointer to the new tree. */
+   if( newnode ) {
+      FreeTree( *node, status );
+      *node = newnode;
+   }
+}
+
+static UnitNode *InvertTree( UnitNode *fwdnode, UnitNode *src, int *status ) {
+/*
+*  Name:
+*     InvertTree
+
+*  Purpose:
+*     Invert a tree of UnitNodes.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     UnitNode *InvertTree( UnitNode *fwdnode, UnitNode *src )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function inverts a tree of UnitNodes. The supplied tree should
+*     have exactly one OP_LDVAR node. This will be the quantity represented
+*     by the node at the head of the returned tree.
+
+*  Parameters:
+*     fwdnode
+*        A pointer to the UnitNode at the head of the tree which is to be
+*        inverted.
+*     src
+*        A pointer to a UnitNode which is to be used as the root of the
+*        inverted tree. That is, the output from this node should form
+*        the (one and only) varying input to the inverted tree. If the
+*        supplied tree is succesfulyl inverted, the tree of which "src"
+*        is the head will be contained within the returned inverted tree.
+*        Therefore "src" only needs to be freed explicitly if this
+*        function fails to invert the supplied tree for any reason. If
+*        this function succeeds, then "src" will be freed as part of
+*        freeing the returned inverted tree.
+
+*  Returned Value:
+*     A pointer to a UnitNode which forms the head of the inverted tree.
+
+*  Algorithm:
+*     The algorithm works through the supplied forward tree, from the head
+*     to the roots. First, the supplied node at the head of the forward
+*     tree is inverted. To be invertable, the supplied head node must have
+*     exactly one varying argument (any other arguments must be fixed,
+*     i.e. not vary). This varying argument becomes the output of the
+*     inverted node. The other (fixed) arguments to the forward node are
+*     also used as arguments to the inverted node. The supplied "src" node
+*     is used as the single varying input to the inverted node. Having
+*     inverted the supplied forward head node, this function is called
+*     recursively to invert the lower parts of the forward tree (i.e. the
+*     part of the forward tree which provided the varying input to node
+*     which has just been inverted).
+
+*  Notes:
+*     - It is assumed that he supplied forward tree has been simplified
+*     using SimplifyTree. This means that the tree contains no nodes with
+*     the following op codes: OP_LOG, OP_SQRT. OP_DIV (SimplifyTree
+*     converts these nodes into OP_LN, OP_POW and OP_MULT nodes).
+*     - A value of NULL will be returned if this function is invoked with
+*     the global error status set, or if it should fail for any reason.
+
+*/
+
+/* Local Variables: */
+   UnitNode *newnode;
+   UnitNode *nextnode;
+   UnitNode *result;
+   UnitNode *node1;
+   Oper fop;
+   int varg;
+
+/* Initialise */
+   result = NULL;
+
+/* Check inherited status. */
+   if( !astOK ) return result;
+
+/* Initiallially, we have no replacement node */
+   newnode = NULL;
+   nextnode = NULL;
+
+/* Save the op code at the head of the forward tree. */
+   fop = fwdnode->opcode;
+
+/* If the head of the forward tree is a OP_EXP node. Inverse of
+   "exp(x)" is "ln(x)". */
+   if( fop == OP_EXP ) {
+      newnode = NewNode( NULL, OP_LN, status );
+      if( astOK ) {
+         newnode->arg[ 0 ] = src;
+         nextnode = fwdnode->arg[ 0 ];
+      }
+
+/* If the head of the forward tree is a OP_LN node. Inverse of
+   "ln(x)" is "exp(x)". */
+   } else if( fop == OP_LN ) {
+      newnode = NewNode( NULL, OP_EXP, status );
+      if( astOK ) {
+         newnode->arg[ 0 ] = src;
+         nextnode = fwdnode->arg[ 0 ];
+      }
+
+/* If the head of the forward tree is a OP_POW node. Inverse of
+   "x**k" is "x**(1/k)" */
+   } else if( fop == OP_POW ) {
+      newnode = NewNode( NULL, OP_POW, status );
+      node1 = NewNode( NULL, OP_LDCON, status );
+      if( astOK ) {
+         node1->con = 1.0/fwdnode->arg[ 1 ]->con;
+         newnode->arg[ 0 ] = src;
+         newnode->arg[ 1 ] = node1;
+         nextnode = fwdnode->arg[ 0 ];
+      }
+
+/* If the head of the forward tree is a OP_MULT node... */
+   } else if( fop == OP_MULT ) {
+
+/* The node is only invertable if it has one constant node and one
+   non-constant node. Get the index of the varying argument. */
+      if( fwdnode->arg[ 0 ]->con != AST__BAD &&
+          fwdnode->arg[ 1 ]->con == AST__BAD ) {
+         varg = 1;
+      } else if( fwdnode->arg[ 0 ]->con == AST__BAD &&
+                 fwdnode->arg[ 1 ]->con != AST__BAD ) {
+         varg = 0;
+      } else {
+         varg = -1;
+      }
+      if( varg != -1 ) {
+
+/* The inverse of "k*x" is "(1/k)*x" (we use MULT nodes instead of DIV
+   nodes to maintain the standardisation implemented by SimplifyTree). */
+         newnode = NewNode( NULL, OP_MULT, status );
+         node1 = NewNode( NULL, OP_LDCON, status );
+         if( astOK ) {
+            node1->con = 1.0/fwdnode->arg[ 1 - varg ]->con;
+            newnode->arg[ 0 ] = node1;
+            newnode->arg[ 1 ] = src;
+            nextnode = fwdnode->arg[ varg ];
+         }
+      }
+
+/* If the head of the forward tree is a OP_LDVAR node, there is nothing
+   left to invert. SO return a pointer to the suppleid source node. */
+   } else if( fop == OP_LDVAR ) {
+      result = src;
+      nextnode = NULL;
+
+/* If the head of the forward tree is any other node (e.g. a OP_LDCON node),
+   the tree cannot be inverted. */
+   } else {
+      nextnode = NULL;
+   }
+
+/* If we managed to invert the node at the head of the supplied tree,
+   continue to invert its varying argument node (if any). */
+   if( nextnode && newnode ) result = InvertTree( nextnode, newnode, status );
+
+/* If the tree could not be inverted, free the newnode. */
+   if( !result ) newnode = FreeTree( newnode, status );
+
+/* If an error has occurred, free any new node. */
+   if( !astOK ) result = FreeTree( result, status );
+
+/* Return the result. */
+   return result;
+
+}
+
+static void LocateUnits( UnitNode *node, UnitNode ***units, int *nunits, int *status ){
+/*
+*  Name:
+*     LocateUnits
+
+*  Purpose:
+*     Locate the units used by a supplied tree of UnitNodes.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     void LocateUnits( UnitNode *node, UnitNode ***units, int *nunits, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function locates the units used by a supplied tree of
+*     UnitNodes.
+
+*  Parameters:
+*     node
+*        A pointer to the UnitNode at the head of the tree to be searched.
+*     units
+*        The address at which is stored a pointer to an array of "*nunits"
+*        elements. Each element of the array holds a pointer to a UnitNode.
+*        The array is extended on exit to hold pointers to the OP_LDVAR nodes
+*        within the supplied tree (i.e. nodes which represent named units,
+*        either known or unknown). A node is only included in the returned
+*        array if no other node for the same unit is already included in the
+*        array. A NULL pointer should be supplied on the first invocation of
+*        this function.
+*     nunits
+*        The address of an integer which holds the number of elements in
+*        the array given by "*units". Updated on exit to included any
+*        elements added to the array. Zero should be supplied on the first
+*        invocation of this function.
+*     status
+*        Pointer to the inherited status variable.
+
+*/
+
+/* Local Variables: */
+   int i;
+   int found;
+
+/* Check the global error status. */
+   if( !astOK ) return;
+
+/* Is the node at the head of the supplied tree an OP_LDVAR node? */
+   if( node->opcode == OP_LDVAR ) {
+
+/* If an array was supplied, see if it already contains a pointer to a node
+   which refers to the same units. */
+      found = 0;
+      if( *units ) {
+         for( i = 0; i < *nunits; i++ ) {
+            if( !strcmp( (*units)[ i ]->name, node->name ) ) {
+               found = 1;
+               break;
+            }
+         }
+      }
+
+/* If not, ensure the array is big enough and add a pointer to the
+   supplied node to the array. */
+      if( !found ) {
+         *units = astGrow( *units, *nunits + 1, sizeof( UnitNode * ) );
+         if( astOK ) (*units)[ (*nunits)++ ] = node;
+      }
+
+/* If the supplied node is not an OP_LDVAR node, call this function
+   recursively to search the argument sub-trees. */
+   } else {
+      for( i = 0; i < node->narg; i++ ) {
+         LocateUnits( node->arg[ i ], units, nunits, status );
+      }
+   }
+}
+
+static const char *MakeExp( UnitNode *tree, int mathmap, int top, int *status ) {
+/*
+*  Name:
+*     MakeExp
+
+*  Purpose:
+*     Make an algebraic expression from a supplied tree of UnitNodes.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     const char *MakeExp( UnitNode *tree, int mathmap, int top, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function produces a string holding an algebraic expression
+*     corresponding to a supplied tree of UnitNodes.
+
+*  Parameters:
+*     tree
+*        A pointer to the UnitNode at the head of the tree to be converted
+*        into an algebraic expression.
+*     mathmap
+*        If zero, format as an axis label expression. If 1, format as a
+*        MathMap expression. If 2, format as a FITS unit string.
+*     top
+*        Should be non-zero for a top-level entry to this function, and
+*        zero for a recursive entry.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     A pointer to the cleaned expression, which should be freed using
+*     astFree when no longer needed.
+
+*  Notes:
+*     - This function returns NULL if it is invoked with the global error
+*     status set, or if it should fail for any reason.
+*/
+
+/* Local Variables: */
+   UnitNode *newtree;
+   UnitNode *sunit;
+   char *a;
+   char *result;
+   char buff[200];
+   const char *arg0;
+   const char *arg1;
+   const char *mtxt;
+   int larg0;
+   int larg1;
+   int lbuff;
+   int mlen;
+   int par;
+   int tlen;
+
+/* Check inherited status. */
+   result = NULL;
+   if( !astOK ) return result;
+
+/* Modify the tree to make the resulting transformation functions more
+   natural to human readers. */
+   newtree = CopyTree( tree, status );
+   ComplicateTree( &newtree, status );
+
+/* If we are producing an axis label... */
+   if( !mathmap ) {
+
+/* Fix all multiplicative constants to 1.0 if they multiply an OP_LDVAR
+   OP_SQRT or OP_POW node. This is on the assumption that the returned label
+   should not include any simple unit scaling (e.g. if the output label would
+   be "2.345*wavelength", we prefer simply to use "wavelength" since a scaled
+   wavelength is still a wavelength - i.e. simple scaling does not change
+   the dimensions of a quantity). */
+      FixConstants( &newtree, 1, status );
+
+/* Simplify the tree again to get rid of the 1.0 terms which may have
+   been introduced by the previous line (but do not re-introduce any
+   standardisations - removing them was the reason for calling ComplicateTree).
+   If this simplication introduces any changes, try fixing multiplicative
+   constants again, and so on, until no more changes occur. */
+      while( SimplifyTree( &newtree, 0, status ) ) {
+         FixConstants( &newtree, 1, status );
+      }
+
+   }
+
+/* Produce a string describing the action performed by the UnitNode at
+   the head of the supplied tree, and then invoke this function recursively
+   to format any arguments of the head node. */
+
+/* Constant valued nodes... just format the constant in a local buffer and
+   then copy the buffer. */
+   if( newtree->con != AST__BAD ) {
+      lbuff = sprintf( buff, "%.*g", AST__DBL_DIG, newtree->con );
+      result = astStore( NULL, buff, lbuff + 1 );
+
+/* "Load Variable Value" nodes - return the variable name. If this is a
+   recursive call to this function, and we are producing a label, append a
+   single space before and after the name. */
+   } else if( newtree->opcode ==  OP_LDVAR ) {
+      tlen = strlen( newtree->name );
+
+      if( !mathmap && !top ){
+         result = astMalloc( tlen + 3 );
+         if( result ) {
+            result[ 0 ] = ' ';
+            memcpy( result + 1, newtree->name, tlen );
+            memcpy( result + tlen + 1, " ", 2 );
+         }
+
+      } else if( mathmap == 2 ) {
+
+         if( newtree->mult ) {
+           mlen = newtree->mult->symlen;
+           mtxt = newtree->mult->sym;
+         } else {
+           mlen = 0;
+           mtxt = NULL;
+         }
+
+         result = astMalloc( tlen + 1 + mlen );
+         if( result ) {
+            if( mtxt ) memcpy( result, mtxt, mlen );
+            memcpy( result + mlen, newtree->name, tlen + 1 );
+         }
+
+      } else {
+         result = astStore( NULL, newtree->name, tlen + 1 );
+      }
+
+/* Single argument functions... place the argument in parentheses after
+   the function name. */
+   } else if( newtree->opcode ==  OP_LOG ) {
+      arg0 = MakeExp( newtree->arg[ 0 ], mathmap, 0, status );
+      larg0 = strlen( arg0 );
+      if( mathmap == 1 ) {
+         result = astMalloc( larg0 + 8 );
+         if( result ) memcpy( result, "log10(", 7 );
+         a = result + 6;
+      } else {
+         result = astMalloc( larg0 + 6 );
+         if( result ) memcpy( result, "log(", 5 );
+         a = result + 4;
+      }
+      if( result ){
+         memcpy( a, arg0, larg0 + 1 );
+         memcpy( a + larg0, ")", 2 );
+      }
+      arg0 = astFree( (void *) arg0 );
+
+   } else if( newtree->opcode ==  OP_LN ) {
+      arg0 = MakeExp( newtree->arg[ 0 ], mathmap, 0, status );
+      larg0 = strlen( arg0 );
+      if( mathmap == 1 ) {
+         result = astMalloc( larg0 + 6 );
+         if( result ) memcpy( result, "log(", 5 );
+         a = result + 4;
+      } else {
+         result = astMalloc( larg0 + 5 );
+         if( result ) memcpy( result, "ln(", 4 );
+         a = result + 3;
+      }
+      if( astOK ){
+         memcpy( a, arg0, larg0 );
+         memcpy( a + larg0, ")", 2 );
+      }
+      arg0 = astFree( (void *) arg0 );
+
+   } else if( newtree->opcode ==  OP_EXP ) {
+      arg0 = MakeExp( newtree->arg[ 0 ], mathmap, 0, status );
+      larg0 = strlen( arg0 );
+      result = astMalloc( larg0 + 6 );
+      if( result ){
+         memcpy( result, "exp(", 5 );
+         memcpy( result + 4, arg0, larg0 );
+         memcpy( result + 4 + larg0, ")", 2 );
+      }
+      arg0 = astFree( (void *) arg0 );
+
+   } else if( newtree->opcode ==  OP_SQRT ) {
+      arg0 = MakeExp( newtree->arg[ 0 ], mathmap, 0, status );
+      larg0 = strlen( arg0 );
+      result = astMalloc( larg0 + 7 );
+      if( result ){
+         memcpy( result, "sqrt(", 6 );
+         memcpy( result + 5, arg0, larg0 );
+         memcpy( result + 5 + larg0, ")", 2 );
+      }
+      arg0 = astFree( (void *) arg0 );
+
+/* POW... the exponent (arg[1]) is always a constant and so does not need
+   to be placed in parentheses. The first argument only needs to be
+   placed in parentheses if it is a two arg node (except we also put it
+   in parentheses if it is an OP_LDVAR node and "mathmap" is zero - this is
+   because such OP_LDVAR nodes will correspond to axis labels which will
+   have spaces before and after them which would look odd if not encloses
+   in parentheses). */
+   } else if( newtree->opcode ==  OP_POW ) {
+
+      arg0 = MakeExp( newtree->arg[ 0 ], mathmap, 0, status );
+      larg0 = strlen( arg0 );
+
+      arg1 = MakeExp( newtree->arg[ 1 ], mathmap, 0, status );
+      larg1 = strlen( arg1 );
+
+      if( newtree->arg[ 0 ]->narg == 2 ||
+          (newtree->arg[ 0 ]->opcode == OP_LDVAR && !mathmap) ) {
+         par = 1;
+         result = astMalloc( larg0 + larg1 + 7 );
+         if( result ) memcpy( result, "(", 2 );
+         a = result + 1;
+      } else {
+         par = 0;
+         result = astMalloc( larg0 + larg1 + 5 );
+         a = result;
+      }
+
+      if( result ) {
+         memcpy( a, arg0, larg0 );
+         a += larg0;
+         if( par ) *(a++) = ')';
+         memcpy( a, "**", 3 );
+         a += 2;
+         memcpy( a, arg1, larg1 );
+         a += larg1;
+         *a = 0;
+      }
+
+      arg0 = astFree( (void *) arg0 );
+      arg1 = astFree( (void *) arg1 );
+
+/* DIV... the first argument (numerator) never needs to be in parentheses.
+   The second argument (denominator) only needs to be placed in parentheses
+   if it is a MULT node. */
+   } else if( newtree->opcode ==  OP_DIV ) {
+
+      if( mathmap == 2 && ( sunit = ModifyPrefix( newtree, status ) ) ) {
+         result = (char *) MakeExp( sunit, mathmap, 0, status );
+         sunit = FreeTree( sunit, status );
+
+      } else {
+         arg0 = MakeExp( newtree->arg[ 0 ], mathmap, 0, status );
+         larg0 = strlen( arg0 );
+
+         arg1 = MakeExp( newtree->arg[ 1 ], mathmap, 0, status );
+         larg1 = strlen( arg1 );
+
+         if( newtree->arg[ 1 ]->opcode == OP_MULT &&
+             strchr( arg1, '*' ) ) {
+            par = 1;
+            result = astMalloc( larg0 + larg1 + 4 );
+         } else {
+            par = 0;
+            result = astMalloc( larg0 + larg1 + 2 );
+         }
+
+         if( result ) {
+            memcpy( result, arg0, larg0 );
+            a = result + larg0;
+            *(a++) = '/';
+            if( par ) *(a++) = '(';
+            memcpy( a, arg1, larg1 );
+            a += larg1;
+            if( par ) *(a++) = ')';
+            *a = 0;
+         }
+
+         arg0 = astFree( (void *) arg0 );
+         arg1 = astFree( (void *) arg1 );
+      }
+
+/* MULT... the second argument never needs to be in parentheses. The first
+   argument only needs to be placed in parentheses if it is a DIV or POW
+   node. */
+   } else if( newtree->opcode ==  OP_MULT ) {
+      if( mathmap == 2 && ( sunit = ModifyPrefix( newtree, status ) ) ) {
+         result = (char *) MakeExp( sunit, mathmap, 0, status );
+         sunit = FreeTree( sunit, status );
+
+      } else {
+         arg0 = MakeExp( newtree->arg[ 0 ], mathmap, 0, status );
+         larg0 = strlen( arg0 );
+
+         arg1 = MakeExp( newtree->arg[ 1 ], mathmap, 0, status );
+         larg1 = strlen( arg1 );
+
+/* If this is a top-level entry and we are producing an axis label, do
+   not include any constant multiplicative terms. */
+         if( top && !mathmap ) {
+            if( newtree->arg[ 0 ]->con != AST__BAD ) arg0 = astFree( (void *) arg0 );
+            if( newtree->arg[ 1 ]->con != AST__BAD ) arg1 = astFree( (void *) arg1 );
+         }
+
+/* If we have two arguments, concatentate them, placing the operands in
+   parentheses if necessary. */
+         if( arg0 && arg1 ) {
+
+            if( ( newtree->arg[ 0 ]->opcode == OP_DIV &&
+                  strchr( arg0, '/' ) ) ||
+                ( newtree->arg[ 0 ]->opcode == OP_POW &&
+                  strstr( arg0, "**" ) ) ) {
+               par = 1;
+               result = astMalloc( larg0 + larg1 + 4 );
+               if( result ) result[ 0 ] = '(';
+               a = result + 1;
+            } else {
+               par = 0;
+               result = astMalloc( larg0 + larg1 + 2 );
+               a = result;
+            }
+
+            if( result ) {
+               memcpy( a, arg0, larg0 );
+               a += larg0;
+               if( par ) *(a++) = ')';
+               *(a++) = '*';
+               memcpy( a, arg1, larg1 );
+               a += larg1;
+               *a = 0;
+            }
+
+            arg0 = astFree( (void *) arg0 );
+            arg1 = astFree( (void *) arg1 );
+
+/* If we do not have two arguments, just return the one we do have. */
+         } else if( arg0 ){
+            result = (char *) arg0;
+
+         } else {
+            result = (char *) arg1;
+         }
+      }
+   }
+
+/* Free the complicated tree. */
+   newtree = FreeTree( newtree, status );
+
+/* Free the returned string if an error has occurred. */
+   if( !astOK ) result = astFree( result );
+
+/* Return the result. */
+   return (const char *) result;
+}
+
+static void MakeKnownUnit( const char *sym, const char *label, const char *exp, int *status ){
+/*
+*  Name:
+*     MakeKnownUnit
+
+*  Purpose:
+*     Create a KnownUnit structure describing a known unit.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     void MakeKnownUnit( const char *sym, const char *label, const char *exp, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function creates a KnownUnit structure decribing a known unit,
+*     and adds it to the head of the linked list of known units stored in
+*     a module variable.
+
+*  Parameters:
+*     sym
+*        A pointer to a string which can be used as a symbol to represent
+*        the new named unit. Once defined, this symbol can be included within
+*        the definition of other derived units. The string should contain
+*        only alphabetical characters (no digits, spaces, punctuation,
+*        etc). Symbols are case sensitive (e.g. "s" is second, but "S" is
+*        Siemens). The string should not include any multiplier prefix.
+*     label
+*        Pointer to a null terminated string containing the label for
+*        the required units. No restriction on content.
+*     exp
+*        This should be a pointer to a null terminated string containing
+*        a definition of the required unit. See the description of the
+*        "in" and "out" parameters for the astUnitMapper function.
+*
+*        A NULL pointer or a blank string may supplied for "exp", which
+*        is interpreted as a request for a new basic unit to be created with
+*        the symbol and label given by the other parameters.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Notes:
+*     -  The supplied symbol and label strings are not copied. The
+*     supplied pointers are simply stored in the returned structure.
+*     Therefore the strings to which the pointers point should not be
+*     modified after this function returned (in fact this function is
+*     always called with literal strings for these arguments).
+*/
+
+/* Local Variables: */
+   KnownUnit *result;
+
+/* Check the global error status. */
+   if( !astOK ) return;
+
+/* Indicate that subsequent memory allocations may never be freed (other
+   than by any AST exit handler). */
+   astBeginPM;
+
+/* Allocate memory for the structure, and check the returned pointer can
+   be used safely. */
+   result = astMalloc( sizeof( KnownUnit ) );
+   if( astOK ) {
+
+/* In case of errors, first nullify the pointer to the next KnownUnit. */
+      result->next = NULL;
+
+/* Store the supplied label and symbol pointers. */
+      result->sym = sym;
+      result->label = label;
+
+/* Store the length of the symbol (without the trailing null character). */
+      result->symlen = strlen( sym );
+
+/* Store the length of the label (without the trailing null character). */
+      result->lablen = strlen( label );
+
+/* Create a tree of UnitNodes describing the unit if an expression was
+   supplied. */
+      result->head = exp ? CreateTree( exp, 1, 0, status ) : NULL;
+
+/* Unit aliases are replaced in use by the KnownUnit pointed to by the
+   "use" component of the structure. Indicate this KnownUnitis not an
+    alias by setting its "use" component NULL. */
+      result->use = NULL;
+   }
+
+/* Mark the end of the section in which memory allocations may never be
+   freed (other than by any AST exit handler). */
+   astEndPM;
+
+/* If an error has occurred, free any returned structure. */
+   if( !astOK ) {
+      result->head = FreeTree( result->head, status );
+      result = astFree( result ) ;
+
+/* Otherwise, add the new KnownUnit to the head of the linked list of
+   known units. */
+   } else {
+      result->next = known_units;
+      known_units = result;
+   }
+
+}
+
+static AstMapping *MakeMapping( UnitNode *tree, int *status ) {
+/*
+*  Name:
+*     MakeMapping
+
+*  Purpose:
+*     Create a new Mapping from a given tree of UnitNodes.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     AstMapping *MakeMapping( UnitNode *tree )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function creates a Mapping with a forward transformation equal
+*     to the transformation described by the tree of UnitNodes. The head
+*     node of the tree corresponds to the output of the Mapping.
+
+*  Parameters:
+*     tree
+*        The UnitNode at the head of the tree to be used. It should have
+*        exactly one OP_LDVAR node, and should have been simplified using
+*        the SimplifyTree function.
+
+*  Returned Value:
+*     A pointer to the Mapping. Its Nin and Nout attributes will both be 1.
+
+*  Notes:
+*     - A value of NULL will be returned if this function is invoked with
+*     the global error status set, or if it should fail for any reason.
+*/
+
+/* Local Variables: */
+   AstMapping *result;
+   UnitNode *inv;
+   UnitNode *src;
+   const char *fwdexp;
+   char *fwdfun;
+   const char *invexp;
+   char *invfun;
+   int lfwd;
+   int linv;
+
+/* Initialise. */
+   result = NULL;
+
+/* Check the inherited status. */
+   if( !astOK ) return result;
+
+/* First see if a UnitMap can be used to represent the Mapping from
+   input units to output units. This will be the case if the supplied tree
+   consists of a aingle OP_LDVAR node (corresponding to the input units). */
+   if( tree->opcode == OP_LDVAR ) {
+      result = (AstMapping *) astUnitMap( 1, "", status );
+
+/* Now see if a UnitMap or ZoomMap can be used to represent the Mapping from
+   input units to output units. This will be the case if the supplied tree
+   consists of a OP_MULT node with one constant argument and on OP_LDVAR
+   argument (corresponding to the input units). The standardisation done by
+   SimplifyTree will have ensured that the constant will be argument 0
+   (and will also have converted "x/k" trees into "(1/k)*x" trees). */
+   } else if( tree->opcode == OP_MULT &&
+              tree->arg[ 0 ]->con != AST__BAD &&
+              tree->arg[ 1 ]->opcode == OP_LDVAR ) {
+
+      if( tree->arg[ 0 ]->con == 1.0 ) {
+         result = (AstMapping *) astUnitMap( 1, "", status );
+      } else {
+         result = (AstMapping *) astZoomMap( 1, tree->arg[ 0 ]->con, "", status );
+      }
+
+/* For other trees we need to create a MathMap. */
+   } else {
+
+/* Format the supplied tree as an algebraic expression, and get its length. */
+      fwdexp = MakeExp( tree, 1, 1, status );
+      lfwd = strlen( fwdexp );
+
+/* The MathMap constructor requires the forward and inverse
+   transformation functions to be specified as equations (i.e. including an
+   equals sign). We use the output variable name "output_units" (the
+   astUnitMapper function creates the supplied tree usign the variable
+   name "input_units" ). */
+      lfwd += 13;
+
+/* Invert the supplied tree and create an algebraic expression from it. */
+      src = NewNode( NULL, OP_LDVAR, status );
+      if( astOK ) src->name = astStore( NULL, "output_units", 13 );
+      inv = InvertTree( tree, src, status );
+      if( !inv ) {
+         src = FreeTree( src, status );
+         astError( AST__BADUN, "MakeMapping(Unit): Failed to invert "
+                   "supplied tree '%s' (internal AST programming error).", status,
+                   fwdexp );
+
+/* If inverted succesfully (which it should be since astUnitMapper should
+   have checked this)... */
+      } else {
+
+/* Format the inverted tree as an algebraic expression, and get its
+   length, adding on extra characters for the variable name ("input_units")
+   and equals sign. */
+         invexp = MakeExp( inv, 1, 1, status );
+         linv = strlen( invexp );
+         linv += 12;
+
+/* Allocate memory for the transformation functions, plus an extra
+   character for the trailing null. */
+         fwdfun = astMalloc( lfwd + 1 );
+         invfun = astMalloc( linv + 1 );
+         if( invfun ) {
+            memcpy( fwdfun, "output_units=", 14 );
+            memcpy( invfun, "input_units=", 13 );
+
+/* Append the expressions following the equals signs. */
+            strcpy( fwdfun + 13, fwdexp );
+            strcpy( invfun + 12, invexp );
+
+/* Create the MathMap. */
+            result = (AstMapping *) astMathMap( 1, 1, 1,
+                                                (const char **) &fwdfun, 1,
+                                                (const char **) &invfun,
+                                                "SimpFI=1,SimpIF=1", status );
+         }
+
+/* Free resources. */
+         inv = FreeTree( inv, status );
+         fwdfun = astFree( fwdfun );
+         invfun = astFree( invfun );
+         invexp = astFree( (void *) invexp );
+      }
+      fwdexp = astFree( (void *) fwdexp );
+   }
+
+/* Free any result if an error occurred. */
+   if( !astOK ) result = astAnnul( result );
+
+/* Return the answer. */
+   return result;
+}
+
+static UnitNode *MakeLabelTree( const char *lab, int nc, int *status ){
+/*
+*  Name:
+*     MakeLabelTree
+
+*  Purpose:
+*     Convert an axis label into a tree of UnitNodes.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     UnitNode *MakeLabelTree( const char *lab, int nc, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function converts an axis label into a tree of UnitNodes.
+*     It is assumed the supplied label represents some "basic" label
+*     modified by the application of one or more single function arguments
+*     and/or exponentiation operators. The (single) OP_LDVAR node in the
+*     returned tree refers to the basic label (it is stored as the "name"
+*     component of UnitNode structure).
+
+*  Parameters:
+*     lab
+*        The label expression.
+*     nc
+*        The number of characters from "lab" to use.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     A pointer to a UnitNode which forms the head of a tree of UnitNodes
+*     representing the supplied label expression.
+
+*  Notes:
+*     -  A NULL value is returned if this function is invoked with the
+*     global error status set or if it should fail for any reason.
+*/
+
+/* Local Variables: */
+   Oper op;
+   UnitNode *result;
+   char buff[ 10 ];
+   const char *c;
+   const char *exp;
+   int depth;
+   int i;
+   int oplen;
+   int n;
+   double con;
+
+/* Initialise */
+   result = NULL;
+   oplen = 0;
+
+/* Check the global error status, and that we have a string. */
+   if ( !astOK || !lab || !nc ) return result;
+
+/* Get a pointer to the first non-blank character, and store the number of
+   characters to examine (this excludes any trailing white space). */
+   exp = lab;
+   while( isspace( *exp ) ) exp++;
+   c = lab + nc - 1;
+   while( c >= exp && isspace( *c ) ) c--;
+   nc = c - exp + 1;
+
+/* Scan through the supplied string looking for the first pow operator at
+   zero depth of nesting within parentheses. */
+   depth = 0;
+   c = exp;
+   i = 0;
+   op = OP_NULL;
+   while( i < nc && *c ){
+
+/* If this character is an opening parenthesis, increment the depth of
+   nesting. */
+      if( *c == '(' ) {
+         depth++;
+
+/* If this character is an closing parenthesis, decrement the depth of
+   nesting. Report an error if it ever goes negative. */
+      } else if( *c == ')' ) {
+         depth--;
+         if( depth < 0 && astOK ) {
+            astError( AST__BADUN, "Missing opening parenthesis." , status);
+            break;
+         }
+
+/* Ignore all other characters unless they are at zero depth of nesting.
+   Also ignore spaces. */
+      } else if( depth == 0 && !isspace( *c ) ) {
+
+/* Compare the next part of the string with each of the "pow" operators. */
+         if( !strncmp( c, "**", 2 ) ) {
+            op = OP_POW;
+            oplen = 2;
+         } else if( *c == '^' ) {
+            op = OP_POW;
+            oplen = 1;
+         }
+
+/* If an operator was found, break out of the loop. */
+         if( op != OP_NULL ) break;
+      }
+
+/* Pass on to check the next character. */
+      i++;
+      c++;
+   }
+
+/* If a "pow" operator was found, the strings on either side of it should be
+   valid unit expressions, in which case we use this routine recursively to
+   create corresponding trees of UnitNodes. */
+   if( op != OP_NULL ) {
+
+/* Create a UnitNode for the operator. */
+      result = NewNode( NULL, op, status );
+      if( astOK ) {
+
+/* Create a tree of unit nodes from the string which precedes the binary
+   operator. Report an error if it cannot be done. */
+        result->arg[ 0 ] = MakeLabelTree( exp, i, status );
+        if( !result->arg[ 0 ] && astOK ) {
+           for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
+           buff[ oplen ] = 0;
+           astError( AST__BADUN, "Missing operand before '%s'.", status, buff );
+        }
+
+/* Create a tree of unit nodes from the string which follows the binary
+   operator. Report an error if it cannot be done. */
+         result->arg[ 1 ] = MakeLabelTree( c + oplen, nc - i - oplen, status );
+         if( !result->arg[ 1 ] && astOK ) {
+             for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
+            buff[ oplen ] = 0;
+            astError( AST__BADUN, "Missing operand after '%s'.", status, buff );
+         }
+      }
+
+/* If no binary operator was found at depth zero, see if the supplied string
+   starts with a function name (the only legal place for a function name
+   given that the string has no binary operators at depth zero). */
+   } else {
+      if( !strncmp( exp, "sqrt(", 5 ) || !strncmp( exp, "SQRT(", 5 ) ) {
+         op = OP_SQRT;
+         oplen = 4;
+      } else if( !strncmp( exp, "exp(", 4 ) || !strncmp( exp, "EXP(", 4 ) ) {
+         op = OP_EXP;
+         oplen = 3;
+      } else if( !strncmp( exp, "ln(", 3 ) || !strncmp( exp, "LN(", 3 ) ) {
+         op = OP_LN;
+         oplen = 2;
+      } else if( !strncmp( exp, "log(", 4 ) || !strncmp( exp, "LOG(", 4 ) ) {
+         op = OP_LOG;
+         oplen = 3;
+      }
+
+/* If a function was found, the string following the function name
+   (including the opening parenthesis) should form a legal units
+   expresssion (all the supported functions take a single argument and
+   so we do not need to worry about comma-separated lists of function
+   arguments). Use this routine recursively to create a tree of UnitNodes
+   from the string which forms the function argument. */
+      if( op != OP_NULL ) {
+
+/* Create a UnitNode for the function. */
+         result = NewNode( NULL, op, status );
+         if( astOK ) {
+
+/* Create a tree of unit nodes from the string which follows the function
+   name. Report an error if it cannot be done. */
+            result->arg[ 0 ] = MakeLabelTree( exp + oplen, nc - oplen, status );
+            if( !result->arg[ 0 ] && astOK ) {
+               for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
+               buff[ oplen ] = 0;
+               astError( AST__BADUN, "Missing argument for '%s'.", status, buff );
+            }
+         }
+
+/* Arrive here if the supplied string does not contain a POW operator
+   or function at depth zero. Check to see if the whole string is contained
+   within parentheses, In which we interpret the contents of the
+   parentheses as a units expression. It is safe simply to check the
+   first and last characters (a string like "(fred)(Harry)" is not a
+   legal possibility since there should be an operator in the middle).*/
+      } else if( nc > 0 && ( exp[ 0 ] == '(' && exp[ nc - 1 ] == ')' ) ) {
+         result = MakeLabelTree( exp + 1, nc - 2, status );
+
+/* Does the string begin with a numerical constant? */
+      } else if( ConStart( exp, &con, &n, status ) == 1 ) {
+
+/* If the entire string was a numerical constant, represent it by a LDCON
+   node. */
+         if( n == nc ) {
+            result = NewNode( NULL, OP_LDCON, status );
+            if( astOK ) result->con = con;
+
+/* If there was anything following the numerical constant, report an
+   error. */
+         } else if( astOK ){
+            astError( AST__BADUN, "Missing operator after "
+                      "numerical string '%.*s'.", status, n, exp );
+         }
+
+/* The only legal possibility left is that the string represents the basic
+   label. Create an OP_LDVAR node for it and store the basic label as
+   the node name, omitting any enclosing white space. */
+      } else {
+         result = NewNode( NULL, OP_LDVAR, status );
+         if( astOK ) {
+            result->name = astStore( NULL, exp, nc + 1 );
+            if( astOK ) ( (char *) result->name)[ nc ] = 0;
+         }
+      }
+   }
+
+/* Free any returned tree if an error has occurred. */
+   if( !astOK ) result = FreeTree( result, status );
+
+/* Return the result. */
+   return result;
+}
+
+static UnitNode *MakeTree( const char *exp, int nc, int lock, int *status ){
+/*
+*  Name:
+*     MakeTree
+
+*  Purpose:
+*     Convert an algebraic units expression into a tree of UnitNodes.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     UnitNode *MakeTree( const char *exp, int nc, int lock, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function converts an algebraic units expression into a tree of
+*     UnitNodes. It is a service routine for CreateTree. The roots of the
+*     returned tree (i.e. the LDVAR nodes) refer to the unit symbols
+*     contained within the supplied expression (i.e. definitions of these
+*     units are not grafted onto the tree in place of the original nodes,
+*     as is done by CreateTree).
+
+*  Parameters:
+*     exp
+*        The units expression. This should not include any leading or
+*        trailing spaces.
+*     nc
+*        The number of characters from "exp" to use.
+*     lock
+*        Use a mutex to guard access to the KnownUnits list?
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     A pointer to a UnitNode which forms the head of a tree of UnitNodes
+*     representing the supplied unit expression.
+
+*  Notes:
+*     -  A NULL value is returned if this function is invoked with the
+*     global error status set or if it should fail for any reason.
+*/
+
+/* Local Variables: */
+   KnownUnit *munit;
+   KnownUnit *unit;
+   Multiplier *mmult;
+   Multiplier *mult;
+   Oper op;
+   UnitNode *result;
+   char buff[ 10 ];
+   char d;
+   const char *c;
+   double con;
+   int depth;
+   int i;
+   int l;
+   int maxlen;
+   int n;
+   int oplen;
+   int plural;
+
+/* Initialise */
+   result = NULL;
+
+/* Check the global error status, and that we have a string. */
+   if ( !astOK || !exp || nc <= 0 ) return result;
+
+/* Scan through the supplied string from the end to the start looking for
+   the last multiplication or division operator at zero depth of nesting
+   within parentheses. We go backwards through the string in order to
+   give the correct priority to multiple division operators (i.e. "a/b/c"
+   needs to be interpreted as "(a/b)/c", not "a/(b/c)"). */
+   op = OP_NULL;
+   oplen = 1;
+   depth = 0;
+   c = exp + nc - 1;
+   i = nc - 1;
+   while( i >= 0 ){
+
+/* If this character is an opening parenthesis, decrement the depth of
+   nesting. Report an error if it ever goes negative. */
+      if( *c == '(' ) {
+         depth--;
+         if( depth < 0 && astOK ) {
+            astError( AST__BADUN, "Missing closing parenthesis." , status);
+            break;
+         }
+
+/* An opening parenthesis at level zero must always be either the first
+   character in the string, or be preceded by the name of a function, or
+   be preceded by an operator. If none of these are true, assume there is
+   an implicit multiplication operator before the parenthesis. */
+         if( depth == 0 && i > 0 ) {
+            d = *( c - 1 );
+            if( d != '*' && d != '/' && d != '^' && d != '.' && d != ' ' &&
+                !EndsWith( c, i + 1, "sqrt(", status ) && !EndsWith( c, i + 1, "exp(", status ) &&
+                !EndsWith( c, i + 1, "ln(", status ) && !EndsWith( c, i + 1, "log(", status ) ) {
+               op = OP_MULT;
+               oplen = 0;
+               break;
+            }
+         }
+
+/* If this character is an closing parenthesis, increment the depth of
+   nesting. */
+      } else if( *c == ')' ) {
+         depth++;
+
+/* A closing parenthesis at level zero must always be either the last
+   character in the string, or be followed by an operator. If neither of
+   these are true, assume there is an implicit multiplication operator. */
+         if( depth == 1 && i < nc - 1 ) {
+            d = *(c+1);
+            if( d != '*' && d != '/' && d != '^' && d != '.' && d != ' ') {
+               op = OP_MULT;
+               oplen = 0;
+
+/* Correct "i" so that it gives the length of the left hand operand of
+   the implicit MULT operator, correct "c" so that it points to the first
+   character in the right hand operand, and leave the loop. */
+               i++;
+               c++;
+               break;
+            }
+         }
+
+/* Ignore all other characters unless they are at zero depth of nesting. */
+      } else if( depth == 0 ) {
+
+/* Compare the next part of the string with each of the multiplication
+   and division operators. */
+         if( *c == '/' ) {
+            op = OP_DIV;
+
+         } else if( *c == ' ' ) {
+            op = OP_MULT;
+
+/* An asterisk is only treated as a multiplication symbol if it does not occur
+   before or after another asterisk. */
+         } else if( *c == '*' ) {
+            if(  c == exp ) {
+               if( *(c+1) != '*' ) op = OP_MULT;
+            } else if( i == nc - 1 ) {
+               if( *(c-1) != '*' ) op = OP_MULT;
+            } else {
+               if( *(c+1) != '*' && *(c-1) != '*' ) op = OP_MULT;
+            }
+
+/* A dot is only treated as a multiplication symbol if it does not occur
+   between two digits. */
+         } else if( *c == '.' ) {
+            if( ( c == exp || !isdigit( *(c-1) ) ) &&
+                ( i == nc - 1 || !isdigit( *(c+1) ) ) ) {
+               op = OP_MULT;
+            }
+         }
+      }
+
+/* If an operator was found, break out of the loop. */
+      if( op != OP_NULL ) break;
+
+/* Pass on to check the next character. */
+      i--;
+      c--;
+   }
+
+/* If a multiplication or division operator was found, the strings on either
+   side of it should be valid unit expressions, in which case we use this
+   routine recursively to create corresponding trees of UnitNodes. */
+   if( op != OP_NULL ) {
+
+/* Create a UnitNode for the binary operator. */
+      result = NewNode( NULL, op, status );
+      if( astOK ) {
+
+/* Create a tree of unit nodes from the string which precedes the binary
+   operator. Report an error if it cannot be done. */
+         result->arg[ 0 ] = MakeTree( exp, i, lock, status );
+         if( !result->arg[ 0 ] && astOK ) {
+            for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
+            buff[ oplen ] = 0;
+            astError( AST__BADUN, "Missing operand before '%s'.", status, buff );
+         }
+
+/* Create a tree of unit nodes from the string which follows the binary
+   operator. Report an error if it cannot be done. */
+         result->arg[ 1 ] = MakeTree( c + oplen, nc - i - oplen, lock, status );
+         if( !result->arg[ 1 ] && astOK ) {
+            for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
+            buff[ oplen ] = 0;
+            astError( AST__BADUN, "Missing operand after '%s'.", status, buff );
+         }
+      }
+
+/* If no multiplication or division operator was found at depth zero, check
+   that the final depth of nesting was zero. Report an error if not. */
+   } else if( depth > 0 && astOK ) {
+      astError( AST__BADUN, "Missing opening parenthesis." , status);
+
+/* Otherwise check for a "Pow" operator at depth zero. */
+   } else {
+
+/* Scan through the supplied string looking for the first pow operator at
+   zero depth of nesting within parentheses. */
+      depth = 0;
+      c = exp;
+      i = 0;
+      while( i < nc && *c ){
+
+/* If this character is an opening parenthesis, increment the depth of
+   nesting. */
+         if( *c == '(' ) {
+            depth++;
+
+/* If this character is an closing parenthesis, decrement the depth of
+   nesting. Report an error if it ever goes negative. */
+         } else if( *c == ')' ) {
+            depth--;
+            if( depth < 0 && astOK ) {
+               astError( AST__BADUN, "Missing opening parenthesis." , status);
+               break;
+            }
+
+/* Ignore all other characters unless they are at zero depth of nesting. */
+         } else if( depth == 0 ) {
+
+/* Compare the next part of the string with each of the "pow" operators. */
+            if( !strncmp( c, "**", 2 ) ) {
+               op = OP_POW;
+               oplen = 2;
+            } else if( *c == '^' ) {
+               op = OP_POW;
+               oplen = 1;
+            }
+
+/* If an operator was found, break out of the loop. */
+            if( op != OP_NULL ) break;
+         }
+
+/* Pass on to check the next character. */
+         i++;
+         c++;
+      }
+
+/* If a "pow" operator was found, the strings on either side of it should be
+   valid unit expressions, in which case we use this routine recursively to
+   create corresponding trees of UnitNodes. */
+      if( op != OP_NULL ) {
+
+/* Create a UnitNode for the operator. */
+         result = NewNode( NULL, op, status );
+         if( astOK ) {
+
+/* Create a tree of unit nodes from the string which precedes the binary
+   operator. Report an error if it cannot be done. */
+            result->arg[ 0 ] = MakeTree( exp, i, lock, status );
+            if( !result->arg[ 0 ] && astOK ) {
+               for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
+               buff[ oplen ] = 0;
+               astError( AST__BADUN, "Missing operand before '%s'.", status, buff );
+            }
+
+/* Create a tree of unit nodes from the string which follows the binary
+   operator. Report an error if it cannot be done. */
+            result->arg[ 1 ] = MakeTree( c + oplen, nc - i - oplen, lock, status );
+            if( !result->arg[ 1 ] && astOK ) {
+                for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
+               buff[ oplen ] = 0;
+               astError( AST__BADUN, "Missing operand after '%s'.", status, buff );
+            }
+         }
+
+/* If no binary operator was found at depth zero, see if the supplied string
+   starts with a function name (the only legal place for a function name
+   given that the string has no binary operators at depth zero). */
+      } else {
+         if( !strncmp( exp, "sqrt(", 5 ) || !strncmp( exp, "SQRT(", 5 ) ) {
+            op = OP_SQRT;
+            oplen = 4;
+         } else if( !strncmp( exp, "exp(", 4 ) || !strncmp( exp, "EXP(", 4 ) ) {
+            op = OP_EXP;
+            oplen = 3;
+         } else if( !strncmp( exp, "ln(", 3 ) || !strncmp( exp, "LN(", 3 ) ) {
+            op = OP_LN;
+            oplen = 2;
+         } else if( !strncmp( exp, "log(", 4 ) || !strncmp( exp, "LOG(", 4 ) ) {
+            op = OP_LOG;
+            oplen = 3;
+         }
+
+/* If a function was found, the string following the function name
+   (including the opening parenthesis) should form a legal units
+   expresssion (all the supported functions take a single argument and
+   so we do not need to worry about comma-separated lists of function
+   arguments). Use this routine recursively to create a tree of UnitNodes
+   from the string which forms the function argument. */
+         if( op != OP_NULL ) {
+
+/* Create a UnitNode for the function. */
+            result = NewNode( NULL, op, status );
+            if( astOK ) {
+
+/* Create a tree of unit nodes from the string which follows the function
+   name. Report an error if it cannot be done. */
+               result->arg[ 0 ] = MakeTree( exp + oplen, nc - oplen, lock, status );
+               if( !result->arg[ 0 ] && astOK ) {
+                  for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
+                  buff[ oplen ] = 0;
+                  astError( AST__BADUN, "Missing argument for '%s'.", status, buff );
+               }
+            }
+
+/* Arrive here if the supplied string does not contain a binary operator
+   or function at depth zero. Check to see if the whole string is contained
+   within parentheses, In which we interpret the contents of the
+   parentheses as a units expression. It is safe simply to check the
+   first and last characters (a string like "(fred)(Harry)" is not a
+   legal possibility since there should be an operator in the middle).*/
+         } else if( exp[ 0 ] == '(' && exp[ nc - 1 ] == ')' ) {
+            result = MakeTree( exp + 1, nc - 2, lock, status );
+
+/* Does the string begin with a numerical constant? */
+         } else if( ConStart( exp, &con, &n, status ) == 1 ) {
+
+/* If the entire string was a numerical constant, represent it by a LDCON
+   node. */
+            if( n == nc ) {
+               result = NewNode( NULL, OP_LDCON, status );
+               if( astOK ) result->con = con;
+
+/* If there was anything following the numerical constant, report an
+   error. */
+            } else if( astOK ){
+               astError( AST__BADUN, "Missing operator after "
+                         "numerical string '%.*s'.", status, n, exp );
+            }
+
+/* Does the string represent one of the named constants? If so represent it
+   by a an appropriate operator. */
+         } else if( nc == 2 && ( !strncmp( exp, "pi", 2 ) ||
+                                 !strncmp( exp, "PI", 2 ) ) ) {
+            result = NewNode( NULL, OP_LDPI, status );
+
+         } else if( nc == 1 && ( !strncmp( exp, "e", 1 ) ||
+                                 !strncmp( exp, "E", 1 ) ) ) {
+            result = NewNode( NULL, OP_LDE, status );
+
+/* The only legal possibility left is that the string represents the name
+   of a basic unit, possibly prefixed by a multiplier character. */
+         } else {
+
+/* See if the string ends with the symbol for any of the known basic
+   units. If it matches more than one basic unit, choose the longest.
+   First ensure descriptions of the known units are  available. */
+            mmult = NULL;
+            plural = 0;
+            while( 1 ) {
+               unit = GetKnownUnits( lock, status );
+
+               maxlen = -1;
+               munit = NULL;
+               while( unit ) {
+                  if( SplitUnit( exp, nc, unit->sym, 1, &mult, &l, status ) ) {
+                     if( l > maxlen ) {
+                        maxlen = l;
+                        munit = unit;
+                        mmult = mult;
+                     }
+                  }
+                  unit = unit->next;
+               }
+
+/* If the above did not produce a match, try matching the unit symbol
+   case insensitive. */
+               if( !munit ) {
+                  unit = GetKnownUnits( lock, status );
+                  while( unit ) {
+                     if( SplitUnit( exp, nc, unit->sym, 0, &mult, &l, status ) ) {
+                        if( l > maxlen ) {
+                           maxlen = l;
+                           munit = unit;
+                           mmult = mult;
+                        }
+                     }
+                     unit = unit->next;
+                  }
+               }
+
+/* If the above did not produce a match, try matching the unit label
+   case insensitive. */
+               if( !munit ) {
+                  unit = GetKnownUnits( lock, status );
+                  while( unit ) {
+                     if( SplitUnit( exp, nc, unit->label, 0, &mult, &l, status ) ) {
+                        if( l > maxlen ) {
+                           maxlen = l;
+                           munit = unit;
+                           mmult = mult;
+                        }
+                     }
+                     unit = unit->next;
+                  }
+               }
+
+/* If we still do not have a match, and if the string ends with "s", try
+   removing the "s" (which could be a plural as in "Angstroms") and
+   trying again. */
+               if( !munit && nc > 1 && !plural &&
+                   ( exp[ nc - 1 ] == 's' || exp[ nc - 1 ] == 'S' ) ) {
+                  plural = 1;
+                  nc--;
+               } else {
+                  break;
+               }
+            }
+            if( plural ) nc++;
+
+/* If a known unit and multiplier combination was found, create an
+   OP_LDVAR node from it. */
+            unit = munit;
+            mult = mmult;
+            if( unit ) {
+
+/* If the unit is an alias for another unit, it will have a non-NULL
+   value for its "use" component.In this case, use the unit for which the
+   identified unit is an alias. */
+               result = NewNode( NULL, OP_LDVAR, status );
+               if( astOK ) {
+                  result->unit = unit->use ? unit->use : unit;
+                  result->mult = mult;
+                  result->name = astStore( NULL, result->unit->sym, result->unit->symlen + 1 );
+               }
+
+/* If no known unit and multiplier combination was found, we assume the
+   string represents a new user-defined basic unit, possibly preceded by a
+   standard multiplier prefix. */
+            } else {
+
+/* Check the string to see if starts with a known multiplier prefix (but
+   do not allow the multiplier to account for the entire string). */
+               mult = GetMultipliers( status );
+               c = exp;
+               while( mult ) {
+                  n = nc - mult->symlen;
+                  if( n > 0 && !strncmp( exp, mult->sym, mult->symlen ) ) {
+                     c += mult->symlen;
+                     break;
+                  }
+                  mult = mult->next;
+               }
+               if( !mult ) n = nc;
+
+/* Check there are no illegal characters in the following string. */
+               for( i = 0; i < n && astOK; i++ ) {
+                  if( !isalpha( c[ i ] ) ) {
+                     astError( AST__BADUN, "Illegal character '%c' found.", status, c[ i ] );
+                     break;
+                  }
+               }
+
+/* If succesfull, create an OP_LDVAR node for th user-defined basic unit. */
+               if( astOK ) {
+                  result = NewNode( NULL, OP_LDVAR, status );
+                  if( astOK ) {
+                     result->mult = mult;
+                     result->name = astStore( NULL, c, n + 1 );
+                     if( astOK ) ( (char *) result->name)[ n ] = 0;
+                  }
+               }
+            }
+         }
+      }
+   }
+
+/* Free any returned tree if an error has occurred. */
+   if( !astOK ) result = FreeTree( result, status );
+
+/* Return the result. */
+   return result;
+}
+
+static void MakeUnitAlias( const char *sym, const char *alias, int *status ){
+/*
+*  Name:
+*     MakeUnitAlias
+
+*  Purpose:
+*     Create a KnownUnit structure describing an alias for a known unit.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     void MakeUnitAlias( const char *sym, const char *alias, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function creates a KnownUnit structure decribing an alias for a
+*     known unit, and adds it to the head of the linked list of known units
+*     stored in a module variable. An alias is a KnownUnit which is
+*     identical to an existing known but which has a different symbol.
+
+*  Parameters:
+*     sym
+*        A pointer to the symbol string of an existing KnwonUnit. The string
+*        should not include any multiplier prefix.
+*     alias
+*        A pointer to the symbol string to use as the alasi for the existing
+*        KnownUnit. The string should not include any multiplier prefix.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Notes:
+*     -  The supplied symbol and label strings are not copied. The
+*     supplied pointers are simply stored in the returned structure.
+*     Therefore the strings to which the pointers point should not be
+*     modified after this function returned (in fact this function is
+*     always called with literal strings for these arguments).
+*/
+
+/* Local Variables: */
+   KnownUnit *unit;
+
+/* Check the global error status. */
+   if( !astOK ) return;
+
+/* Search the existing list of KnownUnits for the specified symbol. */
+   unit = known_units;
+   while( unit ) {
+      if( !strcmp( sym, unit->sym ) ) {
+
+/* Create a new KnownUnit for the alias. It will becomes the head of the
+   known units chain. */
+         MakeKnownUnit( alias, unit->label, NULL, status );
+
+/* Store a pointer to the KnownUnit which is to be used in place of the
+   alias. */
+         known_units->use = unit;
+
+/* Leave the loop. */
+         break;
+      }
+
+/* Move on to check the next existing KnownUnit. */
+      unit = unit->next;
+   }
+
+/* Report an error if the supplied unit was not found. */
+   if( !unit ) {
+      astError( AST__INTER, "MakeUnitAlias(Unit): Cannot find existing "
+                "units \"%s\" to associate with the alias \"%s\" (AST "
+                "internal programming error).", status, sym, alias );
+   }
+}
+
+static UnitNode *ModifyPrefix( UnitNode *old, int *status ) {
+/*
+*  Name:
+*     ModifyPrefix
+
+*  Purpose:
+*     Replace a MULT or DIV node with a LDVAR and suitable multiplier.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     UnitNode *ModifyPrefix( UnitNode *old, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function checks the supplied node. If it is a DIV or MULT node
+*     in which one argument is an LDVAR and the other is a constant, then
+*     its checks to see if the constant can be absorbed into the LDVAR by
+*     changing the multiplier in the LDVAR node. If so, it returns a new
+*     node which is an LDVAR with the modified multiplier. Otherwise it
+*     returns NULL.
+
+*  Parameters:
+*     old
+*        Pointer to an existing UnitNode to be checked.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     A pointer to the new UnitNode.
+
+*  Notes:
+*     - A value of NULL will be returned if this function is invoked with
+*     the global error status set, or if it should fail for any reason.
+*/
+
+/* Local Variables: */
+   Multiplier *mult;
+   Multiplier *mmult;
+   UnitNode *ldcon;
+   UnitNode *ldvar;
+   UnitNode *newtree;
+   UnitNode *result;
+   double con;
+   double cmult;
+   double r;
+   double rmin;
+   int recip;
+   int changed;
+
+/* Initialise. */
+   result = NULL;
+
+/* Check the inherited status. */
+   if( !astOK ) return result;
+
+/* Indicate that we have not yet found any reason to return a changed
+   node. */
+   changed = 0;
+
+/* Check the supplied node is a DIV or MULT node. */
+   if( old->opcode == OP_DIV || old->opcode == OP_MULT ) {
+
+/* Get a copy of the supplied tree which we can modify safely. */
+      newtree = CopyTree( old, status );
+
+/* Identify the LDVAR argument (if any). */
+      if( newtree->arg[ 0 ]->opcode == OP_LDVAR ) {
+         ldvar = newtree->arg[ 0 ];
+
+      } else if( newtree->arg[ 1 ]->opcode == OP_LDVAR ) {
+         ldvar = newtree->arg[ 1 ];
+
+      } else {
+         ldvar = NULL;
+      }
+
+/* Identify the LDCON argument (if any). */
+      if( newtree->arg[ 0 ]->opcode == OP_LDCON ) {
+         ldcon = newtree->arg[ 0 ];
+
+      } else if( newtree->arg[ 1 ]->opcode == OP_LDCON ) {
+         ldcon = newtree->arg[ 1 ];
+
+      } else {
+         ldcon = NULL;
+      }
+
+/* If either was not found, return NULL. */
+      if( !ldvar || !ldcon ) {
+         newtree = FreeTree( newtree, status );
+
+/* Otherwise, extract the multiplier constant. If there is no multiplier, the
+   constant is 1.0. */
+      } else {
+         cmult = ldvar->mult ? ldvar->mult->scale: 1.0;
+
+/* Extract the constant. */
+         con = ldcon->con;
+
+/* Combine the multiplier and the constant. The resulting constant is a
+   factor which is used to multiply the LDVAR quantity. If the original
+   node is a DIV node in which the LDVAR is in the denominator, then
+   flag that we need to reciprocate the new MULT node which represents
+   "constant*LDVAR" before returning. */
+         if( newtree->opcode == OP_MULT ) {
+            con = con*cmult;
+            recip = 0;
+         } else {
+            con = cmult/con;
+            recip = ( ldvar == newtree->arg[ 1 ] );
+         }
+
+/* Find the closest known multiplier to the new constant. */
+         rmin = ( con > 1 ) ? con : 1.0/con;
+         mmult = NULL;
+         mult = GetMultipliers( status );
+         while( mult ) {
+            r = ( con > mult->scale) ? con/mult->scale : mult->scale/con;
+            if( r < rmin ) {
+               mmult = mult;
+               rmin = r;
+            }
+            mult = mult->next;
+         }
+
+/* Modify the constant to take account of the new multiplier chosen
+   above. "mmult" will be NULL if the best multiplier is unity. */
+         if( mmult ) con = con/mmult->scale;
+
+/* If they have changed, associate the chosen multiplier with the LDVAR node,
+   and the constant with the LDCON node. */
+         if( ldvar->mult != mmult ) {
+            ldvar->mult = mmult;
+            changed = 1;
+         }
+
+         if( ldcon->con != con ) {
+            ldcon->con = con;
+            changed = 1;
+         }
+
+/* Unless the node is proportional to the reciprocal of the variable, the
+   new node should be a MULT node (it may originally have been a DIV). */
+         if( !recip ) {
+            if( newtree->opcode != OP_MULT ){
+               newtree->opcode = OP_MULT;
+               changed = 1;
+            }
+
+/* If the constant is 1.0 we can just return the LDVAR node by itself. */
+            if( fabs( con - 1.0 ) < 1.0E-6 ) {
+               result = CopyTree( ldvar, status );
+               newtree = FreeTree( newtree, status );
+               changed = 1;
+
+/* Otherwise return the modified tree containing both LDVAR and LDCON nodes. */
+            } else {
+               result = newtree;
+            }
+
+/* If the node is proportional to the reciprocal of the variable, the
+   new node will already be a DIV node and will have an LDCON as the first
+   argument (numerator) and an LDVAR as the second argument (denominator). */
+         } else {
+
+/* The first argument (the numerator) should be the reciprocal of the constant
+   found above. */
+            ldcon->con = 1.0/ldcon->con;
+            if( !astEQUAL( ldcon->con, old->arg[0]->con ) ) changed = 1;
+
+/* Return the modified tree containing both LDVAR and LDCON nodes. */
+            result = newtree;
+         }
+      }
+   }
+
+/* If the new and old trees are equivalent, then we do not need to return
+   it. */
+   if( !changed && result ) result = FreeTree( result, status );
+
+/* Return the answer. */
+   return result;
+}
+
+
+static UnitNode *NewNode( UnitNode *old, Oper code, int *status ) {
+/*
+*  Name:
+*     NewNode
+
+*  Purpose:
+*     Create and initialise a new UnitNode.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     UnitNode *NewNode( UnitNode *old, Oper code, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function creates and initialises a new UnitNode, or
+*     re-initialises an existing UnitNode to use a different op code.
+
+*  Parameters:
+*     old
+*        Pointer to an existing UnitNode to be modified, or NULL to create
+*        a new UnitNode.
+*     code
+*        The op code for the new UnitNode.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     A pointer to the new UnitNode.
+
+*  Notes:
+*     - A value of NULL will be returned if this function is invoked with
+*     the global error status set, or if it should fail for any reason.
+*/
+
+/* Local Variables: */
+   UnitNode **args;
+   UnitNode *result;
+   int i;
+
+/* Initialise. */
+   result = NULL;
+   args = NULL;
+
+/* Check the inherited status. */
+   if( !astOK ) return result;
+
+/* If an existig UnitNode was supplied, free any memory used to hold
+   pointers to its arguments. */
+   if( old ) {
+      old->arg = astFree( old->arg );
+      result = old;
+
+/* Otherwise, allocate memory for a new structure. */
+   } else {
+      result = astMalloc( sizeof( UnitNode ) );
+   }
+
+/* Check the pointer can be used safely. */
+   if( astOK ) {
+
+/* Initialise the members of the UnitNode structure. */
+      result->opcode = code;
+      result->arg = NULL;
+      result->con = AST__BAD;
+      result->name = NULL;
+      result->unit = NULL;
+      result->mult = NULL;
+      result->narg = 0;
+
+      switch( code ){
+         case OP_LDPI:
+            result->con = PI;
+            break;
+
+         case OP_LDE:
+            result->con = E;
+            break;
+
+         case OP_LOG:
+         case OP_LN:
+         case OP_EXP:
+         case OP_SQRT:
+            result->narg = 1;
+            break;
+
+         case OP_POW:
+         case OP_DIV:
+         case OP_MULT:
+            result->narg = 2;
+            break;
+
+         default:
+            ;
+      }
+
+/* Allocate memory for the UnitNode pointers which will locate the
+   nodes forming the arguments to the new node. */
+      args = astMalloc( (result->narg)*sizeof( UnitNode * ) );
+      if( astOK ) {
+         result->arg = args;
+
+/* Initialise the argument pointers to NULL. */
+         for( i = 0; i < result->narg; i++ ) args[ i ] = NULL;
+      }
+   }
+
+/* Free any result if an error occurred. */
+   if( !astOK ) {
+      args = astFree( args );
+      result = astFree( result );
+   }
+
+/* Return the answer. */
+   return result;
+}
+
+static void RemakeTree( UnitNode **node, int *status ) {
+/*
+*  Name:
+*     RemakeTree
+
+*  Purpose:
+*     Replace derived units within a tree of UnitNodes by basic units.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     void RemakeTree( UnitNode **node, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function searches for LDVAR nodes (i.e. references to unit
+*     symbols) within the given tree, and replaces each such node which
+*     refers to known derived unit with a sub-tree of nodes which
+*     define the derived unit in terms of known basic units.
+
+*  Parameters:
+*     node
+*        The address of a pointer to the UnitNode at the head of the tree
+*        which is to be simplified. On exit the supplied tree is freed and a
+*        pointer to a new tree is placed at the given address.
+*     status
+*        Pointer to the inherited status variable.
+
+*/
+
+/* Local Variables: */
+   KnownUnit *unit;
+   int i;
+   UnitNode *newnode;
+
+/* Check inherited status. */
+   if( !astOK ) return;
+
+/* Initially, we have no replacement node */
+   newnode = NULL;
+
+/* If this is an LDVAR node... */
+   if( (*node)->opcode == OP_LDVAR ) {
+
+/* If the LDVAR node has a multiplier associated with it, we need to
+   introduce a OP_MULT node to perform the scaling. */
+      if( (*node)->mult ) {
+         newnode = NewNode( NULL, OP_MULT, status );
+         if( astOK ) {
+            newnode->arg[0] = NewNode( NULL, OP_LDCON, status );
+            if( astOK ) {
+               newnode->arg[0]->con = 1.0/(*node)->mult->scale;
+
+/* See if the node refers to a known unit. If not, or if the known unit
+   is a basic unit (i.e. not a derived unit) use the supplied node for
+   the second argument of the OP_MULT node (without the multiplier).
+   Otherwise, use a copy of the tree which defines the derived unit. */
+               unit = (*node)->unit;
+               if( unit && unit->head ) {
+                  newnode->arg[1] = CopyTree( unit->head, status );
+               } else {
+                  newnode->arg[1] = CopyTree( *node, status );
+                  if( astOK ) newnode->arg[1]->mult = NULL;
+               }
+            }
+         }
+
+/* If no multiplier is supplied, the replacement node is simply the tree
+   which defines the unscaled unit (if known), or the original node (if
+   unknown). */
+      } else {
+         unit = (*node)->unit;
+         if( unit && unit->head ) newnode = CopyTree( unit->head, status );
+      }
+
+/* If this is not an LDVAR Node, remake the sub-trees which form the
+   arguments of this node. */
+   } else {
+      for( i = 0; i < (*node)->narg; i++ ) {
+         RemakeTree( &((*node)->arg[ i ]), status );
+      }
+   }
+
+/* If an error has occurred, free any new node. */
+   if( !astOK ) newnode = FreeTree( newnode, status );
+
+/* If we have a replacement node, free the supplied tree and return a
+   pointer to the new tree. */
+   if( newnode ) {
+      FreeTree( *node, status );
+      *node = newnode;
+   }
+}
+
+static int ReplaceNode( UnitNode *target, UnitNode *old, UnitNode *new, int *status ) {
+/*
+*  Name:
+*     ReplaceNode
+
+*  Purpose:
+*     Replace a node within a tree of UnitNodes.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     int ReplaceNode( UnitNode *target, UnitNode *old, UnitNode *new, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function replaces a specified node within a tree of UnitNodes
+*     with another given node. The original node is freed if found.
+
+*  Parameters:
+*     target
+*        A pointer to the UnitNode at the head of the tree containing the
+*        node to be replaced.
+*     old
+*        A pointer to the UnitNode to be replaced.
+*     new
+*        A pointer to the UnitNode to replace "old".
+*     status
+*        Pointer to the inherited status variable.
+
+*  Return Value:
+*     Non-zero if the "old" node was found and replaced (in which case
+*     the "old" node will have been freed).
+
+*  Notes:
+*     - It is assumed that the "old" node occurs at most once within the
+*     target tree.
+*     - The node at the head of the target tree is not compared with the
+*     "old" node. It is assumed the called will already have done this.
+*     - A value of zero is returned if an error has already occurred, or
+*     if this function fails for any reason.
+
+*/
+
+/* Local Variables: */
+   int i;
+   int result;
+
+/* Initialise */
+   result = 0;
+
+/* Check inherited status. */
+   if( !astOK ) return result;
+
+/* Loop round the arguments of the node at the head of the target tree.
+   Break out of the loop as soone as the old node is found. */
+   for( i = 0; i < target->narg; i++ ) {
+
+/* If this argument is the node to be replaced, free the old one and store
+   the new one, and then leave the loop. */
+      if( target->arg[ i ] == old ) {
+         FreeTree( old, status );
+         target->arg[ i ] = new;
+         result = 1;
+         break;
+
+/* Otherwise use this function recursively to search for the old node
+   within the current argument. */
+      } else {
+         if( ReplaceNode( target->arg[ i ], old, new, status ) ) break;
+      }
+   }
+
+/* If an error has occurred, return zero. */
+   if( !astOK ) result = 0;
+
+/* Return the answer. */
+   return result;
+}
+
+static int SimplifyTree( UnitNode **node, int std, int *status ) {
+/*
+*  Name:
+*     SimplifyTree
+
+*  Purpose:
+*     Simplify a tree of UnitNodes.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     int SimplifyTree( UnitNode **node, int std, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function simplifies a tree of UnitNodes. It is assumed that
+*     all the OP_LDVAR nodes in the tree refer to the same basic unit.
+*     A primary purpose of this function is to standardise the tree so
+*     that trees which implement equivalent transformations but which
+*     have different structures can be compared (for instance, so that
+*     "2*x" and "x*2" are treated as equal trees). If "std" is non-zero,
+*     reducing the complexity of the tree is only of secondary importance.
+*     This explains why some "simplifications" actually produced trees which
+*     are more complicated.
+
+*  Parameters:
+*     node
+*        The address of a pointer to the UnitNode at the head of the tree
+*        which is to be simplified. On exit the supplied tree is freed and a
+*        pointer to a new tree is placed at the given address.
+*     std
+*        If non-zero, perform standardisations. Otherwise only perform
+*        genuine simplifications.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     Non-zero if some change was made to the tree.
+
+*/
+
+/* Local Variables: */
+   int i;
+   UnitNode *newnode;
+   UnitNode *node1;
+   UnitNode *node2;
+   Oper op;
+   UnitNode **factors;
+   double *powers;
+   int nfactor;
+   double coeff;
+   int result;
+
+/* Initialise */
+   result = 0;
+
+/* Check inherited status. */
+   if( !astOK ) return result;
+
+/* Initiallially, we have no replacement node. */
+   newnode = NULL;
+
+/* First replace any complex constant expressions any corresponding
+   OP_LDCON nodes. */
+   FixConstants( node, 0, status );
+
+/* Simplify the sub-trees corresponding to the arguments of the node at
+   the head of the supplied tree. */
+   for( i = 0; i < (*node)->narg; i++ ) {
+      if( SimplifyTree( &( (*node)->arg[ i ] ), std, status ) ) result = 1;
+   }
+
+/* Now do specific simplifications appropriate to the nature of the node at
+   the head of the tree. */
+   op = (*node)->opcode;
+
+/* Natural log */
+/* =========== */
+/* We standardise argument powers into coefficients of the LN value. */
+   if( op == OP_LN ) {
+
+/* If the argument is a OP_EXP node, they cancel out. Return a copy of the
+   argument of OP_EXP node. */
+      if( (*node)->arg[ 0 ]->opcode == OP_EXP ) {
+         newnode = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
+
+/* If the argument is an OP_POW node, rearrange the nodes to represent
+   k*ln(x) instead of ln(x**k) (note pow nodes always have a constant
+   exponent - this is checked in InvertConstants). SQRT arguments will
+   not occur because they will have been changed into POW nodes when the
+   arguments of the supplied head node were simplified above. */
+      } else if( std && (*node)->arg[ 0 ]->opcode == OP_POW ) {
+         newnode = NewNode( NULL, OP_MULT, status );
+         node1 = CopyTree( (*node)->arg[ 0 ]->arg[ 1 ], status );
+         node2 = NewNode( NULL, OP_LN, status );
+         if( astOK ) {
+            node2->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
+            newnode->arg[ 0 ] = node1;
+            newnode->arg[ 1 ] = node2;
+         }
+      }
+
+/* Common log */
+/* ========== */
+/* We standardise natural logs into common logs. */
+   } else if( op == OP_LOG ) {
+      if( std ) {
+         newnode = NewNode( NULL, OP_DIV, status );
+         node1 = NewNode( NULL, OP_LN, status );
+         node2 = NewNode( NULL, OP_LDCON, status );
+         if( astOK ) {
+            node1->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ], status );
+            node2->con = log( 10.0 );
+            newnode->arg[ 0 ] = node1;
+            newnode->arg[ 1 ] = node2;
+         }
+      }
+
+/* Exponential */
+/* =========== */
+/* We prefer to minimise the number of EXP nodes, so, for instance, we do not
+   change "exp(x*y)" to "exp(x)+exp(y)" (and the code for ADD nodes does
+   the inverse conversion). */
+   } else if( op == OP_EXP ) {
+
+/* If the argument is an OP_LN node, they cancel out. Return a copy of the
+   argument of the OP_LN node. Common log arguments will not occur because
+   they will have been changed into natural logs when the arguments of
+   the supplied head node were simplified above. */
+      if( (*node)->arg[ 0 ]->opcode == OP_LN ) {
+         newnode = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
+      }
+
+/* Square root */
+/* =========== */
+/* We standardise sqrt nodes into pow nodes. */
+   } else if( op == OP_SQRT ) {
+      if( std ) {
+         newnode = NewNode( NULL, OP_POW, status );
+         node1 = CopyTree( (*node)->arg[ 0 ], status );
+         node2 = NewNode( NULL, OP_LDCON, status );
+         if( astOK ) {
+            node2->con = 0.5;
+            newnode->arg[ 0 ] = node1;
+            newnode->arg[ 1 ] = node2;
+         }
+      }
+
+/* Exponentiation */
+/* ============== */
+/* We want to simplfy factors. So, for instance, (x*y)**k is converted to
+   (x**k)*(y**k). */
+   } else if( op == OP_POW ) {
+
+/* If the first argument is an OP_EXP node, then change "(e**x)**k" into
+   "e**(k*x)" */
+      if( (*node)->arg[ 0 ]->opcode == OP_EXP ) {
+         newnode = NewNode( NULL, OP_EXP, status );
+         node1 = NewNode( NULL, OP_MULT, status );
+         if( astOK ) {
+            node1->arg[ 0 ] = CopyTree( (*node)->arg[ 1 ], status );
+            node1->arg[ 1 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
+            newnode->arg[ 0 ] = node1;
+         }
+
+/* "x**0" can be replaced by 1.0 */
+      } else if( (*node)->arg[ 1 ]->con == 0.0 ) {
+         newnode = NewNode( NULL, OP_LDCON, status );
+         if( astOK ) newnode->con = 1.0;
+
+/* "x**1" can be replaced by x */
+      } else if( astEQUAL( (*node)->arg[ 1 ]->con, 1.0 ) ) {
+         newnode = CopyTree( (*node)->arg[ 0 ], status );
+
+/* If the first argument is an OP_POW node, then change "(x**k1)**k2" into
+   "x**(k1*k2)" */
+      } else if( (*node)->arg[ 0 ]->opcode == OP_POW ) {
+         newnode = NewNode( NULL, OP_POW, status );
+         node1 = NewNode( NULL, OP_LDCON, status );
+         if( astOK ) {
+            node1->con = ( (*node)->arg[ 0 ]->arg[ 1 ]->con )*
+                         ( (*node)->arg[ 1 ]->con );
+            newnode->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
+            newnode->arg[ 1 ] = node1;
+         }
+
+/* If the first argument is an OP_MULT node, then change "(x*y)**k" into
+   "(x**(k))*(y**(k))" */
+      } else if( std && (*node)->arg[ 0 ]->opcode == OP_MULT ) {
+         newnode = NewNode( NULL, OP_MULT, status );
+         node1 = NewNode( NULL, OP_POW, status );
+         if( astOK ) {
+            node1->arg[ 1 ] = CopyTree( (*node)->arg[ 1 ], status );
+            node2 = CopyTree( node1, status );
+            node1->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
+            node2->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 1 ], status );
+            newnode->arg[ 0 ] = node1;
+            newnode->arg[ 1 ] = node2;
+         }
+      }
+
+/* Division. */
+/* ========= */
+/* We standardise divisions into corresponding multiplications. */
+   } else if( op == OP_DIV ) {
+
+/* Division by 1 is removed. */
+      if( astEQUAL( (*node)->arg[ 1 ]->con, 1.0 ) ){
+         newnode = CopyTree( (*node)->arg[ 0 ], status );
+
+/* Division by any other constant (except zero) is turned into a
+   multiplication by the reciprocal constant. */
+      } else if( (*node)->arg[ 1 ]->con != AST__BAD ) {
+         if( (*node)->arg[ 1 ]->con != 0.0 ) {
+            newnode = NewNode( NULL, OP_MULT, status );
+            node1 = NewNode( NULL, OP_LDCON, status );
+            if( astOK ) {
+               node1->con = 1.0/(*node)->arg[ 1 ]->con;
+               newnode->arg[ 0 ] = node1;
+               newnode->arg[ 1 ] = CopyTree( (*node)->arg[ 0 ], status );
+            }
+         } else {
+            astError( AST__BADUN, "Simplifying a units expression"
+                      "requires a division by zero." , status);
+         }
+
+/* Other divisions "x/y" are turned into "x*(y**(-1))" */
+      } else if( std ) {
+         newnode = NewNode( NULL, OP_MULT, status );
+         node1 = NewNode( NULL, OP_POW, status );
+         node2 = NewNode( NULL, OP_LDCON, status );
+         if( astOK ) {
+            node2->con = -1.0;
+            node1->arg[ 0 ] = CopyTree( (*node)->arg[ 1 ], status );
+            node1->arg[ 1 ] = node2;
+            newnode->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ], status );
+            newnode->arg[ 1 ] = node1;
+         }
+      }
+
+/* Multiplication */
+/* ============== */
+   } else if( op == OP_MULT ) {
+
+/* If the right hand argument is constant, swap the arguments. */
+      if( (*node)->arg[ 1 ]->con != AST__BAD ) {
+         newnode = NewNode( NULL, OP_MULT, status );
+         if( astOK ) {
+            newnode->arg[ 0 ] = CopyTree( (*node)->arg[ 1 ], status );
+            newnode->arg[ 1 ] = CopyTree( (*node)->arg[ 0 ], status );
+         }
+
+/* Multiplication by zero produces a constant zero. */
+      } else if( (*node)->arg[ 0 ]->con == 0.0 ){
+         newnode = NewNode( NULL, OP_LDCON, status );
+         if( astOK ) newnode->con = 0.0;
+
+/* Multiplication by 1 is removed. */
+      } else if( astEQUAL( (*node)->arg[ 0 ]->con, 1.0 ) ){
+         newnode = CopyTree( (*node)->arg[ 1 ], status );
+
+/* For other MULT nodes, analyse the tree to find a list of all its
+   factors with an associated power for each one, and an overall constant
+   coefficient. */
+      } else if( std ) {
+         FindFactors( (*node), &factors, &powers, &nfactor, &coeff, status );
+
+/* Produce a new tree from these factors. The factors are standardised by
+   ordering them alphabetically (after conversion to a character string). */
+         newnode = CombineFactors( factors, powers, nfactor, coeff, status );
+
+/* Free resources */
+         factors = astFree( factors );
+         powers = astFree( powers );
+
+      }
+   }
+
+/* If we have produced a new node which is identical to the old node,
+   free it. Otherwise, indicate we have made some changes. */
+   if( newnode ) {
+      if( !CmpTree( newnode, *node, 1, status ) ) {
+         newnode = FreeTree( newnode, status );
+      } else {
+         result = 1;
+      }
+   }
+
+/* If an error has occurred, free any new node. */
+   if( !astOK ) newnode = FreeTree( newnode, status );
+
+/* If we have a replacement node, free the supplied tree and return a
+   pointer to the new tree. */
+   if( newnode ) {
+      FreeTree( *node, status );
+      *node = newnode;
+   }
+
+/* If the above produced some change, try re-simplifying the tree. */
+   if( result ) SimplifyTree( node, std, status );
+
+/* Return the result. */
+   return result;
+
+}
+
+static int SplitUnit( const char *str, int ls, const char *u, int cs,
+                      Multiplier **mult, int *l, int *status ) {
+/*
+*  Name:
+*     SplitUnit
+
+*  Purpose:
+*     Split a given string into unit name and multiplier.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     int SplitUnit( const char *str, int ls, const char *u, int cs,
+*                    Multiplier **mult, int *l, int *status  )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     Returns non-zer0 if the supplied string ends with the supplied unit
+*     name or label, and any leading string is a known multiplier.
+
+*  Parameters:
+*     str
+*        The string to test, typically containing a multiplier and a unit
+*        symbol or label.
+*     ls
+*        Number of characters to use from "str" (not including trailing null)
+*     u
+*        Pointer to the unit label or symbol string to be searched for.
+*     cs
+*        If non-zero, the test for "u" is case insensitive.
+*     mult
+*        Address of a location at which to return the multiplier at the
+*        start of the supplied string. NULL is returned if the supplied
+*        string does not match the supplied unit, or if the string
+*        includes no multiplier.
+*     l
+*        Address of an int in which to return the length of "u".
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     Non-zero if "str" ends with "u" and starts with a null string or a
+*     known multiplier string.
+
+*/
+
+/* Local Variables: */
+   int ret;
+   int lm;
+   int lu;
+
+/* Initialise */
+   ret = 0;
+   *mult = NULL;
+   *l = 0;
+
+/* Check inherited status. */
+   if( !astOK ) return ret;
+
+/* Find the number of characters in the supplied unit label or symbol. The
+   difference between lu and ls must be the length of the multiplier. */
+   lu = strlen( u );
+   lm = ls - lu;
+
+/* Make sure "str" is not shorter than "u" */
+   if( lm >= 0 ) {
+
+/* Compare the end of "str" against "u" */
+      if( cs ? !strncmp( str + lm, u, lu ) :
+               !Ustrncmp( str + lm, u, lu, status ) ) {
+         ret = 1;
+
+/* If "str" ends with "u", see if it starts with a known multiplier */
+         if( lm > 0 ) {
+            ret = 0;
+            *mult = GetMultipliers( status );
+            while( *mult ) {
+               if( (*mult)->symlen == lm && !strncmp( str, (*mult)->sym, lm ) ) {
+                  ret = 1;
+                  break;
+               }
+               *mult = (*mult)->next;
+            }
+
+/* If not, try again using case-insensitive matching. */
+            if( !ret ) {
+               *mult = GetMultipliers( status );
+               while( *mult ) {
+                  if( (*mult)->symlen == lm && !Ustrncmp( str, (*mult)->sym, lm, status ) ) {
+                     ret = 1;
+                     break;
+                  }
+                  *mult = (*mult)->next;
+               }
+            }
+
+/* If not, try again using case-insensitive matching against the
+   multiplier label. */
+            if( !ret ) {
+               *mult = GetMultipliers( status );
+               while( *mult ) {
+                  if( (*mult)->lablen == lm && !Ustrncmp( str, (*mult)->label, lm, status ) ) {
+                     ret = 1;
+                     break;
+                  }
+                  *mult = (*mult)->next;
+               }
+            }
+
+         }
+      }
+   }
+
+   *l = lu;
+   return ret;
+}
+
+double astUnitAnalyser_( const char *in, double powers[9], int *status ){
+/*
+*+
+*  Name:
+*     astUnitAnalyser
+
+*  Purpose:
+*     Perform a dimensional analysis of a unti string.
+
+*  Type:
+*     Protected function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     double astUnitAnalyser_( const char *in, double powers[9] )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function parses the supplied units string if possible, and
+*     returns a set of pwoers and a scaling factor which represent the
+*     units string.
+
+*  Parameters:
+*     in
+*        A string representation of the units, for instance "km/h".
+*     powers
+*        An array in which are returned the powers for each of the following
+*        basic units (in the order shown): kilogramme, metre, second, radian,
+*        Kelvin, count, adu, photon, magnitude, pixel. If the supplied unit
+*        does not depend on a given basic unit a value of 0.0 will be returned
+*        in the array. The returns values represent a system of units which is
+*        a scaled form of the supplied units, expressed in the basic units of
+*        m, kg, s, rad, K, count, adu, photon, mag and pixel. For instance, a
+*        returned array of [1,0,-2,0,0,0,0,0,0] would represent "m/s**2".
+
+*  Returned Value:
+*     A scaling factor for the supplied units. The is the value, in the
+*     units represented by the returned powers, which corresponds to a
+*     value of 1.0 in the supplied units.
+
+*  Notes:
+*     -  An error will be reported if the units string cannot be parsed
+*     or refers to units or functions which cannot be analysed in this way.
+*     -  AST__BAD is returned if this function is invoked with the
+*     global error status set or if it should fail for any reason.
+*-
+*/
+
+/* Local Variables: */
+   UnitNode *in_tree;
+   double result;
+
+/* Initialise */
+   result = AST__BAD;
+
+/* Check the global error status. */
+   if ( !astOK ) return result;
+
+/* Parse the input units string, producing a tree of UnitNodes which
+   represents the input units. A pointer to the UnitNode at the head of
+   the tree is returned if succesfull. Report a context message if this
+   fails. */
+   in_tree = CreateTree( in, 1, 1, status );
+   if( in_tree ) {
+
+/* Analyse the tree */
+      if( !DimAnal( in_tree, powers, &result, status ) && astOK ) {
+         result = AST__BAD;
+         astError( AST__BADUN, "astUnitAnalyser: Error analysing input "
+                   "units string '%s' (it may contain unsupported "
+                   "functions or dimensionless units).", status, in );
+      }
+
+/* Free the tree. */
+      in_tree = FreeTree( in_tree, status );
+
+   } else if( astOK ) {
+      astError( AST__BADUN, "astUnitAnalyser: Error parsing input "
+                "units string '%s'.", status, in );
+   }
+
+/* Return the result */
+   return result;
+}
+
+const char *astUnitLabel_( const char *sym, int *status ){
+/*
+*+
+*  Name:
+*     astUnitLabel
+
+*  Purpose:
+*     Return a string label for a given unit symbol.
+
+*  Type:
+*     Protected function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     const char *astUnitLabel( const char *sym )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function returns a pointer to a constant string containing a
+*     descriptive label for the unit specified by the given unit symbol.
+
+*  Parameters:
+*     sym
+*        A string holing a known unit symbol.
+
+*  Returned Value:
+*     A pointer to constant string holding a descriptive label for the
+*     supplied unit. A NULL pointer is returned (without error) if the
+*     supplied unit is unknown.
+
+*  Notes:
+*     -  A NULL pointer is returned if this function is invoked with the
+*     global error status set or if it should fail for any reason.
+*-
+*/
+
+/* Local Variables: */
+   const char *result;
+   KnownUnit *unit;
+
+/* Initialise */
+   result = NULL;
+
+/* Check the global error status. */
+   if ( !astOK ) return result;
+
+/* Ensure descriptions of the known units are available. */
+   unit = GetKnownUnits( 1, status );
+
+/* Loop through the chain of known units looking for a unit with a symbol
+   equal to the supplied string. If found, store a pointer to its label
+   and break out of the loop. */
+   while( unit ) {
+      if( !strcmp( sym, unit->sym ) ) {
+         result = unit->label;
+         break;
+      }
+
+      unit = unit->next;
+   }
+
+/* Return the answer. */
+   return result;
+}
+
+AstMapping *astUnitMapper_( const char *in, const char *out,
+                            const char *in_lab, char **out_lab, int *status ){
+/*
+*+
+*  Name:
+*     astUnitMapper
+
+*  Purpose:
+*     Create a Mapping between two system of units.
+
+*  Type:
+*     Protected function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     AstMapping *astUnitMapper( const char *in, const char *out,
+*                                const char *in_lab, char **out_lab )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function creates a Mapping between two specified system of
+*     units. It also modifes a supplied label (which is typically
+*     the axis label associated with the input units) so that it includes
+*     any functional change implied by the supplied "in" and "out" units.
+
+*  Parameters:
+*     in
+*        A string representation of the input units, for instance "km/h".
+*        See "Unit Representations:" below.
+*     out
+*        A string representation of the output units, for instance "m/s".
+*        See "Unit Representations:" below.
+*     in_lab
+*        A label describing the quantity associated with the input units.
+*        If the "in" string is the Units attribute of an Axis, then
+*        "in_lab" should be the Label of the same Axis. May be supplied
+*        NULL in which case "out_lab" is ignored.
+*     out_lab
+*        The address at which to return a pointer to a label describing the
+*        quantity associated with the output units. For instance, if the
+*        input and output units are "Hz" and "sqrt(Hz)", and the input
+*        label is "Frequency", then the returned output label will be
+*        "sqrt( Frequency )". The returned label is stored in dynamically
+*        allocated memory which should be freed (using astFree) when no longer
+*        needed.
+
+*  Returned Value:
+*     A pointer to a Mapping which can be used to transform values in the
+*     "in" system of units into the "out" system of units. The Mapping
+*     will have 1 input and 1 output.
+
+*  Unit Representations:
+*     The string supplied for "in" and "out" should represent a system of
+*     units following the recommendations of the FITS WCS paper I
+*     "Representation of World Coordinates in FITS" (Greisen & Calabretta).
+*     Various commonly used variants are also allowed.
+*
+*     To summarise, a string describing a system of units should be an
+*     algebraic expression which combines one or more named units. The
+*     following functions and operators may be used within these algebraic
+*     expressions:
+*
+*     - "*": multiplication. A period "." or space " " may also be used
+*       to represent multiplication (a period is only interpreted as a
+*       multiplication operator if it is not positioned between two digits,
+*       and a space is only interpreted as a multiplication operator if it
+*       occurs between two operands).
+*     - "/": division.
+*     - "**": exponentiation. The exponent (i.e. the operand following the
+*       exponentiation operator) must be a constant. The symbol "^" is also
+*       interpreted as an exponentiation operator. Exponentiation is also
+*       implied by an integer following a unit name without any separator
+*       (e.g. "cm2" is "cm^2").
+*     - log(): Common logarithm.
+*     - ln(): Natural logarithm.
+*     - sqrt(): Square root.
+*     - exp(): Exponential.
+*
+*     Function names are case insensitive. White space may be included
+*     within an expression (note that white space between two operands
+*     will be interpreted as a muiltiplication operator as described
+*     above). Parentheses may be used to indicate the order in which
+*     expressions are to be evaluated (normal mathematical precedence is
+*     used otherwise). The following symbols may be used to represent
+*     constants:
+*
+*     - "pi"
+*     - "e"
+*
+*     These symbols are also case in-sensitive.
+*
+*     The above operators and functions are used to combine together one
+*     or more "unit symbols". The following base unit symbols are recognised:
+*
+*     - "m":  metre.
+*     - "g":  gram.
+*     - "s":  second.
+*     - "rad":  radian.
+*     - "sr":  steradian.
+*     - "K":  Kelvin.
+*     - "mol":  mole.
+*     - "cd":  candela.
+*
+*     The following symbols for units derived fro the above basic units are
+*     recognised:
+*
+*     - "sec":  second (1 s)
+*     - "Hz":  Hertz  (1/s).
+*     - "N":  Newton  (kg m/s**2).
+*     - "J":  Joule  (N m).
+*     - "W":  Watt  (J/s).
+*     - "C":  Coulomb  (A s).
+*     - "V":  Volt  (J/C).
+*     - "Pa":  Pascal  (N/m**2).
+*     - "Ohm":  Ohm  (V/A).
+*     - "S":  Siemens  (A/V).
+*     - "F":  Farad  (C/V).
+*     - "Wb":  Weber  (V s).
+*     - "T":  Tesla  (Wb/m**2).
+*     - "H":  Henry  (Wb/A).
+*     - "lm":  lumen  (cd sr).
+*     - "lx":  lux  (lm/m**2).
+*     - "deg":  degree  (pi/180 rad).
+*     - "arcmin":  arc-minute  (1/60 deg).
+*     - "arcsec":  arc-second  (1/3600 deg).
+*     - "mas":  milli-arcsecond  (1/3600000 deg).
+*     - "min":  minute  (60 s).
+*     - "h":  hour  (3600 s).
+*     - "d":  day  (86400 s).
+*     - "yr":  year  (31557600 s).
+*     - "a":  year  (31557600 s).
+*     - "eV":  electron-Volt  (1.60217733E-19 J).
+*     - "erg":  erg  (1.0E-7 J).
+*     - "Ry":  Rydberg  (13.605692 eV).
+*     - "solMass":  solar mass  (1.9891E30 kg).
+*     - "u":  unified atomic mass unit  (1.6605387E-27 kg).
+*     - "solLum":  solar luminosity  (3.8268E26 W).
+*     - "Angstrom":  Angstrom  (1.0E-10 m).
+*     - "Ang":  Angstrom
+*     - "A":  Ampere
+*     - "micron":  micron (1.0E-6 m).
+*     - "solRad":  solar radius  (6.9599E8 m).
+*     - "AU":  astronomical unit  (1.49598E11 m).
+*     - "lyr":  light year  (9.460730E15 m).
+*     - "pc":  parsec  (3.0867E16 m).
+*     - "count":  count.
+*     - "ct":  count.
+*     - "adu":  analogue-to-digital converter unit.
+*     - "photon":  photon.
+*     - "ph":  photon.
+*     - "Jy":  Jansky  (1.0E-26 W /m**2 /Hz).
+*     - "Jan":  Jansky
+*     - "mag":  magnitude.
+*     - "G":  Gauss  (1.0E-4 T).
+*     - "pixel":  pixel.
+*     - "pix":  pixel.
+*     - "barn":  barn  (1.0E-28 m**2).
+*     - "D":  Debye  (1.0E-29/3 C.m).
+*
+*     In addition, any other unknown unit symbol may be used (but of course
+*     no mapping will be possible between unknown units).
+*
+*     Unit symbols may be preceded with a numerical constant (for
+*     instance "1000 m") or a standard multiplier symbol (for instance "km")
+*     to represent some multiple of the unit. The following standard
+*     multipliers are recognised:
+*
+*     - "d":   deci   (1.0E-1)
+*     - "c":   centi  (1.0E-2)
+*     - "m":   milli  (1.0E-3)
+*     - "u":   micro  (1.0E-6)
+*     - "n":   nano   (1.0E-9)
+*     - "p":   pico   (1.0E-12)
+*     - "f":   femto  (1.0E-15)
+*     - "a":   atto   (1.0E-18)
+*     - "z":   zepto  (1.0E-21)
+*     - "y":   yocto  (1.0E-24)
+*     - "da":  deca   (1.0E1)
+*     - "h":   hecto  (1.0E2)
+*     - "k":   kilo   (1.0E3)
+*     - "M":   mega   (1.0E6)
+*     - "G":   giga   (1.0E9)
+*     - "T":   tera   (1.0E12)
+*     - "P":   peta   (1.0E15)
+*     - "E":   exa    (1.0E18)
+*     - "Z":   zetta  (1.0E21)
+*     - "Y":   yotta  (1.0E24)
+
+*  Notes:
+*     -  NULL values are returned without error if the supplied units are
+*     incompatible (for instance, if the input and output units are "kg"
+*     and "m" ).
+*     -  NULL values are returned if this function is invoked with the
+*     global error status set or if it should fail for any reason.
+*-
+*/
+
+/* Local Variables: */
+   AstMapping *result;
+   UnitNode **units;
+   UnitNode *in_tree;
+   UnitNode *intemp;
+   UnitNode *inv;
+   UnitNode *labtree;
+   UnitNode *newtest;
+   UnitNode *out_tree;
+   UnitNode *outtemp;
+   UnitNode *src;
+   UnitNode *testtree;
+   UnitNode *tmp;
+   UnitNode *totaltree;
+   UnitNode *totlabtree;
+   const char *c;
+   const char *exp;
+   int i;
+   int nc;
+   int nunits;
+   int ipass;
+
+/* Initialise */
+   result = NULL;
+   if( in_lab ) *out_lab = NULL;
+
+/* Check the global error status. */
+   if ( !astOK ) return result;
+
+/* A quick check for a common simple case: if the two strings are
+   identical, return a UnitMap.*/
+   if( !strcmp( in, out ) ) {
+      if( in_lab ) *out_lab = astStore( NULL, in_lab, strlen( in_lab ) + 1 );
+      return (AstMapping *) astUnitMap( 1, "", status );
+   }
+
+/* More initialisation. */
+   in_tree = NULL;
+   out_tree = NULL;
+   units = NULL;
+
+/* Parse the input units string, producing a tree of UnitNodes which
+   represents the input units. A pointer to the UnitNode at the head of
+   the tree is returned if succesfull. Report a context message if this
+   fails. The returned tree contains branch nodes which correspond to
+   operators or functions, and leaf nodes which represent constant values
+   or named basic units (m, s, g, K, etc). Each branch node has one or more
+   "arguments" (i.e. child nodes) which are operated on or combined by
+   the branch node in some way to produce the nodes "value". This value
+   is then used as an argument for the node's parent node (if any). If
+   the string supplied by the user refers to any known derived units (e.g. "N",
+   Newton) then each such unit is represented in the returned tree by a
+   complete sub-tree in which the head node corresponds to the derived
+   unit (e.g. "N") and the leaf nodes correspond to the basic units needed
+   to define the derived unit ( for instance, "m", "s" and "g" - metres,
+   seconds and grammes), or numerical constants. Thus every leaf node in the
+   returned tree will be a basic unit (i.e. a unit which is not defined in
+   terms of other units), or a numerical constant. */
+   in_tree = CreateTree( in, 1, 1, status );
+   if( !astOK ) astError( AST__BADUN, "astUnitMapper: Error parsing input "
+                          "units string '%s'.", status, in );
+
+/* Do the same for the output units. */
+   if( astOK ) {
+      out_tree = CreateTree( out, 1, 1, status );
+      if( !astOK ) astError( AST__BADUN, "astUnitMapper: Error parsing output "
+                             "units string '%s'.", status, out );
+   }
+
+/* If a blank string is supplied for both input and output units, then
+   assume a UnitMap is the appropriate Mapping. */
+   if( !in_tree && !out_tree && astOK ) {
+      result = (AstMapping *) astUnitMap( 1, "", status );
+      if( in_lab ) *out_lab = astStore( NULL, in_lab, strlen( in_lab ) + 1 );
+
+/* Otherwise, if we have both input and output trees... */
+   } else if( in_tree && out_tree && astOK ) {
+
+/* Locate all the basic units used within either of these two trees. An
+   array is formed in which each element is a pointer to a UnitNode
+   contained within one of the trees created above. Each basic unit
+   referred to in either tree will have a single entry in this array
+   (even if the unit is referred to more than once). */
+      units = NULL;
+      nunits = 0;
+      LocateUnits( in_tree, &units, &nunits, status );
+      LocateUnits( out_tree, &units, &nunits, status );
+
+/* Due to the simple nature of the simplification process in SimplifyTree,
+   the following alogorithm sometimes fails to find a Mapping form input
+   to output units, but can find a Mapping from output to input units.
+   In this latter case, we can get the required Mapping from input to
+   output  simply by inverting the Mapign from output to input. So try
+   first with the units in the original order. If this fails to find a
+   Mapping, try again with the units swapped, and note that the final
+   Mapping should be inverted before being used. */
+      for( ipass = 0; ipass < 2; ipass++ ){
+         if( ipass == 1 ) {
+            tmp = in_tree;
+            in_tree = out_tree;
+            out_tree = tmp;
+         }
+
+/* We are going to create a new tree of UnitNodes in which the head node
+   corresponds to the requested output units, and which has a single
+   non-constant leaf node corresponding to the input units. Initialise a
+   pointer to this new tree to indicate that it has not yet been created. */
+         testtree = NULL;
+
+/* Loop round each basic unit used in the definition of either the input
+   or the output units (i.e. the elements of the array created above by
+   "LocateUnits"). The unit selected by this loop is referred to as the
+   "current" unit. On each pass through this loop, we create a tree which
+   is a candidate for the final required tree (the "test tree" pointed to
+   by the testtree pointer initialised above). In order for a mapping to
+   be possible between input and output units, the test tree created on
+   each pass through this loop must be equivalent to the test tree for the
+   previous pass (in other words, all the test trees must be equivalent).
+   We break out of the loop (and return a NULL Mapping) as soon as we find
+   a test tree which differs from the previous test tree. */
+         for( i = 0; i < nunits; i++ ) {
+
+/* Create copies of the trees describing the input and output units, in which
+   all units other than the current unit are set to a constant value of 1.
+   This is done by replacing OP_LDVAR nodes (i.e. nodes which "load" the
+   value of a named basic unit) by OP_LDCON nodes (i.e. nodes which load
+   a specified constant value) in the tree copy. */
+            intemp = FixUnits( in_tree, units[ i ], status );
+            outtemp = FixUnits( out_tree, units[ i ], status );
+
+/* Simplify these trees. An important side-effect of this simplification
+   is that trees are "standardised" which allows them to be compared for
+   equivalence. A single mathematical expression can often be represented
+   in many different ways (for instance "A/B" is equivalent to "(B**(-1))*A").
+   Standardisation is a process of forcing all equivalent representations
+   into a single "standard" form. Without standardisation, trees representing
+   the above two expressions would not be considered to be equivalent
+   since thy would contain different nodes and have different structures.
+   As a consequence of this standardisation, the "simplification" performed
+   by SimplifyTree can sometimes actually make the tree more complicated
+   (in terms of the number of nodes in the tree). */
+            SimplifyTree( &intemp, 1, status );
+            SimplifyTree( &outtemp, 1, status );
+
+/* If either of the simplified trees does not depend on the current unit,
+   then the node at the head of the simplified tree will have a constant
+   value (because all the units other than the current unit have been fixed
+   to a constant value of 1.0 above by FixUnits, leaving only the current
+   unit to vary in value). If both simplified trees are constants, then
+   neither tree depends on the current basic unit (i.e. references to the
+   current basic unit cancel out within each string expression - for
+   instance if converting from "m.s.Hz" to "km" and the current unit
+   is "s", then the "s.Hz" term will cause the "s" units to cancel out). In
+   this case ignore this basic unit and pass on to the next. */
+            if( outtemp->con != AST__BAD && intemp->con != AST__BAD ) {
+
+/* If just one simplified tree is constant, then the two units cannot
+   match since one depends on the current basic unit and the other does
+   not. Free any test tree from previous passes and break out of the loop. */
+            } else if( outtemp->con != AST__BAD || intemp->con != AST__BAD ) {
+               intemp = FreeTree( intemp, status );
+               outtemp = FreeTree( outtemp, status );
+               testtree = FreeTree( testtree, status );
+               break;
+
+/* If neither simplified tree is constant, both depend on the current
+   basic unit and so we can continue to see if their dependencies are
+   equivalent. */
+            } else {
+
+/* We are going to create a new tree which is the inverse of the above
+   simplified "intemp" tree. That is, the new tree will have a head node
+   corresponding to the current unit, and a single non-constant leaf node
+   corresponding to the input units. Create an OP_LDVAR node which can be
+   used as the leaf node for this inverted tree. If the input tree is
+   inverted successfully, this root node becomes part of the inverted tree,
+   and so does not need to be freed explicitly (it will be freed when the
+   inverted tree is freed). */
+               src = NewNode( NULL, OP_LDVAR, status );
+               if( astOK ) src->name = astStore( NULL, "input_units", 12 );
+
+/* Now produce the inverted input tree. If the tree cannot be inverted, a
+   null pointer is returned. Check for this. Otherwise a pointer to the
+   UnitNode at the head of the inverted tree is returned. */
+               inv = InvertTree( intemp, src, status );
+               if( inv ) {
+
+/* Concatenate this tree (which goes from "input units" to "current unit")
+   with the simplified output tree (which goes from "current unit" to
+   "output units"), to get a new tree which goes from input units to output
+   units. */
+                  totaltree = ConcatTree( inv, outtemp, status );
+
+/* Simplify this tree. */
+                  SimplifyTree( &totaltree, 1, status );
+
+/* Compare this simplified tree with the tree produced for the previous
+   unit (if any). If they differ, we cannot map between the supplied
+   units so annul the test tree and break out of the loop. If this is the
+   first unit to be tested, use the total tree as the test tree for the
+   next unit. */
+                  if( testtree ) {
+                     if( CmpTree( totaltree, testtree, 0, status ) ) testtree = FreeTree( testtree, status );
+                     totaltree = FreeTree( totaltree, status );
+                     if( !testtree ) break;
+                  } else {
+                     testtree = totaltree;
+                  }
+               }
+
+/* If the input tree was inverted, free the inverted tree. */
+               if( inv ) {
+                  inv = FreeTree( inv, status );
+
+/* If the input tree could not be inverted, we cannot convert between input
+   and output units. Free the node which was created to be the root of the
+   inverted tree (and which has consequently not been incorporated into the
+   inverted tree), free any testtree and break out of the loop. */
+               } else {
+                  src = FreeTree( src, status );
+                  testtree = FreeTree( testtree, status );
+                  break;
+               }
+            }
+
+/* Free the other trees. */
+            intemp = FreeTree( intemp, status );
+            outtemp = FreeTree( outtemp, status );
+
+         }
+
+/* If all the basic units used by either of the supplied system of units
+   produced the same test tree, leave the "swap in and out units" loop. */
+         if( testtree ) break;
+
+      }
+
+/* If the input and output units have been swapped, swap them back to
+   their original order, and invert the test tree (if there is one). */
+      if( ipass > 0 ) {
+         tmp = in_tree;
+         in_tree = out_tree;
+         out_tree = tmp;
+         if( testtree ) {
+            src = NewNode( NULL, OP_LDVAR, status );
+            if( astOK ) src->name = astStore( NULL, "input_units", 12 );
+            newtest = InvertTree( testtree, src, status );
+            FreeTree( testtree, status );
+            testtree = newtest;
+            if( !newtest ) src = FreeTree( src, status );
+         }
+      }
+
+/* If all the basic units used by either of the supplied system of units
+   produced the same test tree, create a Mapping which is equivalent to the
+   test tree and return it. */
+      if( testtree ) {
+         result = MakeMapping( testtree, status );
+
+/* We now go on to produce the output axis label from the supplied input
+   axis label. Get a tree of UnitNodes which describes the supplied label
+   associated with the input axis. The tree will have single OP_LDVAR node
+   corresponding to the basic label (i.e. the label without any single
+   argument functions or exponentiation operators applied). */
+         if( in_lab && astOK ) {
+
+/* Get a pointer to the first non-blank character, and store the number of
+   characters to examine (this excludes any trailing white space). */
+            exp = in_lab;
+            while( isspace( *exp ) ) exp++;
+            c = exp + strlen( exp ) - 1;
+            while( c >= exp && isspace( *c ) ) c--;
+            nc = c - exp + 1;
+
+/* Create the tree. */
+            labtree = MakeLabelTree( exp, nc, status );
+            if( astOK ) {
+
+/* Concatenate this tree (which goes from "basic label" to "input label")
+   with the test tree found above (which goes from "input units" to "output
+   units"), to get a tree which goes from basic label to output label. */
+               totlabtree = ConcatTree( labtree, testtree, status );
+
+/* Simplify this tree. */
+               SimplifyTree( &totlabtree, 1, status );
+
+/* Create the output label from this tree. */
+               *out_lab = (char *) MakeExp( totlabtree, 0, 1, status );
+
+/* Free the trees. */
+               totlabtree = FreeTree( totlabtree, status );
+               labtree = FreeTree( labtree, status );
+
+/* Report a context error if the input label could not be parsed. */
+            } else {
+               astError( AST__BADUN, "astUnitMapper: Error parsing axis "
+                         "label '%s'.", status, in_lab );
+            }
+         }
+
+/* Free the units tree. */
+         testtree = FreeTree( testtree, status );
+
+      }
+   }
+
+/* Free resources. */
+   in_tree = FreeTree( in_tree, status );
+   out_tree = FreeTree( out_tree, status );
+   units = astFree( units );
+
+/* If an error has occurred, annul the returned Mapping. */
+   if( !astOK ) {
+      result = astAnnul( result );
+      if( in_lab ) *out_lab = astFree( *out_lab );
+   }
+
+/* Return the result. */
+   return result;
+
+}
+
+const char *astUnitNormaliser_( const char *in, int *status ){
+/*
+*+
+*  Name:
+*     astUnitNormalizer
+
+*  Purpose:
+*     Normalise a unit string into FITS-WCS format.
+
+*  Type:
+*     Protected function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     const char *astUnitNormaliser( const char *in )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     This function returns a standard FITS-WCS form of the supplied unit
+*     string.
+
+*  Parameters:
+*     in
+*        A string representation of the units, for instance "km/h".
+
+*  Returned Value:
+*     A pointer to a dynamically allocated string holding the normalized
+*     unit string. It should be freed using astFree when no longer needed.
+
+*  Notes:
+*     -  An error will be reported if the units string cannot be parsed.
+*     -  NULL is returned if this function is invoked with the
+*     global error status set or if it should fail for any reason.
+*-
+*/
+
+/* Local Variables: */
+   UnitNode *in_tree;
+   double dval;
+   const char *result;
+
+/* Initialise */
+   result = NULL;
+
+/* Check the global error status. */
+   if ( !astOK ) return result;
+
+/* Parse the input units string, producing a tree of UnitNodes which
+   represents the input units. A pointer to the UnitNode at the head of
+   the tree is returned if succesfull. Report a context message if this
+   fails. */
+   in_tree = CreateTree( in, 0, 1, status );
+   if( in_tree ) {
+
+/* Simplify the units expression, only doing genuine simplifications. */
+      SimplifyTree( &in_tree, 1, status );
+
+/* Invert literal constant unit multipliers. This is because a constant of
+   say 1000 for a unit of "m" means "multiply the value in metres by 1000",
+   but a unit string of "1000 m" means "value in units of 1000 m" (i.e.
+   *divide* the value in metres by 1000). */
+      InvertConstants( &in_tree, status );
+
+/* Convert the tree into string form. */
+      result = MakeExp( in_tree, 2, 1, status );
+
+/* If the result is a constant value, return a blank string. */
+      if( 1 == astSscanf( result, "%lg", &dval ) ) {
+         *((char *) result) = 0;
+      }
+
+/* Free the tree. */
+      in_tree = FreeTree( in_tree, status );
+
+   } else {
+      astError( AST__BADUN, "astUnitNormaliser: Error parsing input "
+                "units string '%s'.", status, in );
+   }
+
+/* Return the result */
+   return result;
+}
+
+static int Ustrcmp( const char *a, const char *b, int *status ){
+/*
+*  Name:
+*     Ustrcmp
+
+*  Purpose:
+*     A case blind version of strcmp.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     int Ustrcmp( const char *a, const char *b, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     Returns 0 if there are no differences between the two strings, and 1
+*     otherwise. Comparisons are case blind.
+
+*  Parameters:
+*     a
+*        Pointer to first string.
+*     b
+*        Pointer to second string.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     Zero if the strings match, otherwise one.
+
+*  Notes:
+*     -  This function does not consider the sign of the difference between
+*     the two strings, whereas "strcmp" does.
+*     -  This function attempts to execute even if an error has occurred.
+
+*/
+
+/* Local Variables: */
+   const char *aa;         /* Pointer to next "a" character */
+   const char *bb;         /* Pointer to next "b" character */
+   int ret;                /* Returned value */
+
+/* Initialise the returned value to indicate that the strings match. */
+   ret = 0;
+
+/* Initialise pointers to the start of each string. */
+   aa = a;
+   bb = b;
+
+/* Loop round each character. */
+   while( 1 ){
+
+/* We leave the loop if either of the strings has been exhausted. */
+      if( !(*aa ) || !(*bb) ){
+
+/* If one of the strings has not been exhausted, indicate that the
+   strings are different. */
+         if( *aa || *bb ) ret = 1;
+
+/* Break out of the loop. */
+         break;
+
+/* If neither string has been exhausted, convert the next characters to
+   upper case and compare them, incrementing the pointers to the next
+   characters at the same time. If they are different, break out of the
+   loop. */
+      } else {
+
+         if( toupper( (int) *(aa++) ) != toupper( (int) *(bb++) ) ){
+            ret = 1;
+            break;
+         }
+
+      }
+
+   }
+
+/* Return the result. */
+   return ret;
+
+}
+
+static int Ustrncmp( const char *a, const char *b, size_t n, int *status ){
+/*
+*  Name:
+*     Ustrncmp
+
+*  Purpose:
+*     A case blind version of strncmp.
+
+*  Type:
+*     Private function.
+
+*  Synopsis:
+*     #include "unit.h"
+*     int Ustrncmp( const char *a, const char *b, size_t n, int *status )
+
+*  Class Membership:
+*     Unit member function.
+
+*  Description:
+*     Returns 0 if there are no differences between the first "n"
+*     characters of the two strings, and 1 otherwise. Comparisons are
+*     case blind.
+
+*  Parameters:
+*     a
+*        Pointer to first string.
+*     b
+*        Pointer to second string.
+*     n
+*        The maximum number of characters to compare.
+*     status
+*        Pointer to the inherited status variable.
+
+*  Returned Value:
+*     Zero if the strings match, otherwise one.
+
+*  Notes:
+*     -  This function does not consider the sign of the difference between
+*     the two strings, whereas "strncmp" does.
+*     -  This function attempts to execute even if an error has occurred.
+
+*/
+
+/* Local Variables: */
+   const char *aa;         /* Pointer to next "a" character */
+   const char *bb;         /* Pointer to next "b" character */
+   int i;                  /* Character index */
+   int ret;                /* Returned value */
+
+
+/* Initialise the returned value to indicate that the strings match. */
+   ret = 0;
+
+/* Check pointer have been supplied. */
+   if( !a || !b ) return ret;
+
+/* Initialise pointers to the start of each string. */
+   aa = a;
+   bb = b;
+
+/* Compare up to "n" characters. */
+   for( i = 0; i < (int) n; i++ ){
+
+/* We leave the loop if either of the strings has been exhausted. */
+      if( !(*aa ) || !(*bb) ){
+
+/* If one of the strings has not been exhausted, indicate that the
+   strings are different. */
+         if( *aa || *bb ) ret = 1;
+
+/* Break out of the loop. */
+         break;
+
+/* If neither string has been exhausted, convert the next characters to
+   upper case and compare them, incrementing the pointers to the next
+   characters at the same time. If they are different, break out of the
+   loop. */
+      } else {
+
+         if( toupper( (int) *(aa++) ) != toupper( (int) *(bb++) ) ){
+            ret = 1;
+            break;
+         }
+
+      }
+
+   }
+
+/* Return the result. */
+   return ret;
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
+/* The rest of this file contains functions which are of use for debugging
+   this module. They are usually commented out.
+
+static const char *DisplayTree( UnitNode *node, int ind ) {
+   int i;
+   char buf[200];
+   const char *result;
+   char *a;
+   const char *arg[ 2 ];
+   int rl;
+   int slen;
+   const opsym[ 100 ];
+
+   result = "";
+
+   for( i = 0; i < ind; i++ ) buf[ i ] = ' ';
+   buf[ ind ] = 0;
+
+   if( !node ) {
+      printf( "%s <null>\n", buf );
+   } else {
+
+      printf( "%s Code: '%s' (%d)\n", buf, OpName( node->opcode ), node->opcode );
+      printf( "%s Narg: %d\n", buf, node->narg );
+      printf( "%s Constant: %g\n", buf, node->con );
+      printf( "%s Name: %s\n", buf, node->name?node->name:"" );
+      printf( "%s Unit: %s\n", buf, node->unit?node->unit->sym:"" );
+      printf( "%s Mult: %s\n", buf, node->mult?node->mult->sym:"" );
+
+      OpSym( node, opsym );
+      slen = strlen( opsym );
+      rl = slen;
+
+      if( node->narg == 0 ) {
+         result = astMalloc( rl + 1 );
+         if( astOK ) strcpy( (char *) result, opsym );
+
+      } else if( node->narg == 1 ) {
+         rl += 2;
+         printf( "%s Arg 0:\n", buf );
+         arg[ 0 ] = DisplayTree( (node->arg)[ 0 ], ind + 2 );
+         rl += strlen( arg[ 0 ] );
+
+         result = astMalloc( rl + 1 );
+         if( astOK ) {
+            a = (char *) result;
+            strcpy( a, opsym );
+            a += slen;
+            *(a++) = '(';
+            strcpy( a, arg[0] );
+            a += strlen( arg[ 0 ] );
+            *(a++) = ')';
+         }
+
+      } else {
+         rl += 4;
+         for( i = 0; i < node->narg; i++ ) {
+            printf( "%s Arg %d:\n", buf, i );
+            arg[ i ] = DisplayTree( (node->arg)[ i ], ind + 2 );
+            rl += strlen( arg[ i ] );
+         }
+
+         result = astMalloc( rl + 1 );
+         if( astOK ) {
+            a = (char *) result;
+            *(a++) = '(';
+            strcpy( a, arg[0] );
+            a += strlen( arg[ 0 ] );
+            *(a++) = ')';
+            strcpy( a, opsym );
+            a += slen;
+            *(a++) = '(';
+            strcpy( a, arg[1] );
+            a += strlen( arg[ 1 ] );
+            *(a++) = ')';
+         }
+      }
+   }
+
+   if( !astOK ) {
+      astFree( (void *) result );
+      result = "";
+   }
+
+   return result;
+}
+
+static const char *OpName( Oper op ) {
+   const char *name;
+
+   if( op ==  OP_LDCON ) {
+      name = "LDCON";
+   } else if( op ==  OP_LDVAR ) {
+      name = "LDVAR";
+   } else if( op ==  OP_LOG ) {
+      name = "LOG";
+   } else if( op ==  OP_LN ) {
+      name = "LN";
+   } else if( op ==  OP_EXP ) {
+      name = "EXP";
+   } else if( op ==  OP_SQRT ) {
+      name = "SQRT";
+   } else if( op ==  OP_POW ) {
+      name = "POW";
+   } else if( op ==  OP_DIV ) {
+      name = "DIV";
+   } else if( op ==  OP_MULT ) {
+      name = "MULT";
+   } else if( op ==  OP_LDPI ) {
+      name = "LDPI";
+   } else if( op ==  OP_LDE ) {
+      name = "LDE";
+   } else if( op ==  OP_NULL ) {
+      name = "NULL";
+   } else {
+      name = "<unknown op code>";
+   }
+
+   return name;
+}
+
+static void OpSym( UnitNode *node, char *buff ) {
+   const char *sym = NULL;
+
+   if( node->con != AST__BAD ) {
+      sprintf( buff, "%g", node->con );
+
+   } else if( node->opcode ==  OP_LDVAR ) {
+      sym = node->name;
+
+   } else if( node->opcode ==  OP_LOG ) {
+      sym = "log";
+
+   } else if( node->opcode ==  OP_LN ) {
+      sym = "ln";
+
+   } else if( node->opcode ==  OP_EXP ) {
+      sym = "exp";
+
+   } else if( node->opcode ==  OP_SQRT ) {
+      sym = "sqrt";
+
+   } else if( node->opcode ==  OP_POW ) {
+      sym = "**";
+
+   } else if( node->opcode ==  OP_DIV ) {
+      sym = "/";
+
+   } else if( node->opcode ==  OP_MULT ) {
+      sym = "*";
+
+   } else if( node->opcode ==  OP_NULL ) {
+      sym = "NULL";
+
+   } else {
+      sym = "<unknown op code>";
+   }
+
+   if( sym ) strcpy( buff, sym );
+}
+
+static const char *TreeExp( UnitNode *node ) {
+   char buff[ 100 ];
+   char buff2[ 100 ];
+
+   if( node->narg == 0 ) {
+      OpSym( node, buff );
+
+   } else if( node->narg == 1 ) {
+      OpSym( node, buff2 );
+      sprintf( buff, "%s(%s)", buff2, TreeExp( node->arg[ 0 ] ) );
+
+   } else if( node->narg == 2 ) {
+      OpSym( node, buff2 );
+      sprintf( buff, "(%s)%s(%s)", TreeExp( node->arg[ 0 ] ), buff2,
+                                   TreeExp( node->arg[ 1 ] ) );
+   }
+
+   return astStore( NULL, buff, strlen( buff ) + 1 );
+}
+
+*/
+
+
+
+