diff options
Diffstat (limited to 'ast/unit.c')
| -rw-r--r-- | ast/unit.c | 6218 |
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 ); +} + +*/ + + + + |