1 /*
2 * Name:
3 * unit.c
4
5 * Purpose:
6 * Implement unit conversion functions.
7
8 * Description:
9 * This file implements the Unit module which is used for identifying
10 * units and transforming between them. It follows the recommendations
11 * for unit handling contained within FITS WCS paper I (Greisen &
12 * Calabretta). All methods have protected access.
13
14 * Methods:
15 * astUnitMapper: Create a Mapping between two systems of units.
16 * astUnitLabel: Returns a label for a given unit symbol.
17
18 * Copyright:
19 * Copyright (C) 1997-2006 Council for the Central Laboratory of the
20 * Research Councils
21
22 * Licence:
23 * This program is free software: you can redistribute it and/or
24 * modify it under the terms of the GNU Lesser General Public
25 * License as published by the Free Software Foundation, either
26 * version 3 of the License, or (at your option) any later
27 * version.
28 *
29 * This program is distributed in the hope that it will be useful,
30 * but WITHOUT ANY WARRANTY; without even the implied warranty of
31 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32 * GNU Lesser General Public License for more details.
33 *
34 * You should have received a copy of the GNU Lesser General
35 * License along with this program. If not, see
36 * <http://www.gnu.org/licenses/>.
37
38 * Authors:
39 * DSB: D.S. Berry (Starlink)
40
41 * History:
42 * 10-DEC-2002 (DSB):
43 * Original version.
44 * 10-FEB-2004 (DSB):
45 * Added debug conditional code to keep track of memory leaks.
46 * 15-JUL-2004 (DSB):
47 * In astUnitMapper: if no Mapping can be found from input to
48 * output units (e.g. because fo the less than perfect simplication
49 * algorithm in SimplifyTree), try finding a Mapping from output to
50 * input units and inverting the result.
51 * 14-DEC-2004 (DSB):
52 * In CreateTree, move the invocation of FixConstants from after
53 * InvertConstants to before InvertConstants. This is because
54 * InvertConstants ignores nodes which contain all constant
55 * arguments. This results in constants not being inverted in
56 * expressions such as "1/60 deg" (because all arguments are
57 * constant in the the "1/60" node).
58 * 18-JAN-2005 (DSB):
59 * Fix memory leaks.
60 * 2-FEB-2005 (DSB):
61 * - Avoid using astStore to allocate more storage than is supplied
62 * in the "data" pointer. This can cause access violations since
63 * astStore will then read beyond the end of the "data" area.
64 * 15-FEB-2005 (DSB):
65 * - Modified CleanExp to fix up some common units mistakes.
66 * 21-FEB-2005 (DSB):
67 * - Modified CleanExp to accept <word><digit> as equivalent to
68 * <word>^<digit>.
69 * - Modified MakeTree to do case insensitive checking if case
70 * sensitive checking failsto produce a match to a multiplier/unit
71 * symbol combination. If this still produces no match, do a case
72 * insensitive match for multiplier/unit label.
73 * 1-MAR-2006 (DSB):
74 * Replace astSetPermMap within DEBUG blocks by astBeginPM/astEndPM.
75 * 6-APR-2006 (DSB):
76 * Modify CleanExp to convert "MJY/STER" to standard form ("MJy/sr").
77 * 7-JUL-2006 (DSB):
78 * Correct initialisation of "word" flag in CleanExp.
79 * 17-MAY-2007 (DSB):
80 * Simplify the units string returned by astUnitNormaliser.
81 * - Fix indexing bug in CombineFactors.
82 * 26-MAY-2007 (DSB):
83 * - Correct error reporting in astUNitNormaliser.
84 * - Fix bug in CleanExp that caused (for example) "-2" to be
85 * converted to "^-2".
86 * 2-DEC-2008 (DSB):
87 * Correct memory allocation bug in CleanExp.
88 * 6-MAY-2011 (DSB):
89 * Include "adu" as basic unit.
90 * 9-MAY-2011 (DSB):
91 * Change "A" to be Ampere (as defined by FITS-WCS paper 1) rather
92 * than "Angstrom".
93 */
94
95 /* Module Macros. */
96 /* ============== */
97 /* Define the astCLASS macro (even although this is not a class
98 implementation) to obtain access to the protected error and memory
99 handling functions. */
100 #define astCLASS
101 #define PI 3.141592653589793
102 #define E 2.718281828459045
103
104 /* Macro which returns the nearest integer to a given floating point
105 value. */
106 #define NINT(x) (int)((x)+(((x)>0.0)?0.5:-0.5))
107
108 /* Macros which return the maximum and minimum of two values. */
109 #define MAX(aa,bb) ((aa)>(bb)?(aa):(bb))
110 #define MIN(aa,bb) ((aa)<(bb)?(aa):(bb))
111
112 /* Macro to check for equality of floating point values. We cannot
113 compare bad values directory because of the danger of floating point
114 exceptions, so bad values are dealt with explicitly. */
115 #define EQUAL(aa,bb) (((aa)==AST__BAD)?(((bb)==AST__BAD)?1:0):(((bb)==AST__BAD)?0:(fabs((aa)-(bb))<=1.0E5*MAX((fabs(aa)+fabs(bb))*DBL_EPSILON,DBL_MIN))))
116
117 /* Macro identifying a character as lower or upper case letter, digit or
118 + or -. */
119 #define ISWORD(c) (isalnum(c)||((c)=='+')||((c)=='-'))
120
121 /* The number of basic dimension quantities used for dimensional analysis.
122 In addition to the usual M, L and T, this includes pseudo-dimensions
123 describing strictly dimensionless quantities such as plane angle,
124 magnitude, etc. */
125 #define NQUANT 10
126
127 /* Include files. */
128 /* ============== */
129 /* Interface definitions. */
130 /* ---------------------- */
131 #include "error.h"
132 #include "memory.h"
133 #include "pointset.h"
134 #include "mapping.h"
135 #include "unitmap.h"
136 #include "zoommap.h"
137 #include "mathmap.h"
138 #include "unit.h"
139
140 /* Error code definitions. */
141 /* ----------------------- */
142 #include "ast_err.h" /* AST error codes */
143
144 /* C header files. */
145 /* --------------- */
146 #include <string.h>
147 #include <stdio.h>
148 #include <ctype.h>
149 #include <limits.h>
150 #include <math.h>
151
152 #ifdef THREAD_SAFE
153 #include <pthread.h>
154 #endif
155
156 /* Module Type Definitions. */
157 /* ======================== */
158
159 /* This declaration enumerates the codes for the operations which may
160 legally be included in a units expression. */
161 typedef enum {
162 OP_LDCON, /* Load constant */
163 OP_LDVAR, /* Load variable */
164 OP_LOG, /* Common logarithm */
165 OP_LN, /* Natural logarithm */
166 OP_EXP, /* Natural exponential */
167 OP_SQRT, /* Square root */
168 OP_POW, /* Exponentiate */
169 OP_DIV, /* Division */
170 OP_MULT, /* Multiplication */
171 OP_LDPI, /* Load constant PI */
172 OP_LDE, /* Load constant E */
173 OP_NULL /* Null operation */
174 } Oper;
175
176 /* A structure describing a standard multiplier prefix. */
177 typedef struct Multiplier {
178 const char *label; /* Multipler label string (null terminated) */
179 const char *sym; /* Multipler symbol string (null terminated) */
180 int symlen; /* Length of symbol (without trailing null ) */
181 int lablen; /* Length of label (without trailing null ) */
182 double scale; /* The scale factor associated with the prefix */
183 struct Multiplier *next; /* Next Multiplier in linked list */
184 } Multiplier;
185
186 /* A structure describing a single node in a tree representation of a
187 single units expression. */
188 typedef struct UnitNode {
189 Oper opcode; /* Code for operation performed by this node */
190 int narg; /* No. of arguments used by the operation */
191 struct UnitNode **arg; /* Array of pointers to argument nodes */
192 double con; /* Constant to be loaded by OP_LDCON operations */
193 struct KnownUnit *unit; /* Known unit referred to by OP_LDVAR nodes */
194 Multiplier *mult; /* Multiplier used by OP_LDVAR nodes */
195 const char *name; /* User-defined unit referred to by OP_LDVAR
196 nodes (no multiplier prefix included) */
197 } UnitNode;
198
199 /* A structure describing a known unit. */
200 typedef struct KnownUnit {
201 const char *sym; /* Unit symbol string (null terminated) */
202 const char *label; /* Unit label string (null terminated) */
203 int symlen; /* Length of symbol (without trailing null ) */
204 int lablen; /* Length of label (without trailing null ) */
205 struct UnitNode *head; /* Head of definition tree (NULL for basic units) */
206 struct KnownUnit *next; /* Next KnownUnit in linked list */
207 struct KnownUnit *use; /* KnownUnit to be used in place of this one */
208 } KnownUnit;
209
210 /* Module Variables. */
211 /* ================= */
212
213 /* A pointer to the KnownUnit structure at the head of a linked list of
214 such structures containing definitions of all known units. */
215 static KnownUnit *known_units = NULL;
216
217 /* An array of pointers to KnownUnits which list the basic quantities
218 used in dimensional analysis. */
219 static KnownUnit *quant_units[ NQUANT ];
220
221 /* A pointer to the Multiplier structure at the head of a linked list of
222 such structures containing definitions of all known multipliers. */
223 static Multiplier *multipliers = NULL;
224
225 /* Set up mutexes */
226 #ifdef THREAD_SAFE
227
228 static pthread_mutex_t mutex1 = PTHREAD_MUTEX_INITIALIZER;
229 #define LOCK_MUTEX1 pthread_mutex_lock( &mutex1 );
230 #define UNLOCK_MUTEX1 pthread_mutex_unlock( &mutex1 );
231
232 static pthread_mutex_t mutex2 = PTHREAD_MUTEX_INITIALIZER;
233 #define LOCK_MUTEX2 pthread_mutex_lock( &mutex2 );
234 #define UNLOCK_MUTEX2 pthread_mutex_unlock( &mutex2 );
235
236 #else
237
238 #define LOCK_MUTEX1
239 #define UNLOCK_MUTEX1
240
241 #define LOCK_MUTEX2
242 #define UNLOCK_MUTEX2
243
244 #endif
245
246 /* Prototypes for Private Functions. */
247 /* ================================= */
248 static AstMapping *MakeMapping( UnitNode *, int * );
249 static KnownUnit *GetKnownUnits( int, int * );
250 static Multiplier *GetMultipliers( int * );
251 static UnitNode *ConcatTree( UnitNode *, UnitNode *, int * );
252 static UnitNode *CopyTree( UnitNode *, int * );
253 static UnitNode *CreateTree( const char *, int, int, int * );
254 static UnitNode *FixUnits( UnitNode *, UnitNode *, int * );
255 static UnitNode *FreeTree( UnitNode *, int * );
256 static UnitNode *MakeTree( const char *, int, int, int * );
257 static UnitNode *MakeLabelTree( const char *, int, int * );
258 static UnitNode *NewNode( UnitNode *, Oper, int * );
259 static UnitNode *CombineFactors( UnitNode **, double *, int, double, int * );
260 static const char *CleanExp( const char *, int * );
261 static int EndsWith( const char *, int, const char *, int * );
262 static int CmpTree( UnitNode *, UnitNode *, int, int * );
263 static void FixConstants( UnitNode **, int, int * );
264 static void InvertConstants( UnitNode **, int * );
265 static UnitNode *InvertTree( UnitNode *, UnitNode *, int * );
266 static void LocateUnits( UnitNode *, UnitNode ***, int *, int * );
267 static void MakeKnownUnit( const char *, const char *, const char *, int * );
268 static void MakeUnitAlias( const char *, const char *, int * );
269 static void RemakeTree( UnitNode **, int * );
270 static int SimplifyTree( UnitNode **, int, int * );
271 static int ComplicateTree( UnitNode **, int * );
272 static int ReplaceNode( UnitNode *, UnitNode *, UnitNode *, int * );
273 static void FindFactors( UnitNode *, UnitNode ***, double **, int *, double *, int * );
274 static const char *MakeExp( UnitNode *, int, int, int * );
275 static int DimAnal( UnitNode *, double[NQUANT], double *, int * );
276 static int Ustrcmp( const char *, const char *, int * );
277 static int Ustrncmp( const char *, const char *, size_t, int * );
278 static int SplitUnit( const char *, int, const char *, int, Multiplier **, int *, int * );
279 static UnitNode *ModifyPrefix( UnitNode *, int * );
280 static int ConStart( const char *, double *, int *, int * );
281
282 /* Debug functions...
283 static const char *DisplayTree( UnitNode *, int );
284 static void OpSym( UnitNode *, char * );
285 static const char *OpName( Oper );
286 static const char *TreeExp( UnitNode * );
287 */
288
289 /* Function implementations. */
290 /* ========================= */
CleanExp(const char * exp,int * status)291 static const char *CleanExp( const char *exp, int *status ) {
292 /*
293 * Name:
294 * CleanExp
295
296 * Purpose:
297 * Produce a clean copy of a units expression.
298
299 * Type:
300 * Private function.
301
302 * Synopsis:
303 * #include "unit.h"
304 * const char *CleanExp( const char *exp )
305
306 * Class Membership:
307 * Unit member function.
308
309 * Description:
310 * This function returns a pointer to dynamic memory containing a
311 * cleaned copy of the supplied units expression. Cleaning consists of
312 * the following operations:
313 * - removal of leading and trailing white space.
314 * - replacement of multiple adjacent spaces by a single space
315 * - removal of spaces adjacent to a parenthesis
316 * - removal of spaces adjacent to a binary operator
317 * - translates various common non-standard units into equivalent
318 * standard units.
319 *
320 * Such carefull handling of spaces is necessary since a space is
321 * recognised by the MakeTree function as a multiplication operator.
322
323 * Parameters:
324 * exp
325 * A pointer to the expression to be cleaned.
326
327 * Returned Value:
328 * A pointer to the cleaned expression, which should be freed using
329 * astFree when no longer needed.
330
331 * Notes:
332 * - This function returns NULL if it is invoked with the global error
333 * status set, or if it should fail for any reason.
334 */
335
336 /* Local Variables: */
337 char **tok;
338 char *p;
339 char *r;
340 char *result;
341 char *s;
342 char *t;
343 char *w;
344 const char *start;
345 int i;
346 int l;
347 int len;
348 char *tt;
349 int ntok;
350 int po;
351 int ps;
352 int word;
353
354 /* Initialise */
355 result = NULL;
356
357 /* Check inherited status */
358 if( !astOK ) return result;
359
360 /* Split the supplied string up into tokens. Each block of contiguous
361 alphanumeric characters is a token. Each contiguous block of
362 non-alphanumerical characters is also a token. The + and - signs are
363 counted as alphanumeric. */
364 start = exp;
365 p = (char *) exp - 1;
366 word = ISWORD( *( p + 1 ) );
367 ntok = 0;
368 tok = NULL;
369 while( *(++p) ){
370 if( word ) {
371 if( !ISWORD( *p ) ) {
372 l = p - start;
373 t = astStore( NULL, start, l + 1 );
374 if( t ) t[ l ] = 0;
375 tok = astGrow( tok, ntok + 1, sizeof( char * ) );
376 if( tok ) tok[ ntok++ ] = t;
377 start = p;
378 word = 0;
379 }
380 } else {
381 if( ISWORD( *p ) ) {
382 l = p - start;
383 t = astStore( NULL, start, l + 1 );
384 if( t ) t[ l ] = 0;
385 tok = astGrow( tok, ntok + 1, sizeof( char * ) );
386 if( tok ) tok[ ntok++ ] = t;
387 start = p;
388 word = 1;
389 }
390 }
391 }
392
393 l = p - start;
394 t = astStore( NULL, start, l + 1 );
395 if( t ) t[ l ] = 0;
396 tok = astGrow( tok, ntok + 1, sizeof( char * ) );
397 if( tok ) tok[ ntok++ ] = t;
398
399 /* Check the tokens for known non-standard unit syntax, and replace with the
400 equivalent standard syntax. Starlink SPLAT has a class called UnitUtilities
401 which has more of these common units mistakes. AST has to be a bit
402 more conservative than SPLAT though because of its wider remit. */
403 len = 0;
404 tt = NULL;
405 for( i = 0; i < ntok; i++ ) {
406 t = tok[ i ];
407 l = strlen( t );
408 tt = astStore( tt, t, l + 1 );
409
410 /* Any alphabetical word followed by a digit is taken as <word>^<digit>.
411 Any alphabetical word followed by a sign and a digit is taken as
412 <word>^<sign><digit>. */
413 if( l > 1 && *t != '-' && *t != '+' &&
414 strcspn( t, "0123456789" ) == l - 1 ) {
415 tok[ i ] = astMalloc( l + 2 );
416 if( tok[ i ] ) {
417 strcpy( tok[ i ], t );
418 w = t + l - 2;
419 if( *w != '+' && *w != '-' ) {
420 tok[ i ][ l - 1 ] = '^';
421 strcpy( tok[ i ] + l, t + l - 1 );
422 } else {
423 tok[ i ][ l - 2 ] = '^';
424 strcpy( tok[ i ] + l - 1, t + l - 2 );
425 }
426 t = astFree( t );
427 }
428 l++;
429
430 /* If the word ends with "micron" change to "(<start>m*1.0E-6)". Should be OK
431 for things like "Kmicron". */
432 } else if( ( s = strstr( t, "micron" ) ) ) {
433 tok[ i ] = astMalloc( s - t + 11 );
434 if( tok[ i ] ) {
435 w = tok[ i ];
436 *(w++) = '(';
437 if( s > t ) {
438 strncpy( w, t, s - t );
439 w += s - t;
440 }
441 strcpy( w, "m*1.0E-6)" );
442 l = s - t + 11;
443 t = astFree( t );
444 }
445
446 /* Convert "STER" to "sr". */
447 } else if( !Ustrcmp( t, "STER", status ) ) {
448 tok[ i ] = astStore( NULL, "sr", 3 );
449 l = 2;
450 t = astFree( t );
451
452 /* If the word ends with "JY" and is preceded by a single character, change
453 to "<start>Jy". Should be OK for things like "MJY". */
454 } else if( l == 3 && !strcmp( t + 1, "JY" ) ) {
455 tok[ i ][ 2 ] = 'y';
456
457 /* If the word begins with "nano" (case-insensitive) change "nano" to
458 "n". Such changes are usually handled by SplitUnit, but we need to
459 handle this as a special case here since scanf seems to read "nan" as
460 a string representation of NaN. */
461 } else if( !Ustrncmp( t, "nano", 4, status ) ) {
462 tok[ i ] = astStore( NULL, t + 3, l - 2 );
463 if( tok[ i ] ) {
464 *(tok[ i ]) = 'n';
465 t = astFree( t );
466 }
467 l -= 3;
468 }
469
470 /* Update the total length of the string. */
471 len += l;
472 }
473 tt = astFree( tt );
474
475 /* Concatentate the tokens into a single string, freeing the individual
476 strings. */
477 result = astMalloc( len + 1 );
478 if( result ) {
479 p = result;
480 for( i = 0; i < ntok; i++ ) {
481 len = strlen( tok[ i ] );
482 memcpy( p, tok[ i ], len );
483 p += len;
484 tok[ i ] = astFree( tok[ i ] );
485 }
486 *p = 0;
487 tok = astFree( tok );
488
489 /* Now do other cleaning.
490 ---------------------- */
491
492 /* Initialise a pointer to the previous character read from the string. */
493 r = result - 1;
494
495 /* Initialise a pointer to the next character to be written to the string. */
496 w = result;
497
498 /* Pretend the previous character written to the string was a space. */
499 ps = 1;
500
501 /* Read all the supplied string, copying it to earlier parts of the
502 string discarding leading spaces and multiple adjacent embedded spaces in
503 the process. */
504 while( *(++r) ) {
505
506 /* If the character read is a space, only write it to the string if the
507 previous character written was not a space (in which case set a flag
508 to indicate that the previous character written to the string is now a
509 space). */
510 if( isspace( *r ) ) {
511 if( !ps ) {
512 *(w++) = *r;
513 ps = 1;
514 }
515
516 /* Write all non-space characters to the string, and clear the flag which
517 indicates if the previous character written to the string was a space. */
518 } else {
519 *(w++) = *r;
520 ps = 0;
521 }
522 }
523
524 /* If the last character written to the string was a space, reduce the
525 length of the string by one since we do not want any trailing spaces. */
526 if( ps ) w--;
527
528 /* Terminate the string. */
529 *w = 0;
530
531 /* We now need to pass through the string again, this time removing any
532 spaces which are adjacent to a binary operator or a parenthesis. */
533 r = result - 1;
534 w = result;
535 ps = 0;
536 po = 0;
537 while( *(++r) ) {
538
539 /* If the current character is a space, only write it if the previous
540 written character was not an operator or parenthesis. */
541 if( isspace( *r ) ) {
542 if( !po ) {
543 *(w++) = *r;
544 po = 1;
545 ps = 1;
546 }
547
548 /* If the current character is an operator or parenthesis, back up one
549 character before writing it out if the previous written character was
550 a space. */
551 } else if( *r == '*' || *r == '/' || *r == '^' || *r == '.' ||
552 *r == ')' || *r == '(' ) {
553 if( ps ) w--;
554 *(w++) = *r;
555 po = 1;
556 ps = 0;
557
558 /* If the current character is not a space and not an operator symbol,
559 just write it out. */
560 } else {
561 *(w++) = *r;
562 po = 0;
563 ps = 0;
564 }
565 }
566
567 /* Terminate the string. */
568 if( ps ) w--;
569 *w = 0;
570
571 }
572
573 /* Return the result. */
574 return (const char *) result;
575 }
576
CmpTree(UnitNode * tree1,UnitNode * tree2,int exact,int * status)577 static int CmpTree( UnitNode *tree1, UnitNode *tree2, int exact, int *status ) {
578 /*
579 * Name:
580 * CmpTree
581
582 * Purpose:
583 * Compares two trees of UnitNodes.
584
585 * Type:
586 * Private function.
587
588 * Synopsis:
589 * #include "unit.h"
590 * int CmpTree( UnitNode *tree1, UnitNode *tree2, int exact, int *status )
591
592 * Class Membership:
593 * Unit member function.
594
595 * Description:
596 * This function returns a zero value if the two trees are
597 * equivalent. This requires the trees to have identical structure
598 * except that, if "exact" is zero, arguments for OP_MULT nodes can
599 * be swapped.
600 *
601 * If the trees are not equivalent then a value of +1 or -1 is returned
602 * depending on whether tree1 should be placed before or after tree2
603 * in a sorted list of trees.
604
605 * Parameters:
606 * tree1
607 * A pointer to the UnitNode at the head of the first tree.
608 * tree2
609 * A pointer to the UnitNode at the head of the second tree.
610 * exact
611 * If non-zero, then OP_MULT nodes must have their arguments the
612 * same way round in order for the OP_MULT nodes to match. Otherwise,
613 * OP_MULT nodes with equivalent arguments match even if the
614 * arguments are swapped.
615 * status
616 * Pointer to the inherited status variable.
617
618 * Returned Value:
619 * Zero if the two trees are equal. +1 if tree1 should be placed before
620 * tree2 in a sorted list of trees. -1 if tree1 should be placed after
621 * tree2 in a sorted list of trees.
622
623 * Notes:
624 * - Zero is returned if an error has already occurred, or
625 * if this function fails for any reason.
626
627 */
628
629 /* Local Variables: */
630 int result;
631 int i;
632 Oper op;
633
634 /* Initialise. */
635 result = 0;
636
637 /* Check inherited status. */
638 if( !astOK ) return result;
639
640 /* If the op codes differ, compare them as integers. */
641 op = tree1->opcode;
642 if( op != tree2->opcode ) {
643 result = ( op > tree2->opcode ) ? 1: -1;
644
645 /* If both supplied nodes are OP_LDVAR nodes, compare the associated names. */
646 } else if( op == OP_LDVAR ){
647 result = strcmp( tree1->name, tree2->name );
648
649 /* If both supplied nodes are constant nodes, compare the constant values. */
650 } else if( tree1->con != AST__BAD ){
651 result = EQUAL( tree1->con, tree2->con ) ? 0 : (
652 ( tree1->con > tree2->con ) ? 1 : -1 );
653
654 /* Otherwise, compare the arguments for the node. */
655 } else {
656 for( i = 0; i < tree1->narg; i++ ) {
657 result = CmpTree( tree1->arg[ i ], tree2->arg[ i ], exact, status );
658 if( result ) break;
659 }
660
661 /* If the head nodes of the two trees are OP_MULT nodes, and the above
662 check determined they are different, this may be just because they
663 have their operands swapped. If "exact" si zero, this is considered an
664 insignificant difference between the two trees which we should ignore.
665 To check for this try comparing the arguments again, this time swapping
666 the arguments of tree2. */
667 if( result && op == OP_MULT && !exact ) {
668 for( i = 0; i < tree1->narg; i++ ) {
669 result = CmpTree( tree1->arg[ i ], tree2->arg[ 1 - i ], 0, status );
670 if( result ) break;
671 }
672 }
673 }
674
675 /* If an error has occurred, return zero. */
676 if( !astOK ) result = 0;
677
678 /* Return the answer. */
679 return result;
680 }
681
CombineFactors(UnitNode ** factors,double * powers,int nfactor,double coeff,int * status)682 static UnitNode *CombineFactors( UnitNode **factors, double *powers,
683 int nfactor, double coeff, int *status ) {
684 /*
685 * Name:
686 * CombineFactors
687
688 * Purpose:
689 * Create a tree which represents the product of the supplied factors.
690
691 * Type:
692 * Private function.
693
694 * Synopsis:
695 * #include "unit.h"
696 * UnitNode *CombineFactors( UnitNode **factors, double *powers,
697 * int nfactor, double coeff, int *status )
698
699 * Class Membership:
700 * Unit member function.
701
702 * Description:
703 * This function createa a tree of UnitNodes which represents the
704 * product of the supplied factors, and the supplied coefficient.
705 * The factors are sorted before being combined, using the sort order
706 * implemented by the CmpTree function.
707
708 * Parameters:
709 * factors
710 * A pointer to an array with "nfactor" elements, each element being
711 * a pointer to a UnitNode which is a factor of the required tree.
712 * On exit, the array is sorted.
713 * powers
714 * A pointer to an array with "nfactor" elements, each element being a
715 * double holding the power of the associated factor in "factors".
716 * On exit, the array reflects the sorting applied to "factors".
717 * nfactor
718 * The number of elements in the "factors" and "powers" arrays.
719 * coeff
720 * The overall coefficient to be applied to the product of the factors.
721 * status
722 * Pointer to the inherited status variable.
723
724 * Returned Value:
725 * A pointer to a UnitNode which is at the head of the new tree.
726
727 * Notes:
728 * - A NULL pointer is returned if an error has already occurred, or
729 * if this function fails for any reason.
730
731 */
732
733 /* Local Variables: */
734 UnitNode *result;
735 int i;
736 int j;
737 int jp;
738 int done;
739 UnitNode *ftmp;
740 UnitNode *node1;
741 UnitNode *node2;
742 UnitNode *pnode;
743 double ptmp;
744
745 /* Initialise. */
746 result = NULL;
747
748 /* Check inherited status. */
749 if( !astOK ) return result;
750
751 /* Sort the supplied list of factors, modifying the powers array
752 correspondingly. A simple bubblesort algorithm is used since there
753 will only be a handfull of factors. */
754 for( i = nfactor - 1; i > 0; i-- ) {
755 done = 1;
756 for( j = 0, jp = 1; j < i; j++, jp++ ) {
757 if( CmpTree( factors[ j ], factors[ jp ], 0, status ) > 0 ) {
758 ftmp = factors[ j ];
759 factors[ j ] = factors[ jp ];
760 factors[ jp ] = ftmp;
761
762 ptmp = powers[ j ];
763 powers[ j ] = powers[ jp ];
764 powers[ jp ] = ptmp;
765
766 done = 0;
767 }
768 }
769 if( done ) break;
770 }
771
772 /* The first root term of the returned tree is the coefficient, unless the
773 coefficient is 1.0, in which case it will be the first factor. */
774 if( coeff != 1.0 ) {
775 node1 = NewNode( NULL, OP_LDCON, status );
776 if( astOK ) node1->con = coeff;
777 } else {
778 node1 = NULL;
779 }
780
781 /* Loop through the factors. */
782 for( i = 0; i < nfactor; i++ ) {
783
784 /* If the power of this factor is zero, we ignore the factor. */
785 if( powers[ i ] != 0.0 ) {
786
787 /* If the power of this factor is one, we use the factor directly. */
788 if( EQUAL( powers[ i ], 1.0 ) ) {
789 node2 = CopyTree( factors[ i ], status );
790
791 /* Otherwise, for non-zero, non-unity powers, we create a POW node for
792 the factor. */
793 } else {
794 node2 = NewNode( NULL, OP_POW, status );
795 pnode = NewNode( NULL, OP_LDCON, status );
796 if( astOK ) {
797 pnode->con = powers[ i ];
798 node2->arg[ 0 ] = CopyTree( factors[ i ], status );
799 node2->arg[ 1 ] = pnode;
800 }
801 }
802
803 /* We now combine node1 and node2 using an OP_MULT node, which becomes
804 the "node1" for the next pass. On the first pass we may have no node1 (if
805 the supplied coefficient was 1.0), in which case we reserve the current
806 node2 as the node1 for the next pass. */
807 if( node1 ) {
808 result = NewNode( NULL, OP_MULT, status );
809 if( astOK ) {
810 result->arg[ 0 ] = node1;
811 result->arg[ 1 ] = node2;
812 node1 = result;
813 }
814 } else {
815 node1 = node2;
816 }
817 }
818 }
819
820 /* Ensure we have a node to return. */
821 if( astOK ) {
822 if( !result ) result = node1;
823 if( !result ) {
824 result = NewNode( NULL, OP_LDCON, status );
825 if( astOK ) result->con = 1.0;
826 }
827 }
828
829 /* If an error has occurred, free any new tree. */
830 if( !astOK ) result = FreeTree( result, status );
831
832 /* Return the answer. */
833 return result;
834 }
835
ComplicateTree(UnitNode ** node,int * status)836 static int ComplicateTree( UnitNode **node, int *status ) {
837 /*
838 * Name:
839 * ComplicateTree
840
841 * Purpose:
842 * Removes standardisations introduced by SimplifyTree.
843
844 * Type:
845 * Private function.
846
847 * Synopsis:
848 * #include "unit.h"
849 * int ComplicateTree( UnitNode **node )
850
851 * Class Membership:
852 * Unit member function.
853
854 * Description:
855 * This function modifies a tree of UnitNodes by removing standardisations
856 * introduced by SimplifyTree. The standardisations removed are ones
857 * which would make the corresponding algebraic expression (as produced
858 * by MakeExp) unnatural to a human reader.
859
860 * Parameters:
861 * node
862 * The address of a pointer to the UnitNode at the head of the tree
863 * which is to be complicated. On exit the supplied tree is freed and
864 * a pointer to a new tree is placed at the given address.
865
866 * Returned Value:
867 * Non-zero if any change was introduced into the tree.
868
869 */
870
871 /* Local Variables: */
872 int i;
873 UnitNode *newnode;
874 UnitNode *node1;
875 UnitNode *node2;
876 UnitNode *node3;
877 Oper op;
878 double con;
879 double fk;
880 int k;
881 int result;
882 double kcon;
883
884 /* Initialise */
885 result = 0;
886
887 /* Check inherited status. */
888 if( !astOK ) return result;
889
890 /* Initiallially, we have no replacement node. */
891 newnode = NULL;
892 node1 = NULL;
893 node3 = NULL;
894
895 /* Complicate the sub-trees corresponding to the arguments of the node at
896 the head of the supplied tree. */
897 for( i = 0; i < (*node)->narg; i++ ) {
898 if( ComplicateTree( &( (*node)->arg[ i ] ), status ) ) result = 1;
899 }
900
901 /* Now undo specific simplifications appropriate to the nature of the node at
902 the head of the tree. */
903 op = (*node)->opcode;
904
905 /* If the head is an OP_MULT node with a constant first argument and
906 a "LN" second argument, rearrange the nodes to represent ln(x**k) instead
907 of k*ln(x). If k is an integer multiple of "0.1/ln(10)" convert the "ln"
908 function into a "log" (base 10) function. Check for "k==1" in which
909 case we do not need a POW node. */
910 if( (*node)->opcode == OP_MULT ) {
911
912 con = (*node)->arg[ 0 ]->con;
913 if( con != AST__BAD && (*node)->arg[ 1 ]->opcode == OP_LN ) {
914 fk = 10.0*con*log( 10.0 );
915 k = NINT(fk);
916 if( EQUAL(fk,k) ) {
917 newnode = NewNode( NULL, OP_LOG, status );
918 con = k/10.0;
919 } else {
920 newnode = NewNode( NULL, OP_LN, status );
921 }
922
923 node2 = CopyTree( (*node)->arg[ 1 ]->arg[ 0 ], status );
924 if( !EQUAL( con, 1.0 ) ){
925 node1 = CopyTree( (*node)->arg[ 0 ], status );
926 node3 = NewNode( NULL, OP_POW, status );
927 }
928
929 if( astOK ) {
930 if( !EQUAL( con, 1.0 ) ){
931 node1->con = con;
932 node3->arg[ 0 ] = node2;
933 node3->arg[ 1 ] = node1;
934 newnode->arg[ 0 ] = node3;
935 } else {
936 newnode->arg[ 0 ] = node2;
937 }
938 }
939
940 /* Replace "(A**-1)*B" with "B/A" */
941 } else if( (*node)->arg[ 0 ]->opcode == OP_POW &&
942 EQUAL( (*node)->arg[ 0 ]->arg[ 1 ]->con, -1.0 )) {
943 newnode = NewNode( NULL, OP_DIV, status );
944 if( astOK ) {
945 newnode->arg[ 0 ] = CopyTree( (*node)->arg[ 1 ], status );
946 newnode->arg[ 1 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
947 }
948
949 /* Replace "B*(A**-1)" with "B/A" */
950 } else if( (*node)->arg[ 1 ]->opcode == OP_POW &&
951 EQUAL( (*node)->arg[ 1 ]->arg[ 1 ]->con, -1.0 )) {
952 newnode = NewNode( NULL, OP_DIV, status );
953 if( astOK ) {
954 newnode->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ], status );
955 newnode->arg[ 1 ] = CopyTree( (*node)->arg[ 1 ]->arg[ 0 ], status );
956 }
957
958 /* Convert (x**k)*(y**k) to (x*y)**k. */
959 } else if( (*node)->arg[ 0 ]->opcode == OP_POW &&
960 (*node)->arg[ 1 ]->opcode == OP_POW &&
961 EQUAL( (*node)->arg[ 0 ]->arg[ 1 ]->con,
962 (*node)->arg[ 1 ]->arg[ 1 ]->con )) {
963 newnode = NewNode( NULL, OP_POW, status );
964 node1 = NewNode( NULL, OP_MULT, status );
965 if( astOK ) {
966 node1->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
967 node1->arg[ 1 ] = CopyTree( (*node)->arg[ 1 ]->arg[ 0 ], status );
968 newnode->arg[ 0 ] = node1;
969 newnode->arg[ 1 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 1 ], status );
970 }
971
972 /* Convert c*sqrt(x) to sqrt((c**2)*x) (if c > 0). */
973 } else if( (kcon=(*node)->arg[ 0 ]->con) != AST__BAD &&
974 kcon > 0.0 && (*node)->arg[ 1 ]->opcode == OP_SQRT ) {
975 newnode = NewNode( NULL, OP_SQRT, status );
976 node1 = NewNode( NULL, OP_MULT, status );
977 node2 = NewNode( NULL, OP_LDCON, status );
978 if( astOK ) {
979 node2->con = kcon*kcon;
980 node1->arg[ 0 ] = node2;
981 node1->arg[ 1 ] = CopyTree( (*node)->arg[ 1 ]->arg[ 0 ], status );
982 newnode->arg[ 0 ] = node1;
983 }
984 }
985
986 /* If the head node is a POW node, replace "x**0.5" by sqrt(x) */
987 } else if( (*node)->opcode == OP_POW ) {
988 if( EQUAL( (*node)->arg[ 1 ]->con, 0.5 ) ) {
989 newnode = NewNode( NULL, OP_SQRT, status );
990 if( astOK ) {
991 newnode->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ], status );
992 }
993 }
994 }
995
996 /* If we have produced a new node which is identical to the old node,
997 free it. Otherwise, indicate we have made some changes. */
998 if( newnode ) {
999 if( !CmpTree( newnode, *node, 1, status ) ) {
1000 newnode = FreeTree( newnode, status );
1001 } else {
1002 result = 1;
1003 }
1004 }
1005
1006 /* If an error has occurred, free any new node. */
1007 if( !astOK ) {
1008 newnode = FreeTree( newnode, status );
1009 result = 0;
1010 }
1011
1012 /* If we have a replacement node, free the supplied tree and return a
1013 pointer to the new tree. */
1014 if( newnode ) {
1015 FreeTree( *node, status );
1016 *node = newnode;
1017 }
1018
1019 /* If the above produced some change, try simplifying (without
1020 re-introducing the standardisation we have just got rid of!) and
1021 then re-complicating the tree. */
1022 if( result ) {
1023 SimplifyTree( node, 0, status );
1024 ComplicateTree( node, status );
1025 }
1026
1027 /* Return the result. */
1028 return result;
1029 }
1030
ConcatTree(UnitNode * tree1,UnitNode * tree2,int * status)1031 static UnitNode *ConcatTree( UnitNode *tree1, UnitNode *tree2, int *status ) {
1032 /*
1033 * Name:
1034 * ConcatTree
1035
1036 * Purpose:
1037 * Concatenate two trees together.
1038
1039 * Type:
1040 * Private function.
1041
1042 * Synopsis:
1043 * #include "unit.h"
1044 * UnitNode *ConcatTree( UnitNode *tree1, UnitNode *tree2, int *status )
1045
1046 * Class Membership:
1047 * Unit member function.
1048
1049 * Description:
1050 * This function a pointer to the head of a new tree of UnitNodes which
1051 * is formed by feeding the output of "tree1" (i.e. the quantity
1052 * represented by the node at the head of tree1) into the (single)
1053 * input of "tree2" (i.e. the single OP_LDVAR Node containined within
1054 * tree2).
1055
1056 * Parameters:
1057 * tree1
1058 * A pointer to the first tree.
1059 * tree2
1060 * A pointer to the second tree. This should have no more than one
1061 * OP_LDVAR node.
1062 * status
1063 * Pointer to the inherited status variable.
1064
1065 * Returned Value:
1066 * A pointer to a UnitNode which is at the head of the new tree.
1067
1068 * Notes:
1069 * - If "tree2" contains zero units, a NULL pointer is returned but no
1070 * error is reported.
1071 * - If "tree2" contains more than one unit, an error is reported
1072 * error is reported.
1073 * - A NULL pointer is returned if an error has already occurred, or
1074 * if this function fails for any reason.
1075
1076 */
1077
1078 /* Local Variables: */
1079 UnitNode *result;
1080 UnitNode **units;
1081 int nunits;
1082
1083 /* Initialise. */
1084 result = NULL;
1085
1086 /* Check inherited status. */
1087 if( !astOK ) return result;
1088
1089 /* Produce a copy of tree2. */
1090 result = CopyTree( tree2, status );
1091
1092 /* Locate the OP_LDVAR node in the copy of tree2. */
1093 units = NULL;
1094 nunits = 0;
1095 LocateUnits( result, &units, &nunits, status );
1096
1097 /* If no OP_LDVAR nodes were found in tree2, we cannot concatenate the
1098 trees. */
1099 if( nunits > 0 ) {
1100
1101 /* Report an error if the number of pointers returned is larger than 1. */
1102 if( nunits > 1 && astOK ) {
1103 astError( AST__INTER, "ConcatTree(unit): tree2 uses %d units - "
1104 "should be 1 (internal AST programming error).", status, nunits );
1105 }
1106
1107 /* Replace the OP_LDVAR node in the copy of tree2 with a copy of tree1. */
1108 if( astOK ) {
1109
1110 /* If the node at the head of the supplied tree2 is the node to be
1111 replaced, just free the tree created earlier and return a copy of
1112 tree1. */
1113 if( units[ 0 ] == result ) {
1114 FreeTree( result, status );
1115 result = CopyTree( tree1, status );
1116
1117 /* Otherwise, search for the node to be replaced and do the substitution
1118 within the tree created earlier. */
1119 } else {
1120 ReplaceNode( result, units[ 0 ], CopyTree( tree1, status ), status );
1121 }
1122 }
1123 }
1124
1125 /* Free resources. */
1126 units = astFree( units );
1127
1128 /* If an error has occurred, free any new tree. */
1129 if( !astOK ) result = FreeTree( result, status );
1130
1131 /* Return the answer. */
1132 return result;
1133 }
1134
ConStart(const char * text,double * val,int * nc,int * status)1135 static int ConStart( const char *text, double *val, int *nc, int *status ) {
1136 /*
1137 * Name:
1138 * ConStart
1139
1140 * Purpose:
1141 * See if the supplied string starts with a literal numeric constant.
1142
1143 * Type:
1144 * Private function.
1145
1146 * Synopsis:
1147 * #include "unit.h"
1148 * int ConStart( const char *text, double *val, int *nc, int *status )
1149
1150 * Class Membership:
1151 * Unit member function.
1152
1153 * Description:
1154 * This function checks if the supplied string starts with a literal
1155 * numeric constant and returns it if it does. It is a wrap-up for scanf
1156 * since scanf has non-standard behaviour on some platforms (e.g. Cygwin
1157 * scanf interprets the character "n" as a floating point number!).
1158
1159 * Parameters:
1160 * text
1161 * The text to check.
1162 * val
1163 * Address of a double to receive any numerical constant read
1164 * from the start of the string. Unity is returned if the string
1165 * does not start with a numerical constant.
1166 * nc
1167 * Address of an int to receive the number of characters used to
1168 * create the value returned in "val". Zero is returned if the
1169 * string does not start with a numerical constant.
1170 * status
1171 * Pointer to the inherited status variable.
1172
1173 * Returned Value:
1174 * Non-zero if the text started with a numerical constant.
1175
1176 */
1177
1178 /* Local Variables: */
1179 int result;
1180 const char *c;
1181
1182 /* Initialise */
1183 *nc = 0;
1184 *val = 1.0;
1185
1186 /* Return zero if no text was supplied */
1187 if( !text ) return 0;
1188
1189 /* Use sscanf to see if the string begin with a numerical constant */
1190 result = astSscanf( text, "%lf%n", val, nc );
1191
1192 /* If so, check that the first non-blank character in the string
1193 is not "N" (interpreted by Cygwin as numerical zero!). */
1194 if( result ) {
1195 c = text;
1196 while( isspace( *c ) ) c++;
1197 if( *c == 'n' || *c == 'N' ) {
1198 result = 0;
1199 *nc = 0;
1200 *val = 1.0;
1201 }
1202 }
1203
1204 /* Return the result. */
1205 return result;
1206 }
1207
CopyTree(UnitNode * tree,int * status)1208 static UnitNode *CopyTree( UnitNode *tree, int *status ) {
1209 /*
1210 * Name:
1211 * CopyTree
1212
1213 * Purpose:
1214 * Create a new tree of UnitNodes containing a copy of a given tree.
1215
1216 * Type:
1217 * Private function.
1218
1219 * Synopsis:
1220 * #include "unit.h"
1221 * UnitNode *CopyTree( UnitNode *tree, int *status )
1222
1223 * Class Membership:
1224 * Unit member function.
1225
1226 * Description:
1227 * This function creates a copy of the supplied tree of UnitNodes.
1228
1229 * Parameters:
1230 * tree
1231 * The UnitNode at the head of the tree to be copied.
1232 * status
1233 * Pointer to the inherited status variable.
1234
1235 * Returned Value:
1236 * A pointer to the UnitNode at the head of the new tree.
1237
1238 * Notes:
1239 * - A value of NULL will be returned if this function is invoked with
1240 * the global error status set, or if it should fail for any reason.
1241 */
1242
1243 /* Local Variables: */
1244 UnitNode **args;
1245 UnitNode *result;
1246 int i;
1247 int narg;
1248
1249 /* Initialise. */
1250 result = NULL;
1251
1252 /* Check the inherited status. */
1253 if( !astOK || !tree ) return result;
1254
1255 /* Create a new node to represent the head of the supplied tree. */
1256 result = astMalloc( sizeof( UnitNode ) );
1257
1258 /* Check pointers can be used safely. */
1259 if( astOK ) {
1260
1261 /* Copy the fields of the supplied node. */
1262 narg = tree->narg;
1263
1264 result->arg = NULL;
1265 result->unit = tree->unit;
1266 result->mult = tree->mult;
1267 result->opcode = tree->opcode;
1268 result->narg = narg;
1269 result->con = tree->con;
1270 result->name = tree->name ? astStore( NULL, tree->name,
1271 strlen( tree->name ) + 1 ) : NULL;
1272
1273 /* Create an array of UnitNode pointers for the arguments. */
1274 args = astMalloc( narg*sizeof( UnitNode * ) );
1275 if( astOK ) {
1276 result->arg = args;
1277
1278 /* Copy the sub-trees headed by the argument nodes. */
1279 for( i = 0; i < narg; i++ ) {
1280 args[ i ] = CopyTree( tree->arg[ i ], status );
1281 }
1282 }
1283 }
1284
1285 /* Free any result if an error occurred. */
1286 if( !astOK ) result = FreeTree( result, status );
1287
1288 /* Return the answer. */
1289 return result;
1290 }
1291
CreateTree(const char * exp,int basic,int lock,int * status)1292 static UnitNode *CreateTree( const char *exp, int basic, int lock, int *status ){
1293 /*
1294 * Name:
1295 * CreateTree
1296
1297 * Purpose:
1298 * Convert an algebraic units expression into a tree of UnitNodes.
1299
1300 * Type:
1301 * Private function.
1302
1303 * Synopsis:
1304 * #include "unit.h"
1305 * UnitNode *CreateTree( const char *exp, int basic, int lock, int *status )
1306
1307 * Class Membership:
1308 * Unit member function.
1309
1310 * Description:
1311 * This function converts the supplied algebraic units expression into
1312 * a tree of UnitNodes. The result tree can optionally be expanded to
1313 * create a tree in which the "roots" (LDVAR nodes) all refer to
1314 * basic units.
1315
1316 * Parameters:
1317 * exp
1318 * The units expression. This should not include any leading or
1319 * trailing spaces.
1320 * basic
1321 * Should the tree created from parsing "exp" be expanded so that
1322 * the leaf nodes of the tree are all basic units?
1323 * lock
1324 * Use a mutex to guard access to the KnownUnits list?
1325 * status
1326 * Pointer to the inherited status variable.
1327
1328 * Returned Value:
1329 * A pointer to a UnitNode which forms the head of a tree of UnitNodes
1330 * representing the supplied unit expression.
1331
1332 * Notes:
1333 * - A NULL value is returned if this function is invoked with the
1334 * global error status set or if it should fail for any reason.
1335 */
1336
1337 /* Local Variables: */
1338 UnitNode *result;
1339 const char *cleanex;
1340
1341 /* Initialise */
1342 result = NULL;
1343
1344 /* Check the global error status, and that we have a string. */
1345 if ( !astOK ) return result;
1346
1347 /* Produce a clean copy of the supplied string. This has no leading
1348 or trailing white space, and any spaces adjacent to operators within
1349 the string are removed (this is needed because spaces are treated as
1350 multiplication symbols). */
1351 cleanex = CleanExp( exp, status );
1352
1353 /* If the string is blank, return the NULL pointer. Otherwise, create a
1354 tree of UnitNodes describing the units. The returned tree has LDVAR
1355 nodes which refer to the unit symbols contained in the supplied string. */
1356 if( cleanex && (*cleanex) ) {
1357 result = MakeTree( cleanex, strlen( cleanex ), lock, status );
1358
1359 /* Replace each subtree which simply combines constants (i.e. which has no
1360 OP_LDVAR nodes) with a single OP_LDCON node. */
1361 FixConstants( &result, 0, status );
1362
1363 /* Invert literal constant unit multipliers. */
1364 InvertConstants( &result, status );
1365
1366 /* Now replace each LDVAR node which refers to a known derived unit with
1367 a sub-tree which defines the derived unit in terms of known basic units.
1368 The LDVAR nodes in the resulting tree all refer to basic units. */
1369 if( basic ) RemakeTree( &result, status );
1370 }
1371
1372 /* Free resources. */
1373 cleanex = astFree( (void *) cleanex );
1374
1375 /* Free any returned tree if an error has occurred. */
1376 if( !astOK ) result = FreeTree( result, status );
1377
1378 /* Return the result. */
1379 return result;
1380 }
1381
DimAnal(UnitNode * node,double powers[NQUANT],double * scale,int * status)1382 static int DimAnal( UnitNode *node, double powers[NQUANT], double *scale, int *status ) {
1383 /*
1384 * Name:
1385 * DimAnal
1386
1387 * Purpose:
1388 * Perform a dimensional analysis of a unit tree.
1389
1390 * Type:
1391 * Protected function.
1392
1393 * Synopsis:
1394 * #include "unit.h"
1395 * int DimAnal( UnitNode *node, double powers[NQUANT], double *scale, int *status )
1396
1397 * Class Membership:
1398 * Unit member function.
1399
1400 * Description:
1401 * This function returns a set of powers and a scaling factor which
1402 * represent the units tree.
1403
1404 * Parameters:
1405 * node
1406 * Pointer to the UnitNode at the head of the unit tree.
1407 * powers
1408 * An array in which are returned the powers for each of the following
1409 * basic units (in the order shown): kilogramme, metre, second, radian,
1410 * Kelvin, count, adu, photon, magnitude, pixel. If the supplied unit
1411 * does not depend on a given basic unit a value of 0.0 will be returned
1412 * in the array. The returns values represent a system of units which is
1413 * a scaled form of the supplied units, expressed in the basic units of
1414 * m, kg, s, rad, K, count, adu, photon, mag and pixel. For instance, a
1415 * returned array of [1,0,-2,0,0,0,0,0,0,0] would represent "m/s**2".
1416 * scale
1417 * Pointer to a location at which to return a scaling factor for the
1418 * supplied units. The is the value, in the units represented by the
1419 * returned powers, which corresponds to a value of 1.0 in the supplied
1420 * units.
1421 * status
1422 * Pointer to the inherited status variable.
1423
1424 * Returned Value:
1425 * Non-zero if the tree was analysed succesfully. Zero otherwise.
1426
1427 * Notes:
1428 * - Zero is returned if this function is invoked with the
1429 * global error status set or if it should fail for any reason.
1430 */
1431
1432 /* Local Variables; */
1433 Oper oper;
1434 int result;
1435 int i;
1436 double p0[ NQUANT ];
1437 double p1[ NQUANT ];
1438 double s0;
1439 double s1;
1440
1441 /* Check inherited status */
1442 if( !astOK ) return 0;
1443
1444 /* Initialise the powers of all dimensions to zero, and set the scaling
1445 factor to unity. */
1446 result = 1;
1447 *scale = 1.0;
1448 for( i = 0; i < NQUANT; i++ ) powers[ i ] = 0.0;
1449
1450 /* Load constant: constant is dimensionaless so leave powers unchanged,
1451 and set the scaling factor. */
1452 oper = node->opcode;
1453 if( oper == OP_LDCON ) {
1454 *scale = 1.0/node->con;
1455
1456 /* Load variable: check it is one of the basic known dimensional
1457 quantities. If so, set the power of the quantity to unity and store
1458 the scale factor. If the unit is "g" modify the scale factor so that
1459 the analysis quantity is "kg". */
1460 } else if( oper == OP_LDVAR ) {
1461 result = 0;
1462 for( i = 0; i < NQUANT; i++ ) {
1463 if( node->unit == quant_units[ i ] ) {
1464 powers[ i ] = 1.0;
1465 *scale = node->mult ? 1.0/node->mult->scale : 1.0;
1466 if( !strcmp( node->unit->sym, "g" ) ) *scale *= 0.001;
1467 result = 1;
1468 break;
1469 }
1470 }
1471
1472 /* How does dimensional analysis handle log or exp units?*/
1473 } else if( oper == OP_LOG ) {
1474 result= 0;
1475
1476 } else if( oper == OP_LN ) {
1477 result= 0;
1478
1479 } else if( oper == OP_EXP ) {
1480 result= 0;
1481
1482 /* Get the powers for the child unit and then multiply each by 0.5 and
1483 take the square root of the scale factor. */
1484 } else if( oper == OP_SQRT ) {
1485 result = DimAnal( node->arg[0], powers, scale, status );
1486 if( result ) {
1487 for( i = 0; i < NQUANT; i++ ) powers[ i ]*= 0.5;
1488 *scale = sqrt( *scale );
1489 }
1490
1491 /* Similarly for pow nodes. */
1492 } else if( oper == OP_POW ) {
1493 result = DimAnal( node->arg[0], powers, scale, status );
1494 if( result ) {
1495 double power = node->arg[1]->con;
1496 for( i = 0; i < NQUANT; i++ ) powers[ i ]*= power;
1497 *scale = pow( *scale, power );
1498 }
1499
1500 /* Binary operators. Analyses the operands dimensions and combine. */
1501 } else if( oper == OP_DIV ) {
1502 if( DimAnal( node->arg[0], p0, &s0, status ) &&
1503 DimAnal( node->arg[1], p1, &s1, status ) ) {
1504 for( i = 0; i < NQUANT; i++ ) powers[ i ] = p0[ i ] - p1[ i ];
1505 *scale = s0/s1;
1506 } else {
1507 result = 0;
1508 }
1509
1510 } else if( oper == OP_MULT ) {
1511 if( DimAnal( node->arg[0], p0, &s0, status ) &&
1512 DimAnal( node->arg[1], p1, &s1, status ) ) {
1513 for( i = 0; i < NQUANT; i++ ) powers[ i ] = p0[ i ] + p1[ i ];
1514 *scale = s0*s1;
1515 } else {
1516 result = 0;
1517 }
1518
1519 /* Named constants are dimensionless */
1520 } else if( oper == OP_LDPI ) {
1521 *scale = 1.0/PI;
1522
1523 } else if( oper == OP_LDE ) {
1524 *scale = 1.0/E;
1525
1526 }
1527
1528 return result;
1529
1530 }
1531
EndsWith(const char * c,int nc,const char * test,int * status)1532 static int EndsWith( const char *c, int nc, const char *test, int *status ){
1533 /*
1534 * Name:
1535 * EndsWith
1536
1537 * Purpose:
1538 * See if a string ends with another string
1539
1540 * Type:
1541 * Private function.
1542
1543 * Synopsis:
1544 * #include "unit.h"
1545 * int EndsWith( const char *c, int nc, const char *test, int *status )
1546
1547 * Class Membership:
1548 * Unit member function.
1549
1550 * Description:
1551 * This function sees if the string given by "c" ends with the string
1552 * given by "test". The comparison is case-insensitive.
1553
1554 * Parameters:
1555 * c
1556 * A pointer to the last character in the string to be tested.
1557 * nc
1558 * The number of characters in the string to be tested.
1559 * test
1560 * A pointer to the string to be tested for.
1561 * status
1562 * Pointer to the inherited status variable.
1563
1564 * Returned Value:
1565 * Non-zero if the string "c" ends with the string "test".
1566
1567 */
1568
1569 /* Local Variables: */
1570 const char *start;
1571 int i;
1572 int result;
1573 int tlen;
1574
1575 /* initialise. */
1576 result = 0;
1577
1578 /* Check inherited status. */
1579 if( !astOK ) return result;
1580
1581 /* Check the string being tested for is not longer than the string being
1582 tested. */
1583 tlen = strlen( test );
1584 if( tlen <= nc ){
1585
1586 /* Get a pointer to where the matching string would start if the string "c"
1587 ends with the required string "test". */
1588 start = c - tlen + 1;
1589
1590 /* Do the comparison. */
1591 result = 1;
1592 for( i = 0; i < tlen; i++ ) {
1593 if( tolower( start[ i ] ) != tolower( test[ i ] ) ) {
1594 result = 0;
1595 break;
1596 }
1597 }
1598 }
1599
1600 /* Return the result. */
1601 return result;
1602
1603 }
1604
FindFactors(UnitNode * node,UnitNode *** factors,double ** powers,int * nfactor,double * coeff,int * status)1605 static void FindFactors( UnitNode *node, UnitNode ***factors, double **powers,
1606 int *nfactor, double *coeff, int *status ){
1607 /*
1608 * Name:
1609 * FindFactors
1610
1611 * Purpose:
1612 * Find the factors within an expression given by a tree of UnitNodes.
1613
1614 * Type:
1615 * Private function.
1616
1617 * Synopsis:
1618 * #include "unit.h"
1619 * void FindFactors( UnitNode *node, UnitNode ***factors, double **powers,
1620 * int *nfactor, double *coeff, int *status )
1621
1622 * Class Membership:
1623 * Unit member function.
1624
1625 * Description:
1626 * This function analyses the supplied tree of UnitNoes and returns
1627 * an array of pointers to nodes within the supplied tree which form
1628 * factors of the tree. The power associated with each factor is also
1629 * returned, together with an overall coefficient for the tree. The
1630 * expression represented by the tree is thus the product of the
1631 * coefficient with each of the factors, each raised to the associated
1632 * power.
1633
1634 * Parameters:
1635 * node
1636 * A pointer to the UnitNode at the head of the tree which is to be
1637 * analysed.
1638 * factors
1639 * The address at which to return a pointer to an array with "*nfactor"
1640 * elements, each element being a pointer to a UnitNode within the
1641 * supplied tree which is a factor of the supplied tree.
1642 * powers
1643 * The address at which to return a pointer to an array with "*nfactor"
1644 * elements, each element being a double holding the power of the
1645 * associated factor in "*factors".
1646 * nfactor
1647 * The address of an int containing the number of elements in the
1648 * returned "*factors" and "*powers" arrays.
1649 * coeff
1650 * The address of a double containing the overall coefficient to be
1651 * applied to the product of the factors.
1652 * status
1653 * Pointer to the inherited status variable.
1654
1655 * Notes:
1656 * - If the supplied node is a constant node, then "*coeff" is
1657 * returned holding the value of the constant and "*nfactor" is returned
1658 * equal to zero ("*factors" and "*powers" are returned holding NULL).
1659 * - If an error has already occurred, or if this function fails, then
1660 * "*factors" and "*powers" are returned holding NULL, "*nfactor" is
1661 * returned holding zero and "*coeff" is returned holding 1.0.
1662
1663 */
1664
1665 /* Local Variables: */
1666 int i;
1667 int j;
1668 int found;
1669 UnitNode **fact1;
1670 double *pow1;
1671 double coeff1;
1672 int nfac1;
1673 double con;
1674
1675 /* Initialise */
1676 *factors = NULL;
1677 *powers = NULL;
1678 *nfactor = 0;
1679 *coeff = 1.0;
1680
1681 /* Check inherited status. */
1682 if( !astOK ) return;
1683
1684 /* If the node at the head of the supplied tree is an OP_MULT node... */
1685 if( node->opcode == OP_MULT ) {
1686
1687 /* Find the factors of the two arguments of the OP_MULT node. */
1688 FindFactors( node->arg[ 0 ], factors, powers, nfactor, coeff, status );
1689 FindFactors( node->arg[ 1 ], &fact1, &pow1, &nfac1, &coeff1, status );
1690
1691 /* Combine the two lists. Loop round the factors of the seocnd argument. */
1692 for( i = 0; i < nfac1; i++ ) {
1693
1694 /* See if there is already an equivalent factor in the returned list of
1695 factors. */
1696 found = 0;
1697 for( j = 0; j < *nfactor; j++ ) {
1698 if( !CmpTree( (*factors)[ j ], fact1[ i ], 0, status ) ){
1699 found = 1;
1700 break;
1701 }
1702 }
1703
1704 /* If so, increment the power of the factor. */
1705 if( found ) {
1706 (*powers)[ j ] += pow1[ i ];
1707
1708 /* Otherwise, add the factor to the end of the returned list. */
1709 } else {
1710 *factors = astGrow( *factors, *nfactor + 1, sizeof( UnitNode *) );
1711 *powers = astGrow( *powers, *nfactor + 1, sizeof( double ) );
1712 if( astOK ) {
1713 (*factors)[ *nfactor ] = fact1[ i ];
1714 (*powers)[ (*nfactor)++ ] = pow1[ i ];
1715 }
1716 }
1717 }
1718
1719 /* Modify the overall coefficient. */
1720 *coeff *= coeff1;
1721
1722 /* Free resources */
1723 fact1 = astFree( fact1 );
1724 pow1 = astFree( pow1 );
1725
1726 /* If the node at the head of the supplied tree is an OP_POW node, */
1727 } else if( node->opcode == OP_POW ) {
1728
1729 /* Find the factors of the first argument. */
1730 FindFactors( node->arg[ 0 ], factors, powers, nfactor, coeff, status );
1731
1732 /* Multiply all the factor powers by the constant exponent of the POW
1733 node. */
1734 con = node->arg[ 1 ]->con;
1735 for( j = 0; j < *nfactor; j++ ) {
1736 (*powers)[ j ] *= con;
1737 }
1738
1739 /* Exponentiate the coefficient. */
1740 if( *coeff >= 0.0 || (int) con == con ) {
1741 *coeff = pow( *coeff, con );
1742 } else {
1743 astError( AST__BADUN, "Simplifying a units expression requires a "
1744 "negative value to be raised to a non-intergal power." , status);
1745 }
1746
1747 /* If the node at the head of the supplied tree is an OP_DIV node, */
1748 } else if( node->opcode == OP_DIV ) {
1749
1750 /* Find the factors of the two arguments of the OP_DIV node. */
1751 FindFactors( node->arg[ 0 ], factors, powers, nfactor, coeff, status );
1752 FindFactors( node->arg[ 1 ], &fact1, &pow1, &nfac1, &coeff1, status );
1753
1754 /* Combine the two lists. Loop round the factors of the second argument
1755 (the denominator). */
1756 for( i = 0; i < nfac1; i++ ) {
1757
1758 /* See if there is already an equivalent factor in the returned list of
1759 factors. */
1760 found = 0;
1761 for( j = 0; j < *nfactor; j++ ) {
1762 if( !CmpTree( (*factors)[ j ], fact1[ i ], 0, status ) ){
1763 found = 1;
1764 break;
1765 }
1766 }
1767
1768 /* If so, decrement the power of the factor. */
1769 if( found ) {
1770 (*powers)[ j ] -= pow1[ i ];
1771
1772 /* Otherwise, add the factor to the end of the returned list, with a
1773 negated power. */
1774 } else {
1775 *factors = astGrow( *factors, *nfactor + 1, sizeof( UnitNode *) );
1776 *powers = astGrow( *powers, *nfactor + 1, sizeof( double ) );
1777 if( astOK ) {
1778 (*factors)[ *nfactor ] = fact1[ i ];
1779 (*powers)[ (*nfactor)++ ] = -pow1[ i ];
1780 }
1781 }
1782 }
1783
1784 /* Modify the overall coefficient. */
1785 if( coeff1 != 0.0 ) {
1786 *coeff /= coeff1;
1787 } else {
1788 astError( AST__BADUN, "Simplifying a units expression"
1789 "requires a division by zero." , status);
1790 }
1791
1792 /* Free resources */
1793 fact1 = astFree( fact1 );
1794 pow1 = astFree( pow1 );
1795
1796 /* If the node at the head of the supplied tree is an OP_SQRT node, */
1797 } else if( node->opcode == OP_SQRT ) {
1798
1799 /* Find the factors of the argument. */
1800 FindFactors( node->arg[ 0 ], factors, powers, nfactor, coeff, status );
1801
1802 /* Multiply all the factor powers by 0.5. */
1803 for( j = 0; j < *nfactor; j++ ) {
1804 (*powers)[ j ] *= 0.5;
1805 }
1806
1807 /* Square root the coefficient. */
1808 if( *coeff >= 0.0 ) {
1809 *coeff = sqrt( *coeff );
1810 } else {
1811 astError( AST__BADUN, "Simplifying a units expression requires "
1812 "the square root of a negative value to be taken." , status);
1813 }
1814
1815 /* If the node at the head of the supplied tree is constant we have no
1816 factors but we have a coeffcient. */
1817 } else if( node->con != AST__BAD ) {
1818 *coeff = node->con;
1819
1820 /* Other nodes have no factors other than themselves, so just return a
1821 pointer to the supplied node. */
1822 } else {
1823 *factors = astMalloc( sizeof( UnitNode *) );
1824 *powers = astMalloc( sizeof( double ) );
1825 if( astOK ) {
1826 *nfactor = 1;
1827 (*factors)[ 0 ] = node;
1828 (*powers)[ 0 ] = 1.0;
1829 *coeff = 1.0;
1830 }
1831 }
1832
1833 /* If an error has occurred, free any returned resources. */
1834 if( !astOK ) {
1835 *factors = astFree( *factors );
1836 *powers = astFree( *powers );
1837 *nfactor = 0;
1838 *coeff = 1.0;
1839 }
1840 }
1841
FixConstants(UnitNode ** node,int unity,int * status)1842 static void FixConstants( UnitNode **node, int unity, int *status ) {
1843 /*
1844 * Name:
1845 * FixConstants
1846
1847 * Purpose:
1848 * Take the reciprocal of all constants in a tree of UnitNodes.
1849
1850 * Type:
1851 * Private function.
1852
1853 * Synopsis:
1854 * #include "unit.h"
1855 * void FixConstants( UnitNode **node, int unity, int *status )
1856
1857 * Class Membership:
1858 * Unit member function.
1859
1860 * Description:
1861 * This function replaces sub-trees which have a constant value by
1862 * a single OP_LDCON node which loads the appropriate constant.
1863
1864 * Parameters:
1865 * node
1866 * The address of a pointer to the UnitNode at the head of the tree
1867 * which is to be fixed. On exit the supplied tree is freed and a
1868 * pointer to a new tree is palced at he given address.
1869 * unity
1870 * If non-zero, then all multiplicative constants are set to 1.0, and
1871 * their original values are forgotten, but only if the other
1872 * argument of the OP_MULT node is an OP_LDVAR, OP_POW or OP_SQRT Node.
1873 * status
1874 * Pointer to the inherited status variable.
1875
1876 */
1877
1878 /* Local Variables: */
1879 int i;
1880 UnitNode *newnode;
1881 int allcon;
1882 Oper op;
1883 double newcon;
1884
1885 /* Check inherited status and pointer. */
1886 if( !astOK || !node || !(*node) ) return;
1887
1888 /* Initiallially, we have no replacement node */
1889 newnode = NULL;
1890 newcon = AST__BAD;
1891
1892 /* There is nothing to fix if the node has no arguments. */
1893 if( (*node)->narg > 0 ) {
1894
1895 /* Note the op code for the node. */
1896 op = (*node)->opcode;
1897
1898 /* Fix up the argument nodes. Also note if all the arguments are
1899 constants. */
1900 allcon = 1;
1901 for( i = 0; i < (*node)->narg; i++ ) {
1902 FixConstants( &( (*node)->arg[ i ] ), unity, status );
1903 if( (*node)->arg[ i ]->con == AST__BAD ) allcon = 0;
1904 }
1905
1906 /* If an OP_MULT nodes within a simplified tree has a constant argument,
1907 it will always be argument zero. If this is an OP_MULT node and arg[0]
1908 is constant and "unity" is non-zero and arg[1] is an OP_LDVAR, OP_POW
1909 or OP_SQRT node, replace the constant value by 1.0. */
1910 if( unity && op == OP_MULT &&
1911 (*node)->arg[ 0 ]->con != AST__BAD &&
1912 ( (*node)->arg[ 1 ]->opcode == OP_LDVAR ||
1913 (*node)->arg[ 1 ]->opcode == OP_SQRT ||
1914 (*node)->arg[ 1 ]->opcode == OP_POW ) ) {
1915 (*node)->arg[ 0 ]->con = 1.0;
1916 }
1917
1918 /* If the arguments of this node are all constants, replace the node by
1919 an OP_LDCON node which loads the resulting constant value. */
1920 if( allcon ) {
1921 if( (*node)->narg > 0 ) {
1922 newnode = NewNode( NULL, OP_LDCON, status );
1923 if( astOK ) {
1924 if( op == OP_LOG ) {
1925 if( (*node)->arg[ 0 ]->con > 0.0 ) {
1926 newcon = log10( (*node)->arg[ 0 ]->con );
1927 } else {
1928 astError( AST__BADUN, "Illegal negative or zero constant "
1929 "value '%g' encountered.", status,
1930 (*node)->arg[ 0 ]->con );
1931 }
1932 } else if( op == OP_LN ){
1933 if( (*node)->arg[ 0 ]->con > 0.0 ) {
1934 newcon = log( (*node)->arg[ 0 ]->con );
1935 } else {
1936 astError( AST__BADUN, "Illegal negative or zero constant value "
1937 "'%g' encountered.", status, (*node)->arg[ 0 ]->con );
1938 }
1939 } else if( op == OP_EXP ){
1940 newcon = exp( (*node)->arg[ 0 ]->con );
1941
1942 } else if( op == OP_SQRT ){
1943 if( (*node)->arg[ 0 ]->con >= 0.0 ) {
1944 newcon = sqrt( (*node)->arg[ 0 ]->con );
1945 } else {
1946 astError( AST__BADUN, "Illegal negative constant value "
1947 "'%g' encountered.", status, (*node)->arg[ 0 ]->con );
1948 }
1949
1950 } else if( op == OP_POW ){
1951 if( (*node)->arg[ 0 ]->con >= 0.0 ||
1952 (int) (*node)->arg[ 1 ]->con == (*node)->arg[ 1 ]->con ) {
1953 newcon = pow( (*node)->arg[ 0 ]->con,
1954 (*node)->arg[ 1 ]->con );
1955 } else {
1956 astError( AST__BADUN, "Illegal negative constant value "
1957 "'%g' encountered.", status, (*node)->arg[ 0 ]->con );
1958 }
1959
1960 } else if( op == OP_DIV ){
1961 if( (*node)->arg[ 1 ]->con != 0.0 ) {
1962 newcon = (*node)->arg[ 0 ]->con / (*node)->arg[ 1 ]->con;
1963 } else {
1964 astError( AST__BADUN, "Illegal zero constant value encountered." , status);
1965 }
1966
1967 } else if( op == OP_MULT ){
1968 newcon = (*node)->arg[ 0 ]->con * (*node)->arg[ 1 ]->con;
1969
1970 }
1971
1972
1973 if( astOK ) newnode->con = newcon;
1974 }
1975 }
1976 }
1977 }
1978
1979 /* If an error has occurred, free any new node. */
1980 if( !astOK ) newnode = FreeTree( newnode, status );
1981
1982 /* If we have a replacement node, free the supplied tree and return a
1983 pointer to the new tree. */
1984 if( newnode ) {
1985 FreeTree( *node, status );
1986 *node = newnode;
1987 }
1988
1989 }
1990
FixUnits(UnitNode * node,UnitNode * test,int * status)1991 static UnitNode *FixUnits( UnitNode *node, UnitNode *test, int *status ) {
1992 /*
1993 * Name:
1994 * FixUnits
1995
1996 * Purpose:
1997 * Assign a constant value to all units except for one.
1998
1999 * Type:
2000 * Private function.
2001
2002 * Synopsis:
2003 * #include "unit.h"
2004 * UnitNode *FixUnits( UnitNode *node, UnitNode *test, int *status )
2005
2006 * Class Membership:
2007 * Unit member function.
2008
2009 * Description:
2010 * This function returns a copy of the supplied tree of UnitNodes. All
2011 * OP_LDVAR nodes within the copy which refer to units which differ
2012 * from those referred to by the supplied test node are replaced by
2013 * OP_LDCON nodes which load the constant value 1.0.
2014
2015 * Parameters:
2016 * node
2017 * A pointer to the UnitNode at the head of the tree to be used.
2018 * test
2019 * A pointer to an OP_LDVAR node which defines the units which are
2020 * *not* to be replaced by a constant value of 1.0.
2021 * status
2022 * Pointer to the inherited status variable.
2023
2024 * Returned Value:
2025 * A pointer to a UnitNode which is at the head of a tree of UnitNodes
2026 * which forms the required copy of th einput tree.
2027
2028 * Notes:
2029 * - A NULL pointer is returned if an error has already occurred, or
2030 * if this function fails for any reason.
2031
2032 */
2033
2034 /* Local Variables: */
2035 int i;
2036 UnitNode *result;
2037
2038 /* Initialise. */
2039 result = NULL;
2040
2041 /* Check inherited status. */
2042 if( !astOK ) return result;
2043
2044 /* Create a complete copy of the supplied tree. */
2045 result = CopyTree( node, status );
2046
2047 /* Is the node at the head of the supplied tree an OP_LDVAR node? */
2048 if( node->opcode == OP_LDVAR ) {
2049
2050 /* Does it refer to a unit which differs from that of the test node? If so
2051 annul the copy created above and return a new OP_LDCON node which loads
2052 the constant value 1.0. */
2053 if( strcmp( test->name, node->name ) ) {
2054 FreeTree( result, status );
2055 result = NewNode( NULL, OP_LDCON, status );
2056 if( astOK ) result->con = 1.0;
2057 }
2058
2059 /* If the supplied node is not an OP_LDVAR node, check each argument of
2060 the head node. */
2061 } else {
2062 for( i = 0; i < node->narg; i++ ) {
2063
2064 /* Free the resources used to hold this argument in the tree copy created
2065 above. */
2066 FreeTree( result->arg[ i ], status );
2067
2068 /* Create a new argument tree by calling this function recursively to
2069 fix units in the argument sub-trees. */
2070 result->arg[ i ] = FixUnits( node->arg[ i ], test, status );
2071 }
2072 }
2073
2074 /* If an error has occurred, free any new tree. */
2075 if( !astOK ) result = FreeTree( result, status );
2076
2077 /* Return the answer. */
2078 return result;
2079 }
2080
FreeTree(UnitNode * node,int * status)2081 static UnitNode *FreeTree( UnitNode *node, int *status ) {
2082 /*
2083 * Name:
2084 * FreeTree
2085
2086 * Purpose:
2087 * Free resources used by a tree of UnitNodes.
2088
2089 * Type:
2090 * Private function.
2091
2092 * Synopsis:
2093 * #include "unit.h"
2094 * UnitNode *FreeTree( UnitNode *node, int *status )
2095
2096 * Class Membership:
2097 * Unit member function.
2098
2099 * Description:
2100 * This function frees the memory used to store a tree of UnitNodes.
2101
2102 * Parameters:
2103 * node
2104 * A pointer to the UnitNode at the head of the tree which is to be
2105 * freed.
2106 * status
2107 * Pointer to the inherited status variable.
2108
2109 * Returned Value:
2110 * A NULL pointer is returned.
2111
2112 * Notes:
2113 * - This function attempts to execute even if it is invoked with
2114 * the global error status set.
2115 */
2116
2117 /* Local Variables: */
2118 int i;
2119
2120 /* Check a node was supplied. */
2121 if( node ) {
2122
2123 /* Recursively free any argument nodes. */
2124 if( node->arg ) {
2125 for( i = 0; i < node->narg; i++ ) {
2126 (node->arg)[ i ] = FreeTree( (node->arg)[ i ], status );
2127 }
2128 node->arg = astFree( node->arg );
2129 }
2130
2131 /* Nullify other pointers for safety. */
2132 node->unit = NULL;
2133 node->mult = NULL;
2134
2135 /* Free the copy of the symbol string (if any). */
2136 node->name = astFree( (char *) node->name );
2137
2138 /* Free the memory holding the node. */
2139 node = astFree( node );
2140 }
2141
2142 /* Return a null pointer. */
2143 return NULL;
2144 }
2145
GetKnownUnits(int lock,int * status)2146 static KnownUnit *GetKnownUnits( int lock, int *status ) {
2147 /*
2148 * Name:
2149 * GetKnownUnits
2150
2151 * Purpose:
2152 * Get a pointer to the head of a linked list of known unit definitions.
2153
2154 * Type:
2155 * Private function.
2156
2157 * Synopsis:
2158 * #include "unit.h"
2159 * KnownUnit *GetKnownUnits( int lock, int *status )
2160
2161 * Class Membership:
2162 * Unit member function.
2163
2164 * Description:
2165 * This function returns a pointer to the head of a linked list of known
2166 * unit definitions. The unit definitions are created as static module
2167 * variables if they have not previously been created.
2168
2169 * Parameters:
2170 * lock
2171 * If non-zero, then lock a mutex prior to accessing the list of
2172 * known units.
2173 * status
2174 * Pointer to the inherited status variable.
2175
2176 * Returned Value:
2177 * A pointer to the first known unit definition.
2178
2179 * Notes:
2180 * - A NULL pointer is returned if it is invoked with the global error
2181 * status set, or if an error occurs.
2182 */
2183
2184 /* Local Variables: */
2185 int iq;
2186 KnownUnit *result;
2187
2188 /* Initialise. */
2189 result = NULL;
2190
2191 /* Check inherited status. */
2192 if( !astOK ) return result;
2193
2194 /* Ensure the known units list is only initialised once. */
2195 if( lock ) {
2196 LOCK_MUTEX1
2197 }
2198
2199 /* If the linked list of KnownUnit structures describing the known units
2200 has not yet been created, create it now. A pointer to the head of the
2201 linked list is put into the static variable "known_units". */
2202 if( !known_units ) {
2203
2204 /* At the same time we store pointers to the units describing the basic
2205 quantities used in dimensional analysis. Initialise th index of the
2206 next such unit. */
2207 iq = 0;
2208
2209 /* Create definitions for the known units. First do all IAU basic units.
2210 We include "g" instead of "kg" because otherwise we would have to
2211 refer to a gramme as a milli-kilogramme. */
2212 MakeKnownUnit( "g", "gram", NULL, status );
2213 quant_units[ iq++ ] = known_units;
2214 MakeKnownUnit( "m", "metre", NULL, status );
2215 quant_units[ iq++ ] = known_units;
2216 MakeKnownUnit( "s", "second", NULL, status );
2217 quant_units[ iq++ ] = known_units;
2218 MakeKnownUnit( "rad", "radian", NULL, status );
2219 quant_units[ iq++ ] = known_units;
2220 MakeKnownUnit( "K", "Kelvin", NULL, status );
2221 quant_units[ iq++ ] = known_units;
2222 MakeKnownUnit( "A", "Ampere", NULL, status );
2223 MakeKnownUnit( "mol", "mole", NULL, status );
2224 MakeKnownUnit( "cd", "candela", NULL, status );
2225
2226 /* Now do all IAU derived units. Unit definitions may only refer to units
2227 which have already been defined. */
2228 MakeKnownUnit( "sr", "steradian", "rad rad", status );
2229 MakeKnownUnit( "Hz", "Hertz", "1/s", status );
2230 MakeKnownUnit( "N", "Newton", "kg m/s**2", status );
2231 MakeKnownUnit( "J", "Joule", "N m", status );
2232 MakeKnownUnit( "W", "Watt", "J/s", status );
2233 MakeKnownUnit( "C", "Coulomb", "A s", status );
2234 MakeKnownUnit( "V", "Volt", "J/C", status );
2235 MakeKnownUnit( "Pa", "Pascal", "N/m**2", status );
2236 MakeKnownUnit( "Ohm", "Ohm", "V/A", status );
2237 MakeKnownUnit( "S", "Siemens", "A/V", status );
2238 MakeKnownUnit( "F", "Farad", "C/V", status );
2239 MakeKnownUnit( "Wb", "Weber", "V s", status );
2240 MakeKnownUnit( "T", "Tesla", "Wb/m**2", status );
2241 MakeKnownUnit( "H", "Henry", "Wb/A", status );
2242 MakeKnownUnit( "lm", "lumen", "cd sr", status );
2243 MakeKnownUnit( "lx", "lux", "lm/m**2", status );
2244
2245 /* Now do additional derived and basic units listed in the FITS-WCS paper. */
2246 MakeKnownUnit( "deg", "degree", "pi/180 rad", status );
2247 MakeKnownUnit( "arcmin", "arc-minute", "1/60 deg", status );
2248 MakeKnownUnit( "arcsec", "arc-second", "1/3600 deg", status );
2249 MakeKnownUnit( "mas", "milli-arcsecond", "1/3600000 deg", status );
2250 MakeKnownUnit( "min", "minute", "60 s", status );
2251 MakeKnownUnit( "h", "hour", "3600 s", status );
2252 MakeKnownUnit( "d", "day", "86400 s", status );
2253 MakeKnownUnit( "yr", "year", "31557600 s", status );
2254 MakeKnownUnit( "a", "year", "31557600 s", status );
2255 MakeKnownUnit( "eV", "electron-Volt", "1.60217733E-19 J", status );
2256 MakeKnownUnit( "erg", "erg", "1.0E-7 J", status );
2257 MakeKnownUnit( "Ry", "Rydberg", "13.605692 eV", status );
2258 MakeKnownUnit( "solMass", "solar mass", "1.9891E30 kg", status );
2259 MakeKnownUnit( "u", "unified atomic mass unit", "1.6605387E-27 kg", status );
2260 MakeKnownUnit( "solLum", "solar luminosity", "3.8268E26 W", status );
2261 MakeKnownUnit( "Angstrom", "Angstrom", "1.0E-10 m", status );
2262 MakeKnownUnit( "micron", "micron", "1.0E-6 m", status );
2263 MakeKnownUnit( "solRad", "solar radius", "6.9599E8 m", status );
2264 MakeKnownUnit( "AU", "astronomical unit", "1.49598E11 m", status );
2265 MakeKnownUnit( "lyr", "light year", "9.460730E15 m", status );
2266 MakeKnownUnit( "pc", "parsec", "3.0867E16 m", status );
2267 MakeKnownUnit( "count", "count", NULL, status );
2268 quant_units[ iq++ ] = known_units;
2269 MakeKnownUnit( "adu", "analogue-to-digital unit", NULL, status );
2270 quant_units[ iq++ ] = known_units;
2271 MakeKnownUnit( "photon", "photon", NULL, status );
2272 quant_units[ iq++ ] = known_units;
2273 MakeKnownUnit( "Jy", "Jansky", "1.0E-26 W /m**2 /Hz", status );
2274 MakeKnownUnit( "mag", "magnitude", NULL, status );
2275 quant_units[ iq++ ] = known_units;
2276 MakeKnownUnit( "G", "Gauss", "1.0E-4 T", status );
2277 MakeKnownUnit( "pixel", "pixel", NULL, status );
2278 quant_units[ iq++ ] = known_units;
2279 MakeKnownUnit( "barn", "barn", "1.0E-28 m**2", status );
2280 MakeKnownUnit( "D", "Debye", "(1.0E-29/3) C.m", status );
2281
2282 if( iq != NQUANT && astOK ) {
2283 astError( AST__INTER, "unit(GetKnownUnits): %d basic quantities "
2284 "noted but this should be %d (internal AST programming "
2285 "error).", status, iq, NQUANT );
2286 }
2287
2288 /* Unit aliases... */
2289 MakeUnitAlias( "Angstrom", "Ang", status );
2290 MakeUnitAlias( "count", "ct", status );
2291 MakeUnitAlias( "photon", "ph", status );
2292 MakeUnitAlias( "Jy", "Jan", status );
2293 MakeUnitAlias( "pixel", "pix", status );
2294 MakeUnitAlias( "s", "sec", status );
2295 MakeUnitAlias( "m", "meter", status );
2296 }
2297
2298 /* If succesful, return the pointer to the head of the list. */
2299 if( astOK ) result = known_units;
2300
2301 /* Allow the next thread to proceed. */
2302 if( lock ) {
2303 UNLOCK_MUTEX1
2304 }
2305
2306 /* Return the result. */
2307 return result;
2308 }
2309
GetMultipliers(int * status)2310 static Multiplier *GetMultipliers( int *status ) {
2311 /*
2312 * Name:
2313 * GetMultiplier
2314
2315 * Purpose:
2316 * Get a pointer to the head of a linked list of multiplier definitions.
2317
2318 * Type:
2319 * Private function.
2320
2321 * Synopsis:
2322 * #include "unit.h"
2323 * Multiplier *Multipliers( void )
2324
2325 * Class Membership:
2326 * Unit member function.
2327
2328 * Description:
2329 * This function returns a pointer to the head of a linked list of known
2330 * multiplier definitions. The multiplier definitions are created as
2331 * static module variables if they have not previously been created.
2332
2333 * Returned Value:
2334 * A pointer to the first known multiplier definition.
2335
2336 * Notes:
2337 * - A NULL pointer is returned if it is invoked with the global error
2338 * status set, or if an error occurs.
2339 */
2340
2341 /* Local Variables: */
2342 Multiplier *result;
2343 Multiplier *mult;
2344
2345 /* Initialise. */
2346 result = NULL;
2347
2348 /* Check inherited status. */
2349 if( !astOK ) return result;
2350
2351 /* Ensure the list is only initialised by one thread. */
2352 LOCK_MUTEX2
2353
2354 /* If the linked list of Multiplier structures describing the known
2355 multipliers has not yet been created, create it now. A pointer to the
2356 head of the linked list is put into the static variable "multipliers". */
2357 if( !multipliers ) {
2358
2359 /* Define a macro to create a multiplier struncture and add it to the
2360 linked list of multiplier structures. */
2361 #define MAKEMULT(s,sl,sc,lab,ll) \
2362 mult = astMalloc( sizeof( Multiplier ) ); \
2363 if( astOK ) { \
2364 mult->sym = s; \
2365 mult->symlen = sl; \
2366 mult->lablen = ll; \
2367 mult->scale = sc; \
2368 mult->label = lab; \
2369 mult->next = multipliers; \
2370 multipliers = mult; \
2371 }
2372
2373 /* Use the above macro to create all the standard multipliers listed in the
2374 FITS WCS paper I. */
2375 MAKEMULT("d",1,1.0E-1,"deci",4)
2376 MAKEMULT("c",1,1.0E-2,"centi",5)
2377 MAKEMULT("m",1,1.0E-3,"milli",5)
2378 MAKEMULT("u",1,1.0E-6,"micro",5)
2379 MAKEMULT("n",1,1.0E-9,"nano",4)
2380 MAKEMULT("p",1,1.0E-12,"pico",4)
2381 MAKEMULT("f",1,1.0E-15,"femto",5)
2382 MAKEMULT("a",1,1.0E-18,"atto",4)
2383 MAKEMULT("z",1,1.0E-21,"zepto",5)
2384 MAKEMULT("y",1,1.0E-24,"yocto",5)
2385 MAKEMULT("da",2,1.0E1,"deca",4)
2386 MAKEMULT("h",1,1.0E2,"hecto",5)
2387 MAKEMULT("k",1,1.0E3,"kilo",4)
2388 MAKEMULT("M",1,1.0E6,"mega",4)
2389 MAKEMULT("G",1,1.0E9,"giga",4)
2390 MAKEMULT("T",1,1.0E12,"tera",4)
2391 MAKEMULT("P",1,1.0E15,"peta",4)
2392 MAKEMULT("E",1,1.0E18,"exa",3)
2393 MAKEMULT("Z",1,1.0E21,"zetta",5)
2394 MAKEMULT("Y",1,1.0E24,"yotta",5)
2395
2396 /* Undefine the macro. */
2397 #undef MAKEMULT
2398
2399 }
2400
2401 /* If succesful, return the pointer to the head of the list. */
2402 if( astOK ) result = multipliers;
2403
2404 /* Allow the next thread to proceed. */
2405 UNLOCK_MUTEX2
2406
2407 /* Return the result. */
2408 return result;
2409 }
2410
InvertConstants(UnitNode ** node,int * status)2411 static void InvertConstants( UnitNode **node, int *status ) {
2412 /*
2413 * Name:
2414 * InvertConstants
2415
2416 * Purpose:
2417 * Take the reciprocal of all constants in a tree of UnitNodes.
2418
2419 * Type:
2420 * Private function.
2421
2422 * Synopsis:
2423 * #include "unit.h"
2424 * void InvertConstants( UnitNode **node, int *status )
2425
2426 * Class Membership:
2427 * Unit member function.
2428
2429 * Description:
2430 * This function replaces constant unit coefficients by their reciprocal.
2431 * This is because a string such as "0.01 m" will be interpreted as
2432 * meaning "multiply a value in metres by 0.01 to get the value in the
2433 * required units", whereas what is actually meant is "use units of
2434 * 0.01 of a metre" which requires us to divide the value in metres by
2435 * 0.01, not multiply it.
2436
2437 * Parameters:
2438 * node
2439 * The address of a pointer to the UnitNode at the head of the tree.
2440 * On exit the supplied tree is freed and a pointer to a new tree is
2441 * placed at the given address.
2442 * status
2443 * Pointer to the inherited status variable.
2444
2445 */
2446
2447 /* Local Variables: */
2448 int i;
2449 UnitNode *newnode;
2450 int allcon;
2451 Oper op;
2452
2453 /* Check inherited status and pointer. */
2454 if( !astOK || !node || !(*node) ) return;
2455
2456 /* Initiallially, we have no replacement node */
2457 newnode = NULL;
2458
2459 /* There is nothing to fix if the node has no arguments. */
2460 if( (*node)->narg > 0 ) {
2461
2462 /* Note the op code for the node. */
2463 op = (*node)->opcode;
2464
2465 /* Fix up the argument nodes. Also note if all the arguments are
2466 constants. */
2467 allcon = 1;
2468 for( i = 0; i < (*node)->narg; i++ ) {
2469 InvertConstants( &( (*node)->arg[ i ] ), status );
2470 if( (*node)->arg[ i ]->con == AST__BAD ) allcon = 0;
2471 }
2472
2473 /* If all nodes are constant, there are no co-efficients to invert. */
2474 if( !allcon ) {
2475
2476 /* Iif this is a multiplication node, see if either of its arguments
2477 is a constant. If so, invert the constant. This is because a string like
2478 "0.01 m" means "each unit is 0.01 of a metre". Therefore, to transform
2479 a value in metres into required units means multiplying the metres
2480 value by 100.0 (i.e the reciprocal of 0.01), not 0.01. */
2481 if( op == OP_MULT ) {
2482 for( i = 0; i < 2; i++ ) {
2483 if( (*node)->arg[ i ]->con != AST__BAD ) {
2484 if( (*node)->arg[ i ]->con != 0.0 ) {
2485
2486 (*node)->arg[ i ]->con = 1.0/(*node)->arg[ i ]->con;
2487 } else {
2488 astError( AST__BADUN, "Illegal zero constant encountered." , status);
2489 }
2490 }
2491 }
2492
2493 /* Likewise, check for division nodes in which the denominator is
2494 constant. */
2495 } else if( op == OP_DIV ) {
2496 if( (*node)->arg[ 1 ]->con != AST__BAD ) {
2497 if( (*node)->arg[ 1 ]->con != 0.0 ) {
2498 (*node)->arg[ 1 ]->con = 1.0/(*node)->arg[ 1 ]->con;
2499 } else {
2500 astError( AST__BADUN, "Illegal zero constant encountered." , status);
2501 }
2502 }
2503
2504 /* If this is a "pow" node check that the second argument is constant
2505 (required by FITS WCS paper I). */
2506 } else if( op == OP_POW ) {
2507 if( (*node)->arg[ 1 ]->con == AST__BAD ) {
2508 astError( AST__BADUN, "Illegal variable exponent." , status);
2509 }
2510 }
2511 }
2512 }
2513
2514 /* If an error has occurred, free any new node. */
2515 if( !astOK ) newnode = FreeTree( newnode, status );
2516
2517 /* If we have a replacement node, free the supplied tree and return a
2518 pointer to the new tree. */
2519 if( newnode ) {
2520 FreeTree( *node, status );
2521 *node = newnode;
2522 }
2523 }
2524
InvertTree(UnitNode * fwdnode,UnitNode * src,int * status)2525 static UnitNode *InvertTree( UnitNode *fwdnode, UnitNode *src, int *status ) {
2526 /*
2527 * Name:
2528 * InvertTree
2529
2530 * Purpose:
2531 * Invert a tree of UnitNodes.
2532
2533 * Type:
2534 * Private function.
2535
2536 * Synopsis:
2537 * #include "unit.h"
2538 * UnitNode *InvertTree( UnitNode *fwdnode, UnitNode *src )
2539
2540 * Class Membership:
2541 * Unit member function.
2542
2543 * Description:
2544 * This function inverts a tree of UnitNodes. The supplied tree should
2545 * have exactly one OP_LDVAR node. This will be the quantity represented
2546 * by the node at the head of the returned tree.
2547
2548 * Parameters:
2549 * fwdnode
2550 * A pointer to the UnitNode at the head of the tree which is to be
2551 * inverted.
2552 * src
2553 * A pointer to a UnitNode which is to be used as the root of the
2554 * inverted tree. That is, the output from this node should form
2555 * the (one and only) varying input to the inverted tree. If the
2556 * supplied tree is succesfulyl inverted, the tree of which "src"
2557 * is the head will be contained within the returned inverted tree.
2558 * Therefore "src" only needs to be freed explicitly if this
2559 * function fails to invert the supplied tree for any reason. If
2560 * this function succeeds, then "src" will be freed as part of
2561 * freeing the returned inverted tree.
2562
2563 * Returned Value:
2564 * A pointer to a UnitNode which forms the head of the inverted tree.
2565
2566 * Algorithm:
2567 * The algorithm works through the supplied forward tree, from the head
2568 * to the roots. First, the supplied node at the head of the forward
2569 * tree is inverted. To be invertable, the supplied head node must have
2570 * exactly one varying argument (any other arguments must be fixed,
2571 * i.e. not vary). This varying argument becomes the output of the
2572 * inverted node. The other (fixed) arguments to the forward node are
2573 * also used as arguments to the inverted node. The supplied "src" node
2574 * is used as the single varying input to the inverted node. Having
2575 * inverted the supplied forward head node, this function is called
2576 * recursively to invert the lower parts of the forward tree (i.e. the
2577 * part of the forward tree which provided the varying input to node
2578 * which has just been inverted).
2579
2580 * Notes:
2581 * - It is assumed that he supplied forward tree has been simplified
2582 * using SimplifyTree. This means that the tree contains no nodes with
2583 * the following op codes: OP_LOG, OP_SQRT. OP_DIV (SimplifyTree
2584 * converts these nodes into OP_LN, OP_POW and OP_MULT nodes).
2585 * - A value of NULL will be returned if this function is invoked with
2586 * the global error status set, or if it should fail for any reason.
2587
2588 */
2589
2590 /* Local Variables: */
2591 UnitNode *newnode;
2592 UnitNode *nextnode;
2593 UnitNode *result;
2594 UnitNode *node1;
2595 Oper fop;
2596 int varg;
2597
2598 /* Initialise */
2599 result = NULL;
2600
2601 /* Check inherited status. */
2602 if( !astOK ) return result;
2603
2604 /* Initiallially, we have no replacement node */
2605 newnode = NULL;
2606 nextnode = NULL;
2607
2608 /* Save the op code at the head of the forward tree. */
2609 fop = fwdnode->opcode;
2610
2611 /* If the head of the forward tree is a OP_EXP node. Inverse of
2612 "exp(x)" is "ln(x)". */
2613 if( fop == OP_EXP ) {
2614 newnode = NewNode( NULL, OP_LN, status );
2615 if( astOK ) {
2616 newnode->arg[ 0 ] = src;
2617 nextnode = fwdnode->arg[ 0 ];
2618 }
2619
2620 /* If the head of the forward tree is a OP_LN node. Inverse of
2621 "ln(x)" is "exp(x)". */
2622 } else if( fop == OP_LN ) {
2623 newnode = NewNode( NULL, OP_EXP, status );
2624 if( astOK ) {
2625 newnode->arg[ 0 ] = src;
2626 nextnode = fwdnode->arg[ 0 ];
2627 }
2628
2629 /* If the head of the forward tree is a OP_POW node. Inverse of
2630 "x**k" is "x**(1/k)" */
2631 } else if( fop == OP_POW ) {
2632 newnode = NewNode( NULL, OP_POW, status );
2633 node1 = NewNode( NULL, OP_LDCON, status );
2634 if( astOK ) {
2635 node1->con = 1.0/fwdnode->arg[ 1 ]->con;
2636 newnode->arg[ 0 ] = src;
2637 newnode->arg[ 1 ] = node1;
2638 nextnode = fwdnode->arg[ 0 ];
2639 }
2640
2641 /* If the head of the forward tree is a OP_MULT node... */
2642 } else if( fop == OP_MULT ) {
2643
2644 /* The node is only invertable if it has one constant node and one
2645 non-constant node. Get the index of the varying argument. */
2646 if( fwdnode->arg[ 0 ]->con != AST__BAD &&
2647 fwdnode->arg[ 1 ]->con == AST__BAD ) {
2648 varg = 1;
2649 } else if( fwdnode->arg[ 0 ]->con == AST__BAD &&
2650 fwdnode->arg[ 1 ]->con != AST__BAD ) {
2651 varg = 0;
2652 } else {
2653 varg = -1;
2654 }
2655 if( varg != -1 ) {
2656
2657 /* The inverse of "k*x" is "(1/k)*x" (we use MULT nodes instead of DIV
2658 nodes to maintain the standardisation implemented by SimplifyTree). */
2659 newnode = NewNode( NULL, OP_MULT, status );
2660 node1 = NewNode( NULL, OP_LDCON, status );
2661 if( astOK ) {
2662 node1->con = 1.0/fwdnode->arg[ 1 - varg ]->con;
2663 newnode->arg[ 0 ] = node1;
2664 newnode->arg[ 1 ] = src;
2665 nextnode = fwdnode->arg[ varg ];
2666 }
2667 }
2668
2669 /* If the head of the forward tree is a OP_LDVAR node, there is nothing
2670 left to invert. SO return a pointer to the suppleid source node. */
2671 } else if( fop == OP_LDVAR ) {
2672 result = src;
2673 nextnode = NULL;
2674
2675 /* If the head of the forward tree is any other node (e.g. a OP_LDCON node),
2676 the tree cannot be inverted. */
2677 } else {
2678 nextnode = NULL;
2679 }
2680
2681 /* If we managed to invert the node at the head of the supplied tree,
2682 continue to invert its varying argument node (if any). */
2683 if( nextnode && newnode ) result = InvertTree( nextnode, newnode, status );
2684
2685 /* If the tree could not be inverted, free the newnode. */
2686 if( !result ) newnode = FreeTree( newnode, status );
2687
2688 /* If an error has occurred, free any new node. */
2689 if( !astOK ) result = FreeTree( result, status );
2690
2691 /* Return the result. */
2692 return result;
2693
2694 }
2695
LocateUnits(UnitNode * node,UnitNode *** units,int * nunits,int * status)2696 static void LocateUnits( UnitNode *node, UnitNode ***units, int *nunits, int *status ){
2697 /*
2698 * Name:
2699 * LocateUnits
2700
2701 * Purpose:
2702 * Locate the units used by a supplied tree of UnitNodes.
2703
2704 * Type:
2705 * Private function.
2706
2707 * Synopsis:
2708 * #include "unit.h"
2709 * void LocateUnits( UnitNode *node, UnitNode ***units, int *nunits, int *status )
2710
2711 * Class Membership:
2712 * Unit member function.
2713
2714 * Description:
2715 * This function locates the units used by a supplied tree of
2716 * UnitNodes.
2717
2718 * Parameters:
2719 * node
2720 * A pointer to the UnitNode at the head of the tree to be searched.
2721 * units
2722 * The address at which is stored a pointer to an array of "*nunits"
2723 * elements. Each element of the array holds a pointer to a UnitNode.
2724 * The array is extended on exit to hold pointers to the OP_LDVAR nodes
2725 * within the supplied tree (i.e. nodes which represent named units,
2726 * either known or unknown). A node is only included in the returned
2727 * array if no other node for the same unit is already included in the
2728 * array. A NULL pointer should be supplied on the first invocation of
2729 * this function.
2730 * nunits
2731 * The address of an integer which holds the number of elements in
2732 * the array given by "*units". Updated on exit to included any
2733 * elements added to the array. Zero should be supplied on the first
2734 * invocation of this function.
2735 * status
2736 * Pointer to the inherited status variable.
2737
2738 */
2739
2740 /* Local Variables: */
2741 int i;
2742 int found;
2743
2744 /* Check the global error status. */
2745 if( !astOK ) return;
2746
2747 /* Is the node at the head of the supplied tree an OP_LDVAR node? */
2748 if( node->opcode == OP_LDVAR ) {
2749
2750 /* If an array was supplied, see if it already contains a pointer to a node
2751 which refers to the same units. */
2752 found = 0;
2753 if( *units ) {
2754 for( i = 0; i < *nunits; i++ ) {
2755 if( !strcmp( (*units)[ i ]->name, node->name ) ) {
2756 found = 1;
2757 break;
2758 }
2759 }
2760 }
2761
2762 /* If not, ensure the array is big enough and add a pointer to the
2763 supplied node to the array. */
2764 if( !found ) {
2765 *units = astGrow( *units, *nunits + 1, sizeof( UnitNode * ) );
2766 if( astOK ) (*units)[ (*nunits)++ ] = node;
2767 }
2768
2769 /* If the supplied node is not an OP_LDVAR node, call this function
2770 recursively to search the argument sub-trees. */
2771 } else {
2772 for( i = 0; i < node->narg; i++ ) {
2773 LocateUnits( node->arg[ i ], units, nunits, status );
2774 }
2775 }
2776 }
2777
MakeExp(UnitNode * tree,int mathmap,int top,int * status)2778 static const char *MakeExp( UnitNode *tree, int mathmap, int top, int *status ) {
2779 /*
2780 * Name:
2781 * MakeExp
2782
2783 * Purpose:
2784 * Make an algebraic expression from a supplied tree of UnitNodes.
2785
2786 * Type:
2787 * Private function.
2788
2789 * Synopsis:
2790 * #include "unit.h"
2791 * const char *MakeExp( UnitNode *tree, int mathmap, int top, int *status )
2792
2793 * Class Membership:
2794 * Unit member function.
2795
2796 * Description:
2797 * This function produces a string holding an algebraic expression
2798 * corresponding to a supplied tree of UnitNodes.
2799
2800 * Parameters:
2801 * tree
2802 * A pointer to the UnitNode at the head of the tree to be converted
2803 * into an algebraic expression.
2804 * mathmap
2805 * If zero, format as an axis label expression. If 1, format as a
2806 * MathMap expression. If 2, format as a FITS unit string.
2807 * top
2808 * Should be non-zero for a top-level entry to this function, and
2809 * zero for a recursive entry.
2810 * status
2811 * Pointer to the inherited status variable.
2812
2813 * Returned Value:
2814 * A pointer to the cleaned expression, which should be freed using
2815 * astFree when no longer needed.
2816
2817 * Notes:
2818 * - This function returns NULL if it is invoked with the global error
2819 * status set, or if it should fail for any reason.
2820 */
2821
2822 /* Local Variables: */
2823 UnitNode *newtree;
2824 UnitNode *sunit;
2825 char *a;
2826 char *result;
2827 char buff[200];
2828 const char *arg0;
2829 const char *arg1;
2830 const char *mtxt;
2831 int larg0;
2832 int larg1;
2833 int lbuff;
2834 int mlen;
2835 int par;
2836 int tlen;
2837
2838 /* Check inherited status. */
2839 result = NULL;
2840 if( !astOK ) return result;
2841
2842 /* Modify the tree to make the resulting transformation functions more
2843 natural to human readers. */
2844 newtree = CopyTree( tree, status );
2845 ComplicateTree( &newtree, status );
2846
2847 /* If we are producing an axis label... */
2848 if( !mathmap ) {
2849
2850 /* Fix all multiplicative constants to 1.0 if they multiply an OP_LDVAR
2851 OP_SQRT or OP_POW node. This is on the assumption that the returned label
2852 should not include any simple unit scaling (e.g. if the output label would
2853 be "2.345*wavelength", we prefer simply to use "wavelength" since a scaled
2854 wavelength is still a wavelength - i.e. simple scaling does not change
2855 the dimensions of a quantity). */
2856 FixConstants( &newtree, 1, status );
2857
2858 /* Simplify the tree again to get rid of the 1.0 terms which may have
2859 been introduced by the previous line (but do not re-introduce any
2860 standardisations - removing them was the reason for calling ComplicateTree).
2861 If this simplication introduces any changes, try fixing multiplicative
2862 constants again, and so on, until no more changes occur. */
2863 while( SimplifyTree( &newtree, 0, status ) ) {
2864 FixConstants( &newtree, 1, status );
2865 }
2866
2867 }
2868
2869 /* Produce a string describing the action performed by the UnitNode at
2870 the head of the supplied tree, and then invoke this function recursively
2871 to format any arguments of the head node. */
2872
2873 /* Constant valued nodes... just format the constant in a local buffer and
2874 then copy the buffer. */
2875 if( newtree->con != AST__BAD ) {
2876 lbuff = sprintf( buff, "%.*g", DBL_DIG, newtree->con );
2877 result = astStore( NULL, buff, lbuff + 1 );
2878
2879 /* "Load Variable Value" nodes - return the variable name. If this is a
2880 recursive call to this function, and we are producing a label, append a
2881 single space before and after the name. */
2882 } else if( newtree->opcode == OP_LDVAR ) {
2883 tlen = strlen( newtree->name );
2884
2885 if( !mathmap && !top ){
2886 result = astMalloc( tlen + 3 );
2887 if( result ) {
2888 result[ 0 ] = ' ';
2889 memcpy( result + 1, newtree->name, tlen );
2890 memcpy( result + tlen + 1, " ", 2 );
2891 }
2892
2893 } else if( mathmap == 2 ) {
2894
2895 if( newtree->mult ) {
2896 mlen = newtree->mult->symlen;
2897 mtxt = newtree->mult->sym;
2898 } else {
2899 mlen = 0;
2900 mtxt = NULL;
2901 }
2902
2903 result = astMalloc( tlen + 1 + mlen );
2904 if( result ) {
2905 if( mtxt ) memcpy( result, mtxt, mlen );
2906 memcpy( result + mlen, newtree->name, tlen + 1 );
2907 }
2908
2909 } else {
2910 result = astStore( NULL, newtree->name, tlen + 1 );
2911 }
2912
2913 /* Single argument functions... place the argument in parentheses after
2914 the function name. */
2915 } else if( newtree->opcode == OP_LOG ) {
2916 arg0 = MakeExp( newtree->arg[ 0 ], mathmap, 0, status );
2917 larg0 = strlen( arg0 );
2918 if( mathmap == 1 ) {
2919 result = astMalloc( larg0 + 8 );
2920 if( result ) memcpy( result, "log10(", 7 );
2921 a = result + 6;
2922 } else {
2923 result = astMalloc( larg0 + 6 );
2924 if( result ) memcpy( result, "log(", 5 );
2925 a = result + 4;
2926 }
2927 if( result ){
2928 memcpy( a, arg0, larg0 + 1 );
2929 memcpy( a + larg0, ")", 2 );
2930 }
2931 arg0 = astFree( (void *) arg0 );
2932
2933 } else if( newtree->opcode == OP_LN ) {
2934 arg0 = MakeExp( newtree->arg[ 0 ], mathmap, 0, status );
2935 larg0 = strlen( arg0 );
2936 if( mathmap == 1 ) {
2937 result = astMalloc( larg0 + 6 );
2938 if( result ) memcpy( result, "log(", 5 );
2939 a = result + 4;
2940 } else {
2941 result = astMalloc( larg0 + 5 );
2942 if( result ) memcpy( result, "ln(", 4 );
2943 a = result + 3;
2944 }
2945 if( astOK ){
2946 memcpy( a, arg0, larg0 );
2947 memcpy( a + larg0, ")", 2 );
2948 }
2949 arg0 = astFree( (void *) arg0 );
2950
2951 } else if( newtree->opcode == OP_EXP ) {
2952 arg0 = MakeExp( newtree->arg[ 0 ], mathmap, 0, status );
2953 larg0 = strlen( arg0 );
2954 result = astMalloc( larg0 + 6 );
2955 if( result ){
2956 memcpy( result, "exp(", 5 );
2957 memcpy( result + 4, arg0, larg0 );
2958 memcpy( result + 4 + larg0, ")", 2 );
2959 }
2960 arg0 = astFree( (void *) arg0 );
2961
2962 } else if( newtree->opcode == OP_SQRT ) {
2963 arg0 = MakeExp( newtree->arg[ 0 ], mathmap, 0, status );
2964 larg0 = strlen( arg0 );
2965 result = astMalloc( larg0 + 7 );
2966 if( result ){
2967 memcpy( result, "sqrt(", 6 );
2968 memcpy( result + 5, arg0, larg0 );
2969 memcpy( result + 5 + larg0, ")", 2 );
2970 }
2971 arg0 = astFree( (void *) arg0 );
2972
2973 /* POW... the exponent (arg[1]) is always a constant and so does not need
2974 to be placed in parentheses. The first argument only needs to be
2975 placed in parentheses if it is a two arg node (except we also put it
2976 in parentheses if it is an OP_LDVAR node and "mathmap" is zero - this is
2977 because such OP_LDVAR nodes will correspond to axis labels which will
2978 have spaces before and after them which would look odd if not encloses
2979 in parentheses). */
2980 } else if( newtree->opcode == OP_POW ) {
2981
2982 arg0 = MakeExp( newtree->arg[ 0 ], mathmap, 0, status );
2983 larg0 = strlen( arg0 );
2984
2985 arg1 = MakeExp( newtree->arg[ 1 ], mathmap, 0, status );
2986 larg1 = strlen( arg1 );
2987
2988 if( newtree->arg[ 0 ]->narg == 2 ||
2989 (newtree->arg[ 0 ]->opcode == OP_LDVAR && !mathmap) ) {
2990 par = 1;
2991 result = astMalloc( larg0 + larg1 + 7 );
2992 if( result ) memcpy( result, "(", 2 );
2993 a = result + 1;
2994 } else {
2995 par = 0;
2996 result = astMalloc( larg0 + larg1 + 5 );
2997 a = result;
2998 }
2999
3000 if( result ) {
3001 memcpy( a, arg0, larg0 );
3002 a += larg0;
3003 if( par ) *(a++) = ')';
3004 memcpy( a, "**", 3 );
3005 a += 2;
3006 memcpy( a, arg1, larg1 );
3007 a += larg1;
3008 *a = 0;
3009 }
3010
3011 arg0 = astFree( (void *) arg0 );
3012 arg1 = astFree( (void *) arg1 );
3013
3014 /* DIV... the first argument (numerator) never needs to be in parentheses.
3015 The second argument (denominator) only needs to be placed in parentheses
3016 if it is a MULT node. */
3017 } else if( newtree->opcode == OP_DIV ) {
3018
3019 if( mathmap == 2 && ( sunit = ModifyPrefix( newtree, status ) ) ) {
3020 result = (char *) MakeExp( sunit, mathmap, 0, status );
3021 sunit = FreeTree( sunit, status );
3022
3023 } else {
3024 arg0 = MakeExp( newtree->arg[ 0 ], mathmap, 0, status );
3025 larg0 = strlen( arg0 );
3026
3027 arg1 = MakeExp( newtree->arg[ 1 ], mathmap, 0, status );
3028 larg1 = strlen( arg1 );
3029
3030 if( newtree->arg[ 1 ]->opcode == OP_MULT &&
3031 strchr( arg1, '*' ) ) {
3032 par = 1;
3033 result = astMalloc( larg0 + larg1 + 4 );
3034 } else {
3035 par = 0;
3036 result = astMalloc( larg0 + larg1 + 2 );
3037 }
3038
3039 if( result ) {
3040 memcpy( result, arg0, larg0 );
3041 a = result + larg0;
3042 *(a++) = '/';
3043 if( par ) *(a++) = '(';
3044 memcpy( a, arg1, larg1 );
3045 a += larg1;
3046 if( par ) *(a++) = ')';
3047 *a = 0;
3048 }
3049
3050 arg0 = astFree( (void *) arg0 );
3051 arg1 = astFree( (void *) arg1 );
3052 }
3053
3054 /* MULT... the second argument never needs to be in parentheses. The first
3055 argument only needs to be placed in parentheses if it is a DIV or POW
3056 node. */
3057 } else if( newtree->opcode == OP_MULT ) {
3058 if( mathmap == 2 && ( sunit = ModifyPrefix( newtree, status ) ) ) {
3059 result = (char *) MakeExp( sunit, mathmap, 0, status );
3060 sunit = FreeTree( sunit, status );
3061
3062 } else {
3063 arg0 = MakeExp( newtree->arg[ 0 ], mathmap, 0, status );
3064 larg0 = strlen( arg0 );
3065
3066 arg1 = MakeExp( newtree->arg[ 1 ], mathmap, 0, status );
3067 larg1 = strlen( arg1 );
3068
3069 /* If this is a top-level entry and we are producing an axis label, do
3070 not include any constant multiplicative terms. */
3071 if( top && !mathmap ) {
3072 if( newtree->arg[ 0 ]->con != AST__BAD ) arg0 = astFree( (void *) arg0 );
3073 if( newtree->arg[ 1 ]->con != AST__BAD ) arg1 = astFree( (void *) arg1 );
3074 }
3075
3076 /* If we have two arguments, concatentate them, placing the operands in
3077 parentheses if necessary. */
3078 if( arg0 && arg1 ) {
3079
3080 if( ( newtree->arg[ 0 ]->opcode == OP_DIV &&
3081 strchr( arg0, '/' ) ) ||
3082 ( newtree->arg[ 0 ]->opcode == OP_POW &&
3083 strstr( arg0, "**" ) ) ) {
3084 par = 1;
3085 result = astMalloc( larg0 + larg1 + 4 );
3086 if( result ) result[ 0 ] = '(';
3087 a = result + 1;
3088 } else {
3089 par = 0;
3090 result = astMalloc( larg0 + larg1 + 2 );
3091 a = result;
3092 }
3093
3094 if( result ) {
3095 memcpy( a, arg0, larg0 );
3096 a += larg0;
3097 if( par ) *(a++) = ')';
3098 *(a++) = '*';
3099 memcpy( a, arg1, larg1 );
3100 a += larg1;
3101 *a = 0;
3102 }
3103
3104 arg0 = astFree( (void *) arg0 );
3105 arg1 = astFree( (void *) arg1 );
3106
3107 /* If we do not have two arguments, just return the one we do have. */
3108 } else if( arg0 ){
3109 result = (char *) arg0;
3110
3111 } else {
3112 result = (char *) arg1;
3113 }
3114 }
3115 }
3116
3117 /* Free the complicated tree. */
3118 newtree = FreeTree( newtree, status );
3119
3120 /* Free the returned string if an error has occurred. */
3121 if( !astOK ) result = astFree( result );
3122
3123 /* Return the result. */
3124 return (const char *) result;
3125 }
3126
MakeKnownUnit(const char * sym,const char * label,const char * exp,int * status)3127 static void MakeKnownUnit( const char *sym, const char *label, const char *exp, int *status ){
3128 /*
3129 * Name:
3130 * MakeKnownUnit
3131
3132 * Purpose:
3133 * Create a KnownUnit structure describing a known unit.
3134
3135 * Type:
3136 * Private function.
3137
3138 * Synopsis:
3139 * #include "unit.h"
3140 * void MakeKnownUnit( const char *sym, const char *label, const char *exp, int *status )
3141
3142 * Class Membership:
3143 * Unit member function.
3144
3145 * Description:
3146 * This function creates a KnownUnit structure decribing a known unit,
3147 * and adds it to the head of the linked list of known units stored in
3148 * a module variable.
3149
3150 * Parameters:
3151 * sym
3152 * A pointer to a string which can be used as a symbol to represent
3153 * the new named unit. Once defined, this symbol can be included within
3154 * the definition of other derived units. The string should contain
3155 * only alphabetical characters (no digits, spaces, punctuation,
3156 * etc). Symbols are case sensitive (e.g. "s" is second, but "S" is
3157 * Siemens). The string should not include any multiplier prefix.
3158 * label
3159 * Pointer to a null terminated string containing the label for
3160 * the required units. No restriction on content.
3161 * exp
3162 * This should be a pointer to a null terminated string containing
3163 * a definition of the required unit. See the description of the
3164 * "in" and "out" parameters for the astUnitMapper function.
3165 *
3166 * A NULL pointer or a blank string may supplied for "exp", which
3167 * is interpreted as a request for a new basic unit to be created with
3168 * the symbol and label given by the other parameters.
3169 * status
3170 * Pointer to the inherited status variable.
3171
3172 * Notes:
3173 * - The supplied symbol and label strings are not copied. The
3174 * supplied pointers are simply stored in the returned structure.
3175 * Therefore the strings to which the pointers point should not be
3176 * modified after this function returned (in fact this function is
3177 * always called with literal strings for these arguments).
3178 */
3179
3180 /* Local Variables: */
3181 KnownUnit *result;
3182
3183 /* Check the global error status. */
3184 if( !astOK ) return;
3185
3186 /* Indicate that subsequent memory allocations may never be freed (other
3187 than by any AST exit handler). */
3188 astBeginPM;
3189
3190 /* Allocate memory for the structure, and check the returned pointer can
3191 be used safely. */
3192 result = astMalloc( sizeof( KnownUnit ) );
3193 if( astOK ) {
3194
3195 /* In case of errors, first nullify the pointer to the next KnownUnit. */
3196 result->next = NULL;
3197
3198 /* Store the supplied label and symbol pointers. */
3199 result->sym = sym;
3200 result->label = label;
3201
3202 /* Store the length of the symbol (without the trailing null character). */
3203 result->symlen = strlen( sym );
3204
3205 /* Store the length of the label (without the trailing null character). */
3206 result->lablen = strlen( label );
3207
3208 /* Create a tree of UnitNodes describing the unit if an expression was
3209 supplied. */
3210 result->head = exp ? CreateTree( exp, 1, 0, status ) : NULL;
3211
3212 /* Unit aliases are replaced in use by the KnownUnit pointed to by the
3213 "use" component of the structure. Indicate this KnownUnitis not an
3214 alias by setting its "use" component NULL. */
3215 result->use = NULL;
3216 }
3217
3218 /* Mark the end of the section in which memory allocations may never be
3219 freed (other than by any AST exit handler). */
3220 astEndPM;
3221
3222 /* If an error has occurred, free any returned structure. */
3223 if( !astOK ) {
3224 result->head = FreeTree( result->head, status );
3225 result = astFree( result ) ;
3226
3227 /* Otherwise, add the new KnownUnit to the head of the linked list of
3228 known units. */
3229 } else {
3230 result->next = known_units;
3231 known_units = result;
3232 }
3233
3234 }
3235
MakeMapping(UnitNode * tree,int * status)3236 static AstMapping *MakeMapping( UnitNode *tree, int *status ) {
3237 /*
3238 * Name:
3239 * MakeMapping
3240
3241 * Purpose:
3242 * Create a new Mapping from a given tree of UnitNodes.
3243
3244 * Type:
3245 * Private function.
3246
3247 * Synopsis:
3248 * #include "unit.h"
3249 * AstMapping *MakeMapping( UnitNode *tree )
3250
3251 * Class Membership:
3252 * Unit member function.
3253
3254 * Description:
3255 * This function creates a Mapping with a forward transformation equal
3256 * to the transformation described by the tree of UnitNodes. The head
3257 * node of the tree corresponds to the output of the Mapping.
3258
3259 * Parameters:
3260 * tree
3261 * The UnitNode at the head of the tree to be used. It should have
3262 * exactly one OP_LDVAR node, and should have been simplified using
3263 * the SimplifyTree function.
3264
3265 * Returned Value:
3266 * A pointer to the Mapping. Its Nin and Nout attributes will both be 1.
3267
3268 * Notes:
3269 * - A value of NULL will be returned if this function is invoked with
3270 * the global error status set, or if it should fail for any reason.
3271 */
3272
3273 /* Local Variables: */
3274 AstMapping *result;
3275 UnitNode *inv;
3276 UnitNode *src;
3277 const char *fwdexp;
3278 char *fwdfun;
3279 const char *invexp;
3280 char *invfun;
3281 int lfwd;
3282 int linv;
3283
3284 /* Initialise. */
3285 result = NULL;
3286
3287 /* Check the inherited status. */
3288 if( !astOK ) return result;
3289
3290 /* First see if a UnitMap can be used to represent the Mapping from
3291 input units to output units. This will be the case if the supplied tree
3292 consists of a aingle OP_LDVAR node (corresponding to the input units). */
3293 if( tree->opcode == OP_LDVAR ) {
3294 result = (AstMapping *) astUnitMap( 1, "", status );
3295
3296 /* Now see if a UnitMap or ZoomMap can be used to represent the Mapping from
3297 input units to output units. This will be the case if the supplied tree
3298 consists of a OP_MULT node with one constant argument and on OP_LDVAR
3299 argument (corresponding to the input units). The standardisation done by
3300 SimplifyTree will have ensured that the constant will be argument 0
3301 (and will also have converted "x/k" trees into "(1/k)*x" trees). */
3302 } else if( tree->opcode == OP_MULT &&
3303 tree->arg[ 0 ]->con != AST__BAD &&
3304 tree->arg[ 1 ]->opcode == OP_LDVAR ) {
3305
3306 if( tree->arg[ 0 ]->con == 1.0 ) {
3307 result = (AstMapping *) astUnitMap( 1, "", status );
3308 } else {
3309 result = (AstMapping *) astZoomMap( 1, tree->arg[ 0 ]->con, "", status );
3310 }
3311
3312 /* For other trees we need to create a MathMap. */
3313 } else {
3314
3315 /* Format the supplied tree as an algebraic expression, and get its length. */
3316 fwdexp = MakeExp( tree, 1, 1, status );
3317 lfwd = strlen( fwdexp );
3318
3319 /* The MathMap constructor requires the forward and inverse
3320 transformation functions to be specified as equations (i.e. including an
3321 equals sign). We use the output variable name "output_units" (the
3322 astUnitMapper function creates the supplied tree usign the variable
3323 name "input_units" ). */
3324 lfwd += 13;
3325
3326 /* Invert the supplied tree and create an algebraic expression from it. */
3327 src = NewNode( NULL, OP_LDVAR, status );
3328 if( astOK ) src->name = astStore( NULL, "output_units", 13 );
3329 inv = InvertTree( tree, src, status );
3330 if( !inv ) {
3331 src = FreeTree( src, status );
3332 astError( AST__BADUN, "MakeMapping(Unit): Failed to invert "
3333 "supplied tree '%s' (internal AST programming error).", status,
3334 fwdexp );
3335
3336 /* If inverted succesfully (which it should be since astUnitMapper should
3337 have checked this)... */
3338 } else {
3339
3340 /* Format the inverted tree as an algebraic expression, and get its
3341 length, adding on extra characters for the variable name ("input_units")
3342 and equals sign. */
3343 invexp = MakeExp( inv, 1, 1, status );
3344 linv = strlen( invexp );
3345 linv += 12;
3346
3347 /* Allocate memory for the transformation functions, plus an extra
3348 character for the trailing null. */
3349 fwdfun = astMalloc( lfwd + 1 );
3350 invfun = astMalloc( linv + 1 );
3351 if( invfun ) {
3352 memcpy( fwdfun, "output_units=", 14 );
3353 memcpy( invfun, "input_units=", 13 );
3354
3355 /* Append the expressions following the equals signs. */
3356 strcpy( fwdfun + 13, fwdexp );
3357 strcpy( invfun + 12, invexp );
3358
3359 /* Create the MathMap. */
3360 result = (AstMapping *) astMathMap( 1, 1, 1,
3361 (const char **) &fwdfun, 1,
3362 (const char **) &invfun,
3363 "SimpFI=1,SimpIF=1", status );
3364 }
3365
3366 /* Free resources. */
3367 inv = FreeTree( inv, status );
3368 fwdfun = astFree( fwdfun );
3369 invfun = astFree( invfun );
3370 invexp = astFree( (void *) invexp );
3371 }
3372 fwdexp = astFree( (void *) fwdexp );
3373 }
3374
3375 /* Free any result if an error occurred. */
3376 if( !astOK ) result = astAnnul( result );
3377
3378 /* Return the answer. */
3379 return result;
3380 }
3381
MakeLabelTree(const char * lab,int nc,int * status)3382 static UnitNode *MakeLabelTree( const char *lab, int nc, int *status ){
3383 /*
3384 * Name:
3385 * MakeLabelTree
3386
3387 * Purpose:
3388 * Convert an axis label into a tree of UnitNodes.
3389
3390 * Type:
3391 * Private function.
3392
3393 * Synopsis:
3394 * #include "unit.h"
3395 * UnitNode *MakeLabelTree( const char *lab, int nc, int *status )
3396
3397 * Class Membership:
3398 * Unit member function.
3399
3400 * Description:
3401 * This function converts an axis label into a tree of UnitNodes.
3402 * It is assumed the supplied label represents some "basic" label
3403 * modified by the application of one or more single function arguments
3404 * and/or exponentiation operators. The (single) OP_LDVAR node in the
3405 * returned tree refers to the basic label (it is stored as the "name"
3406 * component of UnitNode structure).
3407
3408 * Parameters:
3409 * lab
3410 * The label expression.
3411 * nc
3412 * The number of characters from "lab" to use.
3413 * status
3414 * Pointer to the inherited status variable.
3415
3416 * Returned Value:
3417 * A pointer to a UnitNode which forms the head of a tree of UnitNodes
3418 * representing the supplied label expression.
3419
3420 * Notes:
3421 * - A NULL value is returned if this function is invoked with the
3422 * global error status set or if it should fail for any reason.
3423 */
3424
3425 /* Local Variables: */
3426 Oper op;
3427 UnitNode *result;
3428 char buff[ 10 ];
3429 const char *c;
3430 const char *exp;
3431 int depth;
3432 int i;
3433 int oplen;
3434 int n;
3435 double con;
3436
3437 /* Initialise */
3438 result = NULL;
3439 oplen = 0;
3440
3441 /* Check the global error status, and that we have a string. */
3442 if ( !astOK || !lab || !nc ) return result;
3443
3444 /* Get a pointer to the first non-blank character, and store the number of
3445 characters to examine (this excludes any trailing white space). */
3446 exp = lab;
3447 while( isspace( *exp ) ) exp++;
3448 c = lab + nc - 1;
3449 while( c >= exp && isspace( *c ) ) c--;
3450 nc = c - exp + 1;
3451
3452 /* Scan through the supplied string looking for the first pow operator at
3453 zero depth of nesting within parentheses. */
3454 depth = 0;
3455 c = exp;
3456 i = 0;
3457 op = OP_NULL;
3458 while( i < nc && *c ){
3459
3460 /* If this character is an opening parenthesis, increment the depth of
3461 nesting. */
3462 if( *c == '(' ) {
3463 depth++;
3464
3465 /* If this character is an closing parenthesis, decrement the depth of
3466 nesting. Report an error if it ever goes negative. */
3467 } else if( *c == ')' ) {
3468 depth--;
3469 if( depth < 0 && astOK ) {
3470 astError( AST__BADUN, "Missing opening parenthesis." , status);
3471 break;
3472 }
3473
3474 /* Ignore all other characters unless they are at zero depth of nesting.
3475 Also ignore spaces. */
3476 } else if( depth == 0 && !isspace( *c ) ) {
3477
3478 /* Compare the next part of the string with each of the "pow" operators. */
3479 if( !strncmp( c, "**", 2 ) ) {
3480 op = OP_POW;
3481 oplen = 2;
3482 } else if( *c == '^' ) {
3483 op = OP_POW;
3484 oplen = 1;
3485 }
3486
3487 /* If an operator was found, break out of the loop. */
3488 if( op != OP_NULL ) break;
3489 }
3490
3491 /* Pass on to check the next character. */
3492 i++;
3493 c++;
3494 }
3495
3496 /* If a "pow" operator was found, the strings on either side of it should be
3497 valid unit expressions, in which case we use this routine recursively to
3498 create corresponding trees of UnitNodes. */
3499 if( op != OP_NULL ) {
3500
3501 /* Create a UnitNode for the operator. */
3502 result = NewNode( NULL, op, status );
3503 if( astOK ) {
3504
3505 /* Create a tree of unit nodes from the string which precedes the binary
3506 operator. Report an error if it cannot be done. */
3507 result->arg[ 0 ] = MakeLabelTree( exp, i, status );
3508 if( !result->arg[ 0 ] && astOK ) {
3509 for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
3510 buff[ oplen ] = 0;
3511 astError( AST__BADUN, "Missing operand before '%s'.", status, buff );
3512 }
3513
3514 /* Create a tree of unit nodes from the string which follows the binary
3515 operator. Report an error if it cannot be done. */
3516 result->arg[ 1 ] = MakeLabelTree( c + oplen, nc - i - oplen, status );
3517 if( !result->arg[ 1 ] && astOK ) {
3518 for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
3519 buff[ oplen ] = 0;
3520 astError( AST__BADUN, "Missing operand after '%s'.", status, buff );
3521 }
3522 }
3523
3524 /* If no binary operator was found at depth zero, see if the supplied string
3525 starts with a function name (the only legal place for a function name
3526 given that the string has no binary operators at depth zero). */
3527 } else {
3528 if( !strncmp( exp, "sqrt(", 5 ) || !strncmp( exp, "SQRT(", 5 ) ) {
3529 op = OP_SQRT;
3530 oplen = 4;
3531 } else if( !strncmp( exp, "exp(", 4 ) || !strncmp( exp, "EXP(", 4 ) ) {
3532 op = OP_EXP;
3533 oplen = 3;
3534 } else if( !strncmp( exp, "ln(", 3 ) || !strncmp( exp, "LN(", 3 ) ) {
3535 op = OP_LN;
3536 oplen = 2;
3537 } else if( !strncmp( exp, "log(", 4 ) || !strncmp( exp, "LOG(", 4 ) ) {
3538 op = OP_LOG;
3539 oplen = 3;
3540 }
3541
3542 /* If a function was found, the string following the function name
3543 (including the opening parenthesis) should form a legal units
3544 expresssion (all the supported functions take a single argument and
3545 so we do not need to worry about comma-separated lists of function
3546 arguments). Use this routine recursively to create a tree of UnitNodes
3547 from the string which forms the function argument. */
3548 if( op != OP_NULL ) {
3549
3550 /* Create a UnitNode for the function. */
3551 result = NewNode( NULL, op, status );
3552 if( astOK ) {
3553
3554 /* Create a tree of unit nodes from the string which follows the function
3555 name. Report an error if it cannot be done. */
3556 result->arg[ 0 ] = MakeLabelTree( exp + oplen, nc - oplen, status );
3557 if( !result->arg[ 0 ] && astOK ) {
3558 for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
3559 buff[ oplen ] = 0;
3560 astError( AST__BADUN, "Missing argument for '%s'.", status, buff );
3561 }
3562 }
3563
3564 /* Arrive here if the supplied string does not contain a POW operator
3565 or function at depth zero. Check to see if the whole string is contained
3566 within parentheses, In which we interpret the contents of the
3567 parentheses as a units expression. It is safe simply to check the
3568 first and last characters (a string like "(fred)(Harry)" is not a
3569 legal possibility since there should be an operator in the middle).*/
3570 } else if( nc > 0 && ( exp[ 0 ] == '(' && exp[ nc - 1 ] == ')' ) ) {
3571 result = MakeLabelTree( exp + 1, nc - 2, status );
3572
3573 /* Does the string begin with a numerical constant? */
3574 } else if( ConStart( exp, &con, &n, status ) == 1 ) {
3575
3576 /* If the entire string was a numerical constant, represent it by a LDCON
3577 node. */
3578 if( n == nc ) {
3579 result = NewNode( NULL, OP_LDCON, status );
3580 if( astOK ) result->con = con;
3581
3582 /* If there was anything following the numerical constant, report an
3583 error. */
3584 } else if( astOK ){
3585 astError( AST__BADUN, "Missing operator after "
3586 "numerical string '%.*s'.", status, n, exp );
3587 }
3588
3589 /* The only legal possibility left is that the string represents the basic
3590 label. Create an OP_LDVAR node for it and store the basic label as
3591 the node name, omitting any enclosing white space. */
3592 } else {
3593 result = NewNode( NULL, OP_LDVAR, status );
3594 if( astOK ) {
3595 result->name = astStore( NULL, exp, nc + 1 );
3596 if( astOK ) ( (char *) result->name)[ nc ] = 0;
3597 }
3598 }
3599 }
3600
3601 /* Free any returned tree if an error has occurred. */
3602 if( !astOK ) result = FreeTree( result, status );
3603
3604 /* Return the result. */
3605 return result;
3606 }
3607
MakeTree(const char * exp,int nc,int lock,int * status)3608 static UnitNode *MakeTree( const char *exp, int nc, int lock, int *status ){
3609 /*
3610 * Name:
3611 * MakeTree
3612
3613 * Purpose:
3614 * Convert an algebraic units expression into a tree of UnitNodes.
3615
3616 * Type:
3617 * Private function.
3618
3619 * Synopsis:
3620 * #include "unit.h"
3621 * UnitNode *MakeTree( const char *exp, int nc, int lock, int *status )
3622
3623 * Class Membership:
3624 * Unit member function.
3625
3626 * Description:
3627 * This function converts an algebraic units expression into a tree of
3628 * UnitNodes. It is a service routine for CreateTree. The roots of the
3629 * returned tree (i.e. the LDVAR nodes) refer to the unit symbols
3630 * contained within the supplied expression (i.e. definitions of these
3631 * units are not grafted onto the tree in place of the original nodes,
3632 * as is done by CreateTree).
3633
3634 * Parameters:
3635 * exp
3636 * The units expression. This should not include any leading or
3637 * trailing spaces.
3638 * nc
3639 * The number of characters from "exp" to use.
3640 * lock
3641 * Use a mutex to guard access to the KnownUnits list?
3642 * status
3643 * Pointer to the inherited status variable.
3644
3645 * Returned Value:
3646 * A pointer to a UnitNode which forms the head of a tree of UnitNodes
3647 * representing the supplied unit expression.
3648
3649 * Notes:
3650 * - A NULL value is returned if this function is invoked with the
3651 * global error status set or if it should fail for any reason.
3652 */
3653
3654 /* Local Variables: */
3655 KnownUnit *munit;
3656 KnownUnit *unit;
3657 Multiplier *mmult;
3658 Multiplier *mult;
3659 Oper op;
3660 UnitNode *result;
3661 char buff[ 10 ];
3662 char d;
3663 const char *c;
3664 double con;
3665 int depth;
3666 int i;
3667 int l;
3668 int maxlen;
3669 int n;
3670 int oplen;
3671 int plural;
3672
3673 /* Initialise */
3674 result = NULL;
3675
3676 /* Check the global error status, and that we have a string. */
3677 if ( !astOK || !exp || nc <= 0 ) return result;
3678
3679 /* Scan through the supplied string from the end to the start looking for
3680 the last multiplication or division operator at zero depth of nesting
3681 within parentheses. We go backwards through the string in order to
3682 give the correct priority to multiple division operators (i.e. "a/b/c"
3683 needs to be interpreted as "(a/b)/c", not "a/(b/c)"). */
3684 op = OP_NULL;
3685 oplen = 1;
3686 depth = 0;
3687 c = exp + nc - 1;
3688 i = nc - 1;
3689 while( i >= 0 ){
3690
3691 /* If this character is an opening parenthesis, decrement the depth of
3692 nesting. Report an error if it ever goes negative. */
3693 if( *c == '(' ) {
3694 depth--;
3695 if( depth < 0 && astOK ) {
3696 astError( AST__BADUN, "Missing closing parenthesis." , status);
3697 break;
3698 }
3699
3700 /* An opening parenthesis at level zero must always be either the first
3701 character in the string, or be preceded by the name of a function, or
3702 be preceded by an operator. If none of these are true, assume there is
3703 an implicit multiplication operator before the parenthesis. */
3704 if( depth == 0 && i > 0 ) {
3705 d = *( c - 1 );
3706 if( d != '*' && d != '/' && d != '^' && d != '.' && d != ' ' &&
3707 !EndsWith( c, i + 1, "sqrt(", status ) && !EndsWith( c, i + 1, "exp(", status ) &&
3708 !EndsWith( c, i + 1, "ln(", status ) && !EndsWith( c, i + 1, "log(", status ) ) {
3709 op = OP_MULT;
3710 oplen = 0;
3711 break;
3712 }
3713 }
3714
3715 /* If this character is an closing parenthesis, increment the depth of
3716 nesting. */
3717 } else if( *c == ')' ) {
3718 depth++;
3719
3720 /* A closing parenthesis at level zero must always be either the last
3721 character in the string, or be followed by an operator. If neither of
3722 these are true, assume there is an implicit multiplication operator. */
3723 if( depth == 1 && i < nc - 1 ) {
3724 d = *(c+1);
3725 if( d != '*' && d != '/' && d != '^' && d != '.' && d != ' ') {
3726 op = OP_MULT;
3727 oplen = 0;
3728
3729 /* Correct "i" so that it gives the length of the left hand operand of
3730 the implicit MULT operator, correct "c" so that it points to the first
3731 character in the right hand operand, and leave the loop. */
3732 i++;
3733 c++;
3734 break;
3735 }
3736 }
3737
3738 /* Ignore all other characters unless they are at zero depth of nesting. */
3739 } else if( depth == 0 ) {
3740
3741 /* Compare the next part of the string with each of the multiplication
3742 and division operators. */
3743 if( *c == '/' ) {
3744 op = OP_DIV;
3745
3746 } else if( *c == ' ' ) {
3747 op = OP_MULT;
3748
3749 /* An asterisk is only treated as a multiplication symbol if it does not occur
3750 before or after another asterisk. */
3751 } else if( *c == '*' ) {
3752 if( c == exp ) {
3753 if( *(c+1) != '*' ) op = OP_MULT;
3754 } else if( i == nc - 1 ) {
3755 if( *(c-1) != '*' ) op = OP_MULT;
3756 } else {
3757 if( *(c+1) != '*' && *(c-1) != '*' ) op = OP_MULT;
3758 }
3759
3760 /* A dot is only treated as a multiplication symbol if it does not occur
3761 between two digits. */
3762 } else if( *c == '.' ) {
3763 if( ( c == exp || !isdigit( *(c-1) ) ) &&
3764 ( i == nc - 1 || !isdigit( *(c+1) ) ) ) {
3765 op = OP_MULT;
3766 }
3767 }
3768 }
3769
3770 /* If an operator was found, break out of the loop. */
3771 if( op != OP_NULL ) break;
3772
3773 /* Pass on to check the next character. */
3774 i--;
3775 c--;
3776 }
3777
3778 /* If a multiplication or division operator was found, the strings on either
3779 side of it should be valid unit expressions, in which case we use this
3780 routine recursively to create corresponding trees of UnitNodes. */
3781 if( op != OP_NULL ) {
3782
3783 /* Create a UnitNode for the binary operator. */
3784 result = NewNode( NULL, op, status );
3785 if( astOK ) {
3786
3787 /* Create a tree of unit nodes from the string which precedes the binary
3788 operator. Report an error if it cannot be done. */
3789 result->arg[ 0 ] = MakeTree( exp, i, lock, status );
3790 if( !result->arg[ 0 ] && astOK ) {
3791 for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
3792 buff[ oplen ] = 0;
3793 astError( AST__BADUN, "Missing operand before '%s'.", status, buff );
3794 }
3795
3796 /* Create a tree of unit nodes from the string which follows the binary
3797 operator. Report an error if it cannot be done. */
3798 result->arg[ 1 ] = MakeTree( c + oplen, nc - i - oplen, lock, status );
3799 if( !result->arg[ 1 ] && astOK ) {
3800 for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
3801 buff[ oplen ] = 0;
3802 astError( AST__BADUN, "Missing operand after '%s'.", status, buff );
3803 }
3804 }
3805
3806 /* If no multiplication or division operator was found at depth zero, check
3807 that the final depth of nesting was zero. Report an error if not. */
3808 } else if( depth > 0 && astOK ) {
3809 astError( AST__BADUN, "Missing opening parenthesis." , status);
3810
3811 /* Otherwise check for a "Pow" operator at depth zero. */
3812 } else {
3813
3814 /* Scan through the supplied string looking for the first pow operator at
3815 zero depth of nesting within parentheses. */
3816 depth = 0;
3817 c = exp;
3818 i = 0;
3819 while( i < nc && *c ){
3820
3821 /* If this character is an opening parenthesis, increment the depth of
3822 nesting. */
3823 if( *c == '(' ) {
3824 depth++;
3825
3826 /* If this character is an closing parenthesis, decrement the depth of
3827 nesting. Report an error if it ever goes negative. */
3828 } else if( *c == ')' ) {
3829 depth--;
3830 if( depth < 0 && astOK ) {
3831 astError( AST__BADUN, "Missing opening parenthesis." , status);
3832 break;
3833 }
3834
3835 /* Ignore all other characters unless they are at zero depth of nesting. */
3836 } else if( depth == 0 ) {
3837
3838 /* Compare the next part of the string with each of the "pow" operators. */
3839 if( !strncmp( c, "**", 2 ) ) {
3840 op = OP_POW;
3841 oplen = 2;
3842 } else if( *c == '^' ) {
3843 op = OP_POW;
3844 oplen = 1;
3845 }
3846
3847 /* If an operator was found, break out of the loop. */
3848 if( op != OP_NULL ) break;
3849 }
3850
3851 /* Pass on to check the next character. */
3852 i++;
3853 c++;
3854 }
3855
3856 /* If a "pow" operator was found, the strings on either side of it should be
3857 valid unit expressions, in which case we use this routine recursively to
3858 create corresponding trees of UnitNodes. */
3859 if( op != OP_NULL ) {
3860
3861 /* Create a UnitNode for the operator. */
3862 result = NewNode( NULL, op, status );
3863 if( astOK ) {
3864
3865 /* Create a tree of unit nodes from the string which precedes the binary
3866 operator. Report an error if it cannot be done. */
3867 result->arg[ 0 ] = MakeTree( exp, i, lock, status );
3868 if( !result->arg[ 0 ] && astOK ) {
3869 for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
3870 buff[ oplen ] = 0;
3871 astError( AST__BADUN, "Missing operand before '%s'.", status, buff );
3872 }
3873
3874 /* Create a tree of unit nodes from the string which follows the binary
3875 operator. Report an error if it cannot be done. */
3876 result->arg[ 1 ] = MakeTree( c + oplen, nc - i - oplen, lock, status );
3877 if( !result->arg[ 1 ] && astOK ) {
3878 for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
3879 buff[ oplen ] = 0;
3880 astError( AST__BADUN, "Missing operand after '%s'.", status, buff );
3881 }
3882 }
3883
3884 /* If no binary operator was found at depth zero, see if the supplied string
3885 starts with a function name (the only legal place for a function name
3886 given that the string has no binary operators at depth zero). */
3887 } else {
3888 if( !strncmp( exp, "sqrt(", 5 ) || !strncmp( exp, "SQRT(", 5 ) ) {
3889 op = OP_SQRT;
3890 oplen = 4;
3891 } else if( !strncmp( exp, "exp(", 4 ) || !strncmp( exp, "EXP(", 4 ) ) {
3892 op = OP_EXP;
3893 oplen = 3;
3894 } else if( !strncmp( exp, "ln(", 3 ) || !strncmp( exp, "LN(", 3 ) ) {
3895 op = OP_LN;
3896 oplen = 2;
3897 } else if( !strncmp( exp, "log(", 4 ) || !strncmp( exp, "LOG(", 4 ) ) {
3898 op = OP_LOG;
3899 oplen = 3;
3900 }
3901
3902 /* If a function was found, the string following the function name
3903 (including the opening parenthesis) should form a legal units
3904 expresssion (all the supported functions take a single argument and
3905 so we do not need to worry about comma-separated lists of function
3906 arguments). Use this routine recursively to create a tree of UnitNodes
3907 from the string which forms the function argument. */
3908 if( op != OP_NULL ) {
3909
3910 /* Create a UnitNode for the function. */
3911 result = NewNode( NULL, op, status );
3912 if( astOK ) {
3913
3914 /* Create a tree of unit nodes from the string which follows the function
3915 name. Report an error if it cannot be done. */
3916 result->arg[ 0 ] = MakeTree( exp + oplen, nc - oplen, lock, status );
3917 if( !result->arg[ 0 ] && astOK ) {
3918 for( i = 0; i < oplen; i++ ) buff[ i ] = c[ i ];
3919 buff[ oplen ] = 0;
3920 astError( AST__BADUN, "Missing argument for '%s'.", status, buff );
3921 }
3922 }
3923
3924 /* Arrive here if the supplied string does not contain a binary operator
3925 or function at depth zero. Check to see if the whole string is contained
3926 within parentheses, In which we interpret the contents of the
3927 parentheses as a units expression. It is safe simply to check the
3928 first and last characters (a string like "(fred)(Harry)" is not a
3929 legal possibility since there should be an operator in the middle).*/
3930 } else if( exp[ 0 ] == '(' && exp[ nc - 1 ] == ')' ) {
3931 result = MakeTree( exp + 1, nc - 2, lock, status );
3932
3933 /* Does the string begin with a numerical constant? */
3934 } else if( ConStart( exp, &con, &n, status ) == 1 ) {
3935
3936 /* If the entire string was a numerical constant, represent it by a LDCON
3937 node. */
3938 if( n == nc ) {
3939 result = NewNode( NULL, OP_LDCON, status );
3940 if( astOK ) result->con = con;
3941
3942 /* If there was anything following the numerical constant, report an
3943 error. */
3944 } else if( astOK ){
3945 astError( AST__BADUN, "Missing operator after "
3946 "numerical string '%.*s'.", status, n, exp );
3947 }
3948
3949 /* Does the string represent one of the named constants? If so represent it
3950 by a an appropriate operator. */
3951 } else if( nc == 2 && ( !strncmp( exp, "pi", 2 ) ||
3952 !strncmp( exp, "PI", 2 ) ) ) {
3953 result = NewNode( NULL, OP_LDPI, status );
3954
3955 } else if( nc == 1 && ( !strncmp( exp, "e", 1 ) ||
3956 !strncmp( exp, "E", 1 ) ) ) {
3957 result = NewNode( NULL, OP_LDE, status );
3958
3959 /* The only legal possibility left is that the string represents the name
3960 of a basic unit, possibly prefixed by a multiplier character. */
3961 } else {
3962
3963 /* See if the string ends with the symbol for any of the known basic
3964 units. If it matches more than one basic unit, choose the longest.
3965 First ensure descriptions of the known units are available. */
3966 mmult = NULL;
3967 plural = 0;
3968 while( 1 ) {
3969 unit = GetKnownUnits( lock, status );
3970
3971 maxlen = -1;
3972 munit = NULL;
3973 while( unit ) {
3974 if( SplitUnit( exp, nc, unit->sym, 1, &mult, &l, status ) ) {
3975 if( l > maxlen ) {
3976 maxlen = l;
3977 munit = unit;
3978 mmult = mult;
3979 }
3980 }
3981 unit = unit->next;
3982 }
3983
3984 /* If the above did not produce a match, try matching the unit symbol
3985 case insensitive. */
3986 if( !munit ) {
3987 unit = GetKnownUnits( lock, status );
3988 while( unit ) {
3989 if( SplitUnit( exp, nc, unit->sym, 0, &mult, &l, status ) ) {
3990 if( l > maxlen ) {
3991 maxlen = l;
3992 munit = unit;
3993 mmult = mult;
3994 }
3995 }
3996 unit = unit->next;
3997 }
3998 }
3999
4000 /* If the above did not produce a match, try matching the unit label
4001 case insensitive. */
4002 if( !munit ) {
4003 unit = GetKnownUnits( lock, status );
4004 while( unit ) {
4005 if( SplitUnit( exp, nc, unit->label, 0, &mult, &l, status ) ) {
4006 if( l > maxlen ) {
4007 maxlen = l;
4008 munit = unit;
4009 mmult = mult;
4010 }
4011 }
4012 unit = unit->next;
4013 }
4014 }
4015
4016 /* If we still do not have a match, and if the string ends with "s", try
4017 removing the "s" (which could be a plural as in "Angstroms") and
4018 trying again. */
4019 if( !munit && nc > 1 && !plural &&
4020 ( exp[ nc - 1 ] == 's' || exp[ nc - 1 ] == 'S' ) ) {
4021 plural = 1;
4022 nc--;
4023 } else {
4024 break;
4025 }
4026 }
4027 if( plural ) nc++;
4028
4029 /* If a known unit and multiplier combination was found, create an
4030 OP_LDVAR node from it. */
4031 unit = munit;
4032 mult = mmult;
4033 if( unit ) {
4034
4035 /* If the unit is an alias for another unit, it will have a non-NULL
4036 value for its "use" component.In this case, use the unit for which the
4037 identified unit is an alias. */
4038 result = NewNode( NULL, OP_LDVAR, status );
4039 if( astOK ) {
4040 result->unit = unit->use ? unit->use : unit;
4041 result->mult = mult;
4042 result->name = astStore( NULL, result->unit->sym, result->unit->symlen + 1 );
4043 }
4044
4045 /* If no known unit and multiplier combination was found, we assume the
4046 string represents a new user-defined basic unit, possibly preceded by a
4047 standard multiplier prefix. */
4048 } else {
4049
4050 /* Check the string to see if starts with a known multiplier prefix (but
4051 do not allow the multiplier to account for the entire string). */
4052 mult = GetMultipliers( status );
4053 c = exp;
4054 while( mult ) {
4055 n = nc - mult->symlen;
4056 if( n > 0 && !strncmp( exp, mult->sym, mult->symlen ) ) {
4057 c += mult->symlen;
4058 break;
4059 }
4060 mult = mult->next;
4061 }
4062 if( !mult ) n = nc;
4063
4064 /* Check there are no illegal characters in the following string. */
4065 for( i = 0; i < n && astOK; i++ ) {
4066 if( !isalpha( c[ i ] ) ) {
4067 astError( AST__BADUN, "Illegal character '%c' found.", status, c[ i ] );
4068 break;
4069 }
4070 }
4071
4072 /* If succesfull, create an OP_LDVAR node for th user-defined basic unit. */
4073 if( astOK ) {
4074 result = NewNode( NULL, OP_LDVAR, status );
4075 if( astOK ) {
4076 result->mult = mult;
4077 result->name = astStore( NULL, c, n + 1 );
4078 if( astOK ) ( (char *) result->name)[ n ] = 0;
4079 }
4080 }
4081 }
4082 }
4083 }
4084 }
4085
4086 /* Free any returned tree if an error has occurred. */
4087 if( !astOK ) result = FreeTree( result, status );
4088
4089 /* Return the result. */
4090 return result;
4091 }
4092
MakeUnitAlias(const char * sym,const char * alias,int * status)4093 static void MakeUnitAlias( const char *sym, const char *alias, int *status ){
4094 /*
4095 * Name:
4096 * MakeUnitAlias
4097
4098 * Purpose:
4099 * Create a KnownUnit structure describing an alias for a known unit.
4100
4101 * Type:
4102 * Private function.
4103
4104 * Synopsis:
4105 * #include "unit.h"
4106 * void MakeUnitAlias( const char *sym, const char *alias, int *status )
4107
4108 * Class Membership:
4109 * Unit member function.
4110
4111 * Description:
4112 * This function creates a KnownUnit structure decribing an alias for a
4113 * known unit, and adds it to the head of the linked list of known units
4114 * stored in a module variable. An alias is a KnownUnit which is
4115 * identical to an existing known but which has a different symbol.
4116
4117 * Parameters:
4118 * sym
4119 * A pointer to the symbol string of an existing KnwonUnit. The string
4120 * should not include any multiplier prefix.
4121 * alias
4122 * A pointer to the symbol string to use as the alasi for the existing
4123 * KnownUnit. The string should not include any multiplier prefix.
4124 * status
4125 * Pointer to the inherited status variable.
4126
4127 * Notes:
4128 * - The supplied symbol and label strings are not copied. The
4129 * supplied pointers are simply stored in the returned structure.
4130 * Therefore the strings to which the pointers point should not be
4131 * modified after this function returned (in fact this function is
4132 * always called with literal strings for these arguments).
4133 */
4134
4135 /* Local Variables: */
4136 KnownUnit *unit;
4137
4138 /* Check the global error status. */
4139 if( !astOK ) return;
4140
4141 /* Search the existing list of KnownUnits for the specified symbol. */
4142 unit = known_units;
4143 while( unit ) {
4144 if( !strcmp( sym, unit->sym ) ) {
4145
4146 /* Create a new KnownUnit for the alias. It will becomes the head of the
4147 known units chain. */
4148 MakeKnownUnit( alias, unit->label, NULL, status );
4149
4150 /* Store a pointer to the KnownUnit which is to be used in place of the
4151 alias. */
4152 known_units->use = unit;
4153
4154 /* Leave the loop. */
4155 break;
4156 }
4157
4158 /* Move on to check the next existing KnownUnit. */
4159 unit = unit->next;
4160 }
4161
4162 /* Report an error if the supplied unit was not found. */
4163 if( !unit ) {
4164 astError( AST__INTER, "MakeUnitAlias(Unit): Cannot find existing "
4165 "units \"%s\" to associate with the alias \"%s\" (AST "
4166 "internal programming error).", status, sym, alias );
4167 }
4168 }
4169
ModifyPrefix(UnitNode * old,int * status)4170 static UnitNode *ModifyPrefix( UnitNode *old, int *status ) {
4171 /*
4172 * Name:
4173 * ModifyPrefix
4174
4175 * Purpose:
4176 * Replace a MULT or DIV node with a LDVAR and suitable multiplier.
4177
4178 * Type:
4179 * Private function.
4180
4181 * Synopsis:
4182 * #include "unit.h"
4183 * UnitNode *ModifyPrefix( UnitNode *old, int *status )
4184
4185 * Class Membership:
4186 * Unit member function.
4187
4188 * Description:
4189 * This function checks the supplied node. If it is a DIV or MULT node
4190 * in which one argument is an LDVAR and the other is a constant, then
4191 * its checks to see if the constant can be absorbed into the LDVAR by
4192 * changing the multiplier in the LDVAR node. If so, it returns a new
4193 * node which is an LDVAR with the modified multiplier. Otherwise it
4194 * returns NULL.
4195
4196 * Parameters:
4197 * old
4198 * Pointer to an existing UnitNode to be checked.
4199 * status
4200 * Pointer to the inherited status variable.
4201
4202 * Returned Value:
4203 * A pointer to the new UnitNode.
4204
4205 * Notes:
4206 * - A value of NULL will be returned if this function is invoked with
4207 * the global error status set, or if it should fail for any reason.
4208 */
4209
4210 /* Local Variables: */
4211 Multiplier *mult;
4212 Multiplier *mmult;
4213 UnitNode *ldcon;
4214 UnitNode *ldvar;
4215 UnitNode *newtree;
4216 UnitNode *result;
4217 double con;
4218 double cmult;
4219 double r;
4220 double rmin;
4221 int recip;
4222 int changed;
4223
4224 /* Initialise. */
4225 result = NULL;
4226
4227 /* Check the inherited status. */
4228 if( !astOK ) return result;
4229
4230 /* Indicate that we have not yet found any reason to return a changed
4231 node. */
4232 changed = 0;
4233
4234 /* Check the supplied node is a DIV or MULT node. */
4235 if( old->opcode == OP_DIV || old->opcode == OP_MULT ) {
4236
4237 /* Get a copy of the supplied tree which we can modify safely. */
4238 newtree = CopyTree( old, status );
4239
4240 /* Identify the LDVAR argument (if any). */
4241 if( newtree->arg[ 0 ]->opcode == OP_LDVAR ) {
4242 ldvar = newtree->arg[ 0 ];
4243
4244 } else if( newtree->arg[ 1 ]->opcode == OP_LDVAR ) {
4245 ldvar = newtree->arg[ 1 ];
4246
4247 } else {
4248 ldvar = NULL;
4249 }
4250
4251 /* Identify the LDCON argument (if any). */
4252 if( newtree->arg[ 0 ]->opcode == OP_LDCON ) {
4253 ldcon = newtree->arg[ 0 ];
4254
4255 } else if( newtree->arg[ 1 ]->opcode == OP_LDCON ) {
4256 ldcon = newtree->arg[ 1 ];
4257
4258 } else {
4259 ldcon = NULL;
4260 }
4261
4262 /* If either was not found, return NULL. */
4263 if( !ldvar || !ldcon ) {
4264 newtree = FreeTree( newtree, status );
4265
4266 /* Otherwise, extract the multiplier constant. If there is no multiplier, the
4267 constant is 1.0. */
4268 } else {
4269 cmult = ldvar->mult ? ldvar->mult->scale: 1.0;
4270
4271 /* Extract the constant. */
4272 con = ldcon->con;
4273
4274 /* Combine the multiplier and the constant. The resulting constant is a
4275 factor which is used to multiply the LDVAR quantity. If the original
4276 node is a DIV node in which the LDVAR is in the denominator, then
4277 flag that we need to reciprocate the new MULT node which represents
4278 "constant*LDVAR" before returning. */
4279 if( newtree->opcode == OP_MULT ) {
4280 con = con*cmult;
4281 recip = 0;
4282 } else {
4283 con = cmult/con;
4284 recip = ( ldvar == newtree->arg[ 1 ] );
4285 }
4286
4287 /* Find the closest known multiplier to the new constant. */
4288 rmin = ( con > 1 ) ? con : 1.0/con;
4289 mmult = NULL;
4290 mult = GetMultipliers( status );
4291 while( mult ) {
4292 r = ( con > mult->scale) ? con/mult->scale : mult->scale/con;
4293 if( r < rmin ) {
4294 mmult = mult;
4295 rmin = r;
4296 }
4297 mult = mult->next;
4298 }
4299
4300 /* Modify the constant to take account of the new multiplier chosen
4301 above. "mmult" will be NULL if the best multiplier is unity. */
4302 if( mmult ) con = con/mmult->scale;
4303
4304 /* If they have changed, associate the chosen multiplier with the LDVAR node,
4305 and the constant with the LDCON node. */
4306 if( ldvar->mult != mmult ) {
4307 ldvar->mult = mmult;
4308 changed = 1;
4309 }
4310
4311 if( ldcon->con != con ) {
4312 ldcon->con = con;
4313 changed = 1;
4314 }
4315
4316 /* Unless the node is proportional to the reciprocal of the variable, the
4317 new node should be a MULT node (it may originally have been a DIV). */
4318 if( !recip ) {
4319 if( newtree->opcode != OP_MULT ){
4320 newtree->opcode = OP_MULT;
4321 changed = 1;
4322 }
4323
4324 /* If the constant is 1.0 we can just return the LDVAR node by itself. */
4325 if( fabs( con - 1.0 ) < 1.0E-6 ) {
4326 result = CopyTree( ldvar, status );
4327 newtree = FreeTree( newtree, status );
4328 changed = 1;
4329
4330 /* Otherwise return the modified tree containing both LDVAR and LDCON nodes. */
4331 } else {
4332 result = newtree;
4333 }
4334
4335 /* If the node is proportional to the reciprocal of the variable, the
4336 new node will already be a DIV node and will have an LDCON as the first
4337 argument (numerator) and an LDVAR as the second argument (denominator). */
4338 } else {
4339
4340 /* The first argument (the numerator) should be the reciprocal of the constant
4341 found above. */
4342 ldcon->con = 1.0/ldcon->con;
4343 if( !EQUAL( ldcon->con, old->arg[0]->con ) ) changed = 1;
4344
4345 /* Return the modified tree containing both LDVAR and LDCON nodes. */
4346 result = newtree;
4347 }
4348 }
4349 }
4350
4351 /* If the new and old trees are equivalent, then we do not need to return
4352 it. */
4353 if( !changed && result ) result = FreeTree( result, status );
4354
4355 /* Return the answer. */
4356 return result;
4357 }
4358
4359
NewNode(UnitNode * old,Oper code,int * status)4360 static UnitNode *NewNode( UnitNode *old, Oper code, int *status ) {
4361 /*
4362 * Name:
4363 * NewNode
4364
4365 * Purpose:
4366 * Create and initialise a new UnitNode.
4367
4368 * Type:
4369 * Private function.
4370
4371 * Synopsis:
4372 * #include "unit.h"
4373 * UnitNode *NewNode( UnitNode *old, Oper code, int *status )
4374
4375 * Class Membership:
4376 * Unit member function.
4377
4378 * Description:
4379 * This function creates and initialises a new UnitNode, or
4380 * re-initialises an existing UnitNode to use a different op code.
4381
4382 * Parameters:
4383 * old
4384 * Pointer to an existing UnitNode to be modified, or NULL to create
4385 * a new UnitNode.
4386 * code
4387 * The op code for the new UnitNode.
4388 * status
4389 * Pointer to the inherited status variable.
4390
4391 * Returned Value:
4392 * A pointer to the new UnitNode.
4393
4394 * Notes:
4395 * - A value of NULL will be returned if this function is invoked with
4396 * the global error status set, or if it should fail for any reason.
4397 */
4398
4399 /* Local Variables: */
4400 UnitNode **args;
4401 UnitNode *result;
4402 int i;
4403
4404 /* Initialise. */
4405 result = NULL;
4406 args = NULL;
4407
4408 /* Check the inherited status. */
4409 if( !astOK ) return result;
4410
4411 /* If an existig UnitNode was supplied, free any memory used to hold
4412 pointers to its arguments. */
4413 if( old ) {
4414 old->arg = astFree( old->arg );
4415 result = old;
4416
4417 /* Otherwise, allocate memory for a new structure. */
4418 } else {
4419 result = astMalloc( sizeof( UnitNode ) );
4420 }
4421
4422 /* Check the pointer can be used safely. */
4423 if( astOK ) {
4424
4425 /* Initialise the members of the UnitNode structure. */
4426 result->opcode = code;
4427 result->arg = NULL;
4428 result->con = AST__BAD;
4429 result->name = NULL;
4430 result->unit = NULL;
4431 result->mult = NULL;
4432 result->narg = 0;
4433
4434 switch( code ){
4435 case OP_LDPI:
4436 result->con = PI;
4437 break;
4438
4439 case OP_LDE:
4440 result->con = E;
4441 break;
4442
4443 case OP_LOG:
4444 case OP_LN:
4445 case OP_EXP:
4446 case OP_SQRT:
4447 result->narg = 1;
4448 break;
4449
4450 case OP_POW:
4451 case OP_DIV:
4452 case OP_MULT:
4453 result->narg = 2;
4454 break;
4455
4456 default:
4457 ;
4458 }
4459
4460 /* Allocate memory for the UnitNode pointers which will locate the
4461 nodes forming the arguments to the new node. */
4462 args = astMalloc( (result->narg)*sizeof( UnitNode * ) );
4463 if( astOK ) {
4464 result->arg = args;
4465
4466 /* Initialise the argument pointers to NULL. */
4467 for( i = 0; i < result->narg; i++ ) args[ i ] = NULL;
4468 }
4469 }
4470
4471 /* Free any result if an error occurred. */
4472 if( !astOK ) {
4473 args = astFree( args );
4474 result = astFree( result );
4475 }
4476
4477 /* Return the answer. */
4478 return result;
4479 }
4480
RemakeTree(UnitNode ** node,int * status)4481 static void RemakeTree( UnitNode **node, int *status ) {
4482 /*
4483 * Name:
4484 * RemakeTree
4485
4486 * Purpose:
4487 * Replace derived units within a tree of UnitNodes by basic units.
4488
4489 * Type:
4490 * Private function.
4491
4492 * Synopsis:
4493 * #include "unit.h"
4494 * void RemakeTree( UnitNode **node, int *status )
4495
4496 * Class Membership:
4497 * Unit member function.
4498
4499 * Description:
4500 * This function searches for LDVAR nodes (i.e. references to unit
4501 * symbols) within the given tree, and replaces each such node which
4502 * refers to known derived unit with a sub-tree of nodes which
4503 * define the derived unit in terms of known basic units.
4504
4505 * Parameters:
4506 * node
4507 * The address of a pointer to the UnitNode at the head of the tree
4508 * which is to be simplified. On exit the supplied tree is freed and a
4509 * pointer to a new tree is placed at the given address.
4510 * status
4511 * Pointer to the inherited status variable.
4512
4513 */
4514
4515 /* Local Variables: */
4516 KnownUnit *unit;
4517 int i;
4518 UnitNode *newnode;
4519
4520 /* Check inherited status. */
4521 if( !astOK ) return;
4522
4523 /* Initially, we have no replacement node */
4524 newnode = NULL;
4525
4526 /* If this is an LDVAR node... */
4527 if( (*node)->opcode == OP_LDVAR ) {
4528
4529 /* If the LDVAR node has a multiplier associated with it, we need to
4530 introduce a OP_MULT node to perform the scaling. */
4531 if( (*node)->mult ) {
4532 newnode = NewNode( NULL, OP_MULT, status );
4533 if( astOK ) {
4534 newnode->arg[0] = NewNode( NULL, OP_LDCON, status );
4535 if( astOK ) {
4536 newnode->arg[0]->con = 1.0/(*node)->mult->scale;
4537
4538 /* See if the node refers to a known unit. If not, or if the known unit
4539 is a basic unit (i.e. not a derived unit) use the supplied node for
4540 the second argument of the OP_MULT node (without the multiplier).
4541 Otherwise, use a copy of the tree which defines the derived unit. */
4542 unit = (*node)->unit;
4543 if( unit && unit->head ) {
4544 newnode->arg[1] = CopyTree( unit->head, status );
4545 } else {
4546 newnode->arg[1] = CopyTree( *node, status );
4547 if( astOK ) newnode->arg[1]->mult = NULL;
4548 }
4549 }
4550 }
4551
4552 /* If no multiplier is supplied, the replacement node is simply the tree
4553 which defines the unscaled unit (if known), or the original node (if
4554 unknown). */
4555 } else {
4556 unit = (*node)->unit;
4557 if( unit && unit->head ) newnode = CopyTree( unit->head, status );
4558 }
4559
4560 /* If this is not an LDVAR Node, remake the sub-trees which form the
4561 arguments of this node. */
4562 } else {
4563 for( i = 0; i < (*node)->narg; i++ ) {
4564 RemakeTree( &((*node)->arg[ i ]), status );
4565 }
4566 }
4567
4568 /* If an error has occurred, free any new node. */
4569 if( !astOK ) newnode = FreeTree( newnode, status );
4570
4571 /* If we have a replacement node, free the supplied tree and return a
4572 pointer to the new tree. */
4573 if( newnode ) {
4574 FreeTree( *node, status );
4575 *node = newnode;
4576 }
4577 }
4578
ReplaceNode(UnitNode * target,UnitNode * old,UnitNode * new,int * status)4579 static int ReplaceNode( UnitNode *target, UnitNode *old, UnitNode *new, int *status ) {
4580 /*
4581 * Name:
4582 * ReplaceNode
4583
4584 * Purpose:
4585 * Replace a node within a tree of UnitNodes.
4586
4587 * Type:
4588 * Private function.
4589
4590 * Synopsis:
4591 * #include "unit.h"
4592 * int ReplaceNode( UnitNode *target, UnitNode *old, UnitNode *new, int *status )
4593
4594 * Class Membership:
4595 * Unit member function.
4596
4597 * Description:
4598 * This function replaces a specified node within a tree of UnitNodes
4599 * with another given node. The original node is freed if found.
4600
4601 * Parameters:
4602 * target
4603 * A pointer to the UnitNode at the head of the tree containing the
4604 * node to be replaced.
4605 * old
4606 * A pointer to the UnitNode to be replaced.
4607 * new
4608 * A pointer to the UnitNode to replace "old".
4609 * status
4610 * Pointer to the inherited status variable.
4611
4612 * Return Value:
4613 * Non-zero if the "old" node was found and replaced (in which case
4614 * the "old" node will have been freed).
4615
4616 * Notes:
4617 * - It is assumed that the "old" node occurs at most once within the
4618 * target tree.
4619 * - The node at the head of the target tree is not compared with the
4620 * "old" node. It is assumed the called will already have done this.
4621 * - A value of zero is returned if an error has already occurred, or
4622 * if this function fails for any reason.
4623
4624 */
4625
4626 /* Local Variables: */
4627 int i;
4628 int result;
4629
4630 /* Initialise */
4631 result = 0;
4632
4633 /* Check inherited status. */
4634 if( !astOK ) return result;
4635
4636 /* Loop round the arguments of the node at the head of the target tree.
4637 Break out of the loop as soone as the old node is found. */
4638 for( i = 0; i < target->narg; i++ ) {
4639
4640 /* If this argument is the node to be replaced, free the old one and store
4641 the new one, and then leave the loop. */
4642 if( target->arg[ i ] == old ) {
4643 FreeTree( old, status );
4644 target->arg[ i ] = new;
4645 result = 1;
4646 break;
4647
4648 /* Otherwise use this function recursively to search for the old node
4649 within the current argument. */
4650 } else {
4651 if( ReplaceNode( target->arg[ i ], old, new, status ) ) break;
4652 }
4653 }
4654
4655 /* If an error has occurred, return zero. */
4656 if( !astOK ) result = 0;
4657
4658 /* Return the answer. */
4659 return result;
4660 }
4661
SimplifyTree(UnitNode ** node,int std,int * status)4662 static int SimplifyTree( UnitNode **node, int std, int *status ) {
4663 /*
4664 * Name:
4665 * SimplifyTree
4666
4667 * Purpose:
4668 * Simplify a tree of UnitNodes.
4669
4670 * Type:
4671 * Private function.
4672
4673 * Synopsis:
4674 * #include "unit.h"
4675 * int SimplifyTree( UnitNode **node, int std, int *status )
4676
4677 * Class Membership:
4678 * Unit member function.
4679
4680 * Description:
4681 * This function simplifies a tree of UnitNodes. It is assumed that
4682 * all the OP_LDVAR nodes in the tree refer to the same basic unit.
4683 * A primary purpose of this function is to standardise the tree so
4684 * that trees which implement equivalent transformations but which
4685 * have different structures can be compared (for instance, so that
4686 * "2*x" and "x*2" are treated as equal trees). If "std" is non-zero,
4687 * reducing the complexity of the tree is only of secondary importance.
4688 * This explains why some "simplifications" actually produced trees which
4689 * are more complicated.
4690
4691 * Parameters:
4692 * node
4693 * The address of a pointer to the UnitNode at the head of the tree
4694 * which is to be simplified. On exit the supplied tree is freed and a
4695 * pointer to a new tree is placed at the given address.
4696 * std
4697 * If non-zero, perform standardisations. Otherwise only perform
4698 * genuine simplifications.
4699 * status
4700 * Pointer to the inherited status variable.
4701
4702 * Returned Value:
4703 * Non-zero if some change was made to the tree.
4704
4705 */
4706
4707 /* Local Variables: */
4708 int i;
4709 UnitNode *newnode;
4710 UnitNode *node1;
4711 UnitNode *node2;
4712 Oper op;
4713 UnitNode **factors;
4714 double *powers;
4715 int nfactor;
4716 double coeff;
4717 int result;
4718
4719 /* Initialise */
4720 result = 0;
4721
4722 /* Check inherited status. */
4723 if( !astOK ) return result;
4724
4725 /* Initiallially, we have no replacement node. */
4726 newnode = NULL;
4727
4728 /* First replace any complex constant expressions any corresponding
4729 OP_LDCON nodes. */
4730 FixConstants( node, 0, status );
4731
4732 /* Simplify the sub-trees corresponding to the arguments of the node at
4733 the head of the supplied tree. */
4734 for( i = 0; i < (*node)->narg; i++ ) {
4735 if( SimplifyTree( &( (*node)->arg[ i ] ), std, status ) ) result = 1;
4736 }
4737
4738 /* Now do specific simplifications appropriate to the nature of the node at
4739 the head of the tree. */
4740 op = (*node)->opcode;
4741
4742 /* Natural log */
4743 /* =========== */
4744 /* We standardise argument powers into coefficients of the LN value. */
4745 if( op == OP_LN ) {
4746
4747 /* If the argument is a OP_EXP node, they cancel out. Return a copy of the
4748 argument of OP_EXP node. */
4749 if( (*node)->arg[ 0 ]->opcode == OP_EXP ) {
4750 newnode = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
4751
4752 /* If the argument is an OP_POW node, rearrange the nodes to represent
4753 k*ln(x) instead of ln(x**k) (note pow nodes always have a constant
4754 exponent - this is checked in InvertConstants). SQRT arguments will
4755 not occur because they will have been changed into POW nodes when the
4756 arguments of the supplied head node were simplified above. */
4757 } else if( std && (*node)->arg[ 0 ]->opcode == OP_POW ) {
4758 newnode = NewNode( NULL, OP_MULT, status );
4759 node1 = CopyTree( (*node)->arg[ 0 ]->arg[ 1 ], status );
4760 node2 = NewNode( NULL, OP_LN, status );
4761 if( astOK ) {
4762 node2->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
4763 newnode->arg[ 0 ] = node1;
4764 newnode->arg[ 1 ] = node2;
4765 }
4766 }
4767
4768 /* Common log */
4769 /* ========== */
4770 /* We standardise natural logs into common logs. */
4771 } else if( op == OP_LOG ) {
4772 if( std ) {
4773 newnode = NewNode( NULL, OP_DIV, status );
4774 node1 = NewNode( NULL, OP_LN, status );
4775 node2 = NewNode( NULL, OP_LDCON, status );
4776 if( astOK ) {
4777 node1->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ], status );
4778 node2->con = log( 10.0 );
4779 newnode->arg[ 0 ] = node1;
4780 newnode->arg[ 1 ] = node2;
4781 }
4782 }
4783
4784 /* Exponential */
4785 /* =========== */
4786 /* We prefer to minimise the number of EXP nodes, so, for instance, we do not
4787 change "exp(x*y)" to "exp(x)+exp(y)" (and the code for ADD nodes does
4788 the inverse conversion). */
4789 } else if( op == OP_EXP ) {
4790
4791 /* If the argument is an OP_LN node, they cancel out. Return a copy of the
4792 argument of the OP_LN node. Common log arguments will not occur because
4793 they will have been changed into natural logs when the arguments of
4794 the supplied head node were simplified above. */
4795 if( (*node)->arg[ 0 ]->opcode == OP_LN ) {
4796 newnode = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
4797 }
4798
4799 /* Square root */
4800 /* =========== */
4801 /* We standardise sqrt nodes into pow nodes. */
4802 } else if( op == OP_SQRT ) {
4803 if( std ) {
4804 newnode = NewNode( NULL, OP_POW, status );
4805 node1 = CopyTree( (*node)->arg[ 0 ], status );
4806 node2 = NewNode( NULL, OP_LDCON, status );
4807 if( astOK ) {
4808 node2->con = 0.5;
4809 newnode->arg[ 0 ] = node1;
4810 newnode->arg[ 1 ] = node2;
4811 }
4812 }
4813
4814 /* Exponentiation */
4815 /* ============== */
4816 /* We want to simplfy factors. So, for instance, (x*y)**k is converted to
4817 (x**k)*(y**k). */
4818 } else if( op == OP_POW ) {
4819
4820 /* If the first argument is an OP_EXP node, then change "(e**x)**k" into
4821 "e**(k*x)" */
4822 if( (*node)->arg[ 0 ]->opcode == OP_EXP ) {
4823 newnode = NewNode( NULL, OP_EXP, status );
4824 node1 = NewNode( NULL, OP_MULT, status );
4825 if( astOK ) {
4826 node1->arg[ 0 ] = CopyTree( (*node)->arg[ 1 ], status );
4827 node1->arg[ 1 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
4828 newnode->arg[ 0 ] = node1;
4829 }
4830
4831 /* "x**0" can be replaced by 1.0 */
4832 } else if( (*node)->arg[ 1 ]->con == 0.0 ) {
4833 newnode = NewNode( NULL, OP_LDCON, status );
4834 if( astOK ) newnode->con = 1.0;
4835
4836 /* "x**1" can be replaced by x */
4837 } else if( EQUAL( (*node)->arg[ 1 ]->con, 1.0 ) ) {
4838 newnode = CopyTree( (*node)->arg[ 0 ], status );
4839
4840 /* If the first argument is an OP_POW node, then change "(x**k1)**k2" into
4841 "x**(k1*k2)" */
4842 } else if( (*node)->arg[ 0 ]->opcode == OP_POW ) {
4843 newnode = NewNode( NULL, OP_POW, status );
4844 node1 = NewNode( NULL, OP_LDCON, status );
4845 if( astOK ) {
4846 node1->con = ( (*node)->arg[ 0 ]->arg[ 1 ]->con )*
4847 ( (*node)->arg[ 1 ]->con );
4848 newnode->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
4849 newnode->arg[ 1 ] = node1;
4850 }
4851
4852 /* If the first argument is an OP_MULT node, then change "(x*y)**k" into
4853 "(x**(k))*(y**(k))" */
4854 } else if( std && (*node)->arg[ 0 ]->opcode == OP_MULT ) {
4855 newnode = NewNode( NULL, OP_MULT, status );
4856 node1 = NewNode( NULL, OP_POW, status );
4857 if( astOK ) {
4858 node1->arg[ 1 ] = CopyTree( (*node)->arg[ 1 ], status );
4859 node2 = CopyTree( node1, status );
4860 node1->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 0 ], status );
4861 node2->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ]->arg[ 1 ], status );
4862 newnode->arg[ 0 ] = node1;
4863 newnode->arg[ 1 ] = node2;
4864 }
4865 }
4866
4867 /* Division. */
4868 /* ========= */
4869 /* We standardise divisions into corresponding multiplications. */
4870 } else if( op == OP_DIV ) {
4871
4872 /* Division by 1 is removed. */
4873 if( EQUAL( (*node)->arg[ 1 ]->con, 1.0 ) ){
4874 newnode = CopyTree( (*node)->arg[ 0 ], status );
4875
4876 /* Division by any other constant (except zero) is turned into a
4877 multiplication by the reciprocal constant. */
4878 } else if( (*node)->arg[ 1 ]->con != AST__BAD ) {
4879 if( (*node)->arg[ 1 ]->con != 0.0 ) {
4880 newnode = NewNode( NULL, OP_MULT, status );
4881 node1 = NewNode( NULL, OP_LDCON, status );
4882 if( astOK ) {
4883 node1->con = 1.0/(*node)->arg[ 1 ]->con;
4884 newnode->arg[ 0 ] = node1;
4885 newnode->arg[ 1 ] = CopyTree( (*node)->arg[ 0 ], status );
4886 }
4887 } else {
4888 astError( AST__BADUN, "Simplifying a units expression"
4889 "requires a division by zero." , status);
4890 }
4891
4892 /* Other divisions "x/y" are turned into "x*(y**(-1))" */
4893 } else if( std ) {
4894 newnode = NewNode( NULL, OP_MULT, status );
4895 node1 = NewNode( NULL, OP_POW, status );
4896 node2 = NewNode( NULL, OP_LDCON, status );
4897 if( astOK ) {
4898 node2->con = -1.0;
4899 node1->arg[ 0 ] = CopyTree( (*node)->arg[ 1 ], status );
4900 node1->arg[ 1 ] = node2;
4901 newnode->arg[ 0 ] = CopyTree( (*node)->arg[ 0 ], status );
4902 newnode->arg[ 1 ] = node1;
4903 }
4904 }
4905
4906 /* Multiplication */
4907 /* ============== */
4908 } else if( op == OP_MULT ) {
4909
4910 /* If the right hand argument is constant, swap the arguments. */
4911 if( (*node)->arg[ 1 ]->con != AST__BAD ) {
4912 newnode = NewNode( NULL, OP_MULT, status );
4913 if( astOK ) {
4914 newnode->arg[ 0 ] = CopyTree( (*node)->arg[ 1 ], status );
4915 newnode->arg[ 1 ] = CopyTree( (*node)->arg[ 0 ], status );
4916 }
4917
4918 /* Multiplication by zero produces a constant zero. */
4919 } else if( (*node)->arg[ 0 ]->con == 0.0 ){
4920 newnode = NewNode( NULL, OP_LDCON, status );
4921 if( astOK ) newnode->con = 0.0;
4922
4923 /* Multiplication by 1 is removed. */
4924 } else if( EQUAL( (*node)->arg[ 0 ]->con, 1.0 ) ){
4925 newnode = CopyTree( (*node)->arg[ 1 ], status );
4926
4927 /* For other MULT nodes, analyse the tree to find a list of all its
4928 factors with an associated power for each one, and an overall constant
4929 coefficient. */
4930 } else if( std ) {
4931 FindFactors( (*node), &factors, &powers, &nfactor, &coeff, status );
4932
4933 /* Produce a new tree from these factors. The factors are standardised by
4934 ordering them alphabetically (after conversion to a character string). */
4935 newnode = CombineFactors( factors, powers, nfactor, coeff, status );
4936
4937 /* Free resources */
4938 factors = astFree( factors );
4939 powers = astFree( powers );
4940
4941 }
4942 }
4943
4944 /* If we have produced a new node which is identical to the old node,
4945 free it. Otherwise, indicate we have made some changes. */
4946 if( newnode ) {
4947 if( !CmpTree( newnode, *node, 1, status ) ) {
4948 newnode = FreeTree( newnode, status );
4949 } else {
4950 result = 1;
4951 }
4952 }
4953
4954 /* If an error has occurred, free any new node. */
4955 if( !astOK ) newnode = FreeTree( newnode, status );
4956
4957 /* If we have a replacement node, free the supplied tree and return a
4958 pointer to the new tree. */
4959 if( newnode ) {
4960 FreeTree( *node, status );
4961 *node = newnode;
4962 }
4963
4964 /* If the above produced some change, try re-simplifying the tree. */
4965 if( result ) SimplifyTree( node, std, status );
4966
4967 /* Return the result. */
4968 return result;
4969
4970 }
4971
SplitUnit(const char * str,int ls,const char * u,int cs,Multiplier ** mult,int * l,int * status)4972 static int SplitUnit( const char *str, int ls, const char *u, int cs,
4973 Multiplier **mult, int *l, int *status ) {
4974 /*
4975 * Name:
4976 * SplitUnit
4977
4978 * Purpose:
4979 * Split a given string into unit name and multiplier.
4980
4981 * Type:
4982 * Private function.
4983
4984 * Synopsis:
4985 * #include "unit.h"
4986 * int SplitUnit( const char *str, int ls, const char *u, int cs,
4987 * Multiplier **mult, int *l, int *status )
4988
4989 * Class Membership:
4990 * Unit member function.
4991
4992 * Description:
4993 * Returns non-zer0 if the supplied string ends with the supplied unit
4994 * name or label, and any leading string is a known multiplier.
4995
4996 * Parameters:
4997 * str
4998 * The string to test, typically containing a multiplier and a unit
4999 * symbol or label.
5000 * ls
5001 * Number of characters to use from "str" (not including trailing null)
5002 * u
5003 * Pointer to the unit label or symbol string to be searched for.
5004 * cs
5005 * If non-zero, the test for "u" is case insensitive.
5006 * mult
5007 * Address of a location at which to return the multiplier at the
5008 * start of the supplied string. NULL is returned if the supplied
5009 * string does not match the supplied unit, or if the string
5010 * includes no multiplier.
5011 * l
5012 * Address of an int in which to return the length of "u".
5013 * status
5014 * Pointer to the inherited status variable.
5015
5016 * Returned Value:
5017 * Non-zero if "str" ends with "u" and starts with a null string or a
5018 * known multiplier string.
5019
5020 */
5021
5022 /* Local Variables: */
5023 int ret;
5024 int lm;
5025 int lu;
5026
5027 /* Initialise */
5028 ret = 0;
5029 *mult = NULL;
5030 *l = 0;
5031
5032 /* Check inherited status. */
5033 if( !astOK ) return ret;
5034
5035 /* Find the number of characters in the supplied unit label or symbol. The
5036 difference between lu and ls must be the length of the multiplier. */
5037 lu = strlen( u );
5038 lm = ls - lu;
5039
5040 /* Make sure "str" is not shorter than "u" */
5041 if( lm >= 0 ) {
5042
5043 /* Compare the end of "str" against "u" */
5044 if( cs ? !strncmp( str + lm, u, lu ) :
5045 !Ustrncmp( str + lm, u, lu, status ) ) {
5046 ret = 1;
5047
5048 /* If "str" ends with "u", see if it starts with a known multiplier */
5049 if( lm > 0 ) {
5050 ret = 0;
5051 *mult = GetMultipliers( status );
5052 while( *mult ) {
5053 if( (*mult)->symlen == lm && !strncmp( str, (*mult)->sym, lm ) ) {
5054 ret = 1;
5055 break;
5056 }
5057 *mult = (*mult)->next;
5058 }
5059
5060 /* If not, try again using case-insensitive matching. */
5061 if( !ret ) {
5062 *mult = GetMultipliers( status );
5063 while( *mult ) {
5064 if( (*mult)->symlen == lm && !Ustrncmp( str, (*mult)->sym, lm, status ) ) {
5065 ret = 1;
5066 break;
5067 }
5068 *mult = (*mult)->next;
5069 }
5070 }
5071
5072 /* If not, try again using case-insensitive matching against the
5073 multiplier label. */
5074 if( !ret ) {
5075 *mult = GetMultipliers( status );
5076 while( *mult ) {
5077 if( (*mult)->lablen == lm && !Ustrncmp( str, (*mult)->label, lm, status ) ) {
5078 ret = 1;
5079 break;
5080 }
5081 *mult = (*mult)->next;
5082 }
5083 }
5084
5085 }
5086 }
5087 }
5088
5089 *l = lu;
5090 return ret;
5091 }
5092
astUnitAnalyser_(const char * in,double powers[9],int * status)5093 double astUnitAnalyser_( const char *in, double powers[9], int *status ){
5094 /*
5095 *+
5096 * Name:
5097 * astUnitAnalyser
5098
5099 * Purpose:
5100 * Perform a dimensional analysis of a unti string.
5101
5102 * Type:
5103 * Protected function.
5104
5105 * Synopsis:
5106 * #include "unit.h"
5107 * double astUnitAnalyser_( const char *in, double powers[9] )
5108
5109 * Class Membership:
5110 * Unit member function.
5111
5112 * Description:
5113 * This function parses the supplied units string if possible, and
5114 * returns a set of pwoers and a scaling factor which represent the
5115 * units string.
5116
5117 * Parameters:
5118 * in
5119 * A string representation of the units, for instance "km/h".
5120 * powers
5121 * An array in which are returned the powers for each of the following
5122 * basic units (in the order shown): kilogramme, metre, second, radian,
5123 * Kelvin, count, adu, photon, magnitude, pixel. If the supplied unit
5124 * does not depend on a given basic unit a value of 0.0 will be returned
5125 * in the array. The returns values represent a system of units which is
5126 * a scaled form of the supplied units, expressed in the basic units of
5127 * m, kg, s, rad, K, count, adu, photon, mag and pixel. For instance, a
5128 * returned array of [1,0,-2,0,0,0,0,0,0] would represent "m/s**2".
5129
5130 * Returned Value:
5131 * A scaling factor for the supplied units. The is the value, in the
5132 * units represented by the returned powers, which corresponds to a
5133 * value of 1.0 in the supplied units.
5134
5135 * Notes:
5136 * - An error will be reported if the units string cannot be parsed
5137 * or refers to units or functions which cannot be analysed in this way.
5138 * - AST__BAD is returned if this function is invoked with the
5139 * global error status set or if it should fail for any reason.
5140 *-
5141 */
5142
5143 /* Local Variables: */
5144 UnitNode *in_tree;
5145 double result;
5146
5147 /* Initialise */
5148 result = AST__BAD;
5149
5150 /* Check the global error status. */
5151 if ( !astOK ) return result;
5152
5153 /* Parse the input units string, producing a tree of UnitNodes which
5154 represents the input units. A pointer to the UnitNode at the head of
5155 the tree is returned if succesfull. Report a context message if this
5156 fails. */
5157 in_tree = CreateTree( in, 1, 1, status );
5158 if( in_tree ) {
5159
5160 /* Analyse the tree */
5161 if( !DimAnal( in_tree, powers, &result, status ) && astOK ) {
5162 result = AST__BAD;
5163 astError( AST__BADUN, "astUnitAnalyser: Error analysing input "
5164 "units string '%s' (it may contain unsupported "
5165 "functions or dimensionless units).", status, in );
5166 }
5167
5168 /* Free the tree. */
5169 in_tree = FreeTree( in_tree, status );
5170
5171 } else if( astOK ) {
5172 astError( AST__BADUN, "astUnitAnalyser: Error parsing input "
5173 "units string '%s'.", status, in );
5174 }
5175
5176 /* Return the result */
5177 return result;
5178 }
5179
astUnitLabel_(const char * sym,int * status)5180 const char *astUnitLabel_( const char *sym, int *status ){
5181 /*
5182 *+
5183 * Name:
5184 * astUnitLabel
5185
5186 * Purpose:
5187 * Return a string label for a given unit symbol.
5188
5189 * Type:
5190 * Protected function.
5191
5192 * Synopsis:
5193 * #include "unit.h"
5194 * const char *astUnitLabel( const char *sym )
5195
5196 * Class Membership:
5197 * Unit member function.
5198
5199 * Description:
5200 * This function returns a pointer to a constant string containing a
5201 * descriptive label for the unit specified by the given unit symbol.
5202
5203 * Parameters:
5204 * sym
5205 * A string holing a known unit symbol.
5206
5207 * Returned Value:
5208 * A pointer to constant string holding a descriptive label for the
5209 * supplied unit. A NULL pointer is returned (without error) if the
5210 * supplied unit is unknown.
5211
5212 * Notes:
5213 * - A NULL pointer is returned if this function is invoked with the
5214 * global error status set or if it should fail for any reason.
5215 *-
5216 */
5217
5218 /* Local Variables: */
5219 const char *result;
5220 KnownUnit *unit;
5221
5222 /* Initialise */
5223 result = NULL;
5224
5225 /* Check the global error status. */
5226 if ( !astOK ) return result;
5227
5228 /* Ensure descriptions of the known units are available. */
5229 unit = GetKnownUnits( 1, status );
5230
5231 /* Loop through the chain of known units looking for a unit with a symbol
5232 equal to the supplied string. If found, store a pointer to its label
5233 and break out of the loop. */
5234 while( unit ) {
5235 if( !strcmp( sym, unit->sym ) ) {
5236 result = unit->label;
5237 break;
5238 }
5239
5240 unit = unit->next;
5241 }
5242
5243 /* Return the answer. */
5244 return result;
5245 }
5246
astUnitMapper_(const char * in,const char * out,const char * in_lab,char ** out_lab,int * status)5247 AstMapping *astUnitMapper_( const char *in, const char *out,
5248 const char *in_lab, char **out_lab, int *status ){
5249 /*
5250 *+
5251 * Name:
5252 * astUnitMapper
5253
5254 * Purpose:
5255 * Create a Mapping between two system of units.
5256
5257 * Type:
5258 * Protected function.
5259
5260 * Synopsis:
5261 * #include "unit.h"
5262 * AstMapping *astUnitMapper( const char *in, const char *out,
5263 * const char *in_lab, char **out_lab )
5264
5265 * Class Membership:
5266 * Unit member function.
5267
5268 * Description:
5269 * This function creates a Mapping between two specified system of
5270 * units. It also modifes a supplied label (which is typically
5271 * the axis label associated with the input units) so that it includes
5272 * any functional change implied by the supplied "in" and "out" units.
5273
5274 * Parameters:
5275 * in
5276 * A string representation of the input units, for instance "km/h".
5277 * See "Unit Representations:" below.
5278 * out
5279 * A string representation of the output units, for instance "m/s".
5280 * See "Unit Representations:" below.
5281 * in_lab
5282 * A label describing the quantity associated with the input units.
5283 * If the "in" string is the Units attribute of an Axis, then
5284 * "in_lab" should be the Label of the same Axis. May be supplied
5285 * NULL in which case "out_lab" is ignored.
5286 * out_lab
5287 * The address at which to return a pointer to a label describing the
5288 * quantity associated with the output units. For instance, if the
5289 * input and output units are "Hz" and "sqrt(Hz)", and the input
5290 * label is "Frequency", then the returned output label will be
5291 * "sqrt( Frequency )". The returned label is stored in dynamically
5292 * allocated memory which should be freed (using astFree) when no longer
5293 * needed.
5294
5295 * Returned Value:
5296 * A pointer to a Mapping which can be used to transform values in the
5297 * "in" system of units into the "out" system of units. The Mapping
5298 * will have 1 input and 1 output.
5299
5300 * Unit Representations:
5301 * The string supplied for "in" and "out" should represent a system of
5302 * units following the recommendations of the FITS WCS paper I
5303 * "Representation of World Coordinates in FITS" (Greisen & Calabretta).
5304 * Various commonly used variants are also allowed.
5305 *
5306 * To summarise, a string describing a system of units should be an
5307 * algebraic expression which combines one or more named units. The
5308 * following functions and operators may be used within these algebraic
5309 * expressions:
5310 *
5311 * - "*": multiplication. A period "." or space " " may also be used
5312 * to represent multiplication (a period is only interpreted as a
5313 * multiplication operator if it is not positioned between two digits,
5314 * and a space is only interpreted as a multiplication operator if it
5315 * occurs between two operands).
5316 * - "/": division.
5317 * - "**": exponentiation. The exponent (i.e. the operand following the
5318 * exponentiation operator) must be a constant. The symbol "^" is also
5319 * interpreted as an exponentiation operator. Exponentiation is also
5320 * implied by an integer following a unit name without any separator
5321 * (e.g. "cm2" is "cm^2").
5322 * - log(): Common logarithm.
5323 * - ln(): Natural logarithm.
5324 * - sqrt(): Square root.
5325 * - exp(): Exponential.
5326 *
5327 * Function names are case insensitive. White space may be included
5328 * within an expression (note that white space between two operands
5329 * will be interpreted as a muiltiplication operator as described
5330 * above). Parentheses may be used to indicate the order in which
5331 * expressions are to be evaluated (normal mathematical precedence is
5332 * used otherwise). The following symbols may be used to represent
5333 * constants:
5334 *
5335 * - "pi"
5336 * - "e"
5337 *
5338 * These symbols are also case in-sensitive.
5339 *
5340 * The above operators and functions are used to combine together one
5341 * or more "unit symbols". The following base unit symbols are recognised:
5342 *
5343 * - "m": metre.
5344 * - "g": gram.
5345 * - "s": second.
5346 * - "rad": radian.
5347 * - "sr": steradian.
5348 * - "K": Kelvin.
5349 * - "mol": mole.
5350 * - "cd": candela.
5351 *
5352 * The following symbols for units derived fro the above basic units are
5353 * recognised:
5354 *
5355 * - "sec": second (1 s)
5356 * - "Hz": Hertz (1/s).
5357 * - "N": Newton (kg m/s**2).
5358 * - "J": Joule (N m).
5359 * - "W": Watt (J/s).
5360 * - "C": Coulomb (A s).
5361 * - "V": Volt (J/C).
5362 * - "Pa": Pascal (N/m**2).
5363 * - "Ohm": Ohm (V/A).
5364 * - "S": Siemens (A/V).
5365 * - "F": Farad (C/V).
5366 * - "Wb": Weber (V s).
5367 * - "T": Tesla (Wb/m**2).
5368 * - "H": Henry (Wb/A).
5369 * - "lm": lumen (cd sr).
5370 * - "lx": lux (lm/m**2).
5371 * - "deg": degree (pi/180 rad).
5372 * - "arcmin": arc-minute (1/60 deg).
5373 * - "arcsec": arc-second (1/3600 deg).
5374 * - "mas": milli-arcsecond (1/3600000 deg).
5375 * - "min": minute (60 s).
5376 * - "h": hour (3600 s).
5377 * - "d": day (86400 s).
5378 * - "yr": year (31557600 s).
5379 * - "a": year (31557600 s).
5380 * - "eV": electron-Volt (1.60217733E-19 J).
5381 * - "erg": erg (1.0E-7 J).
5382 * - "Ry": Rydberg (13.605692 eV).
5383 * - "solMass": solar mass (1.9891E30 kg).
5384 * - "u": unified atomic mass unit (1.6605387E-27 kg).
5385 * - "solLum": solar luminosity (3.8268E26 W).
5386 * - "Angstrom": Angstrom (1.0E-10 m).
5387 * - "Ang": Angstrom
5388 * - "A": Ampere
5389 * - "micron": micron (1.0E-6 m).
5390 * - "solRad": solar radius (6.9599E8 m).
5391 * - "AU": astronomical unit (1.49598E11 m).
5392 * - "lyr": light year (9.460730E15 m).
5393 * - "pc": parsec (3.0867E16 m).
5394 * - "count": count.
5395 * - "ct": count.
5396 * - "adu": analogue-to-digital converter unit.
5397 * - "photon": photon.
5398 * - "ph": photon.
5399 * - "Jy": Jansky (1.0E-26 W /m**2 /Hz).
5400 * - "Jan": Jansky
5401 * - "mag": magnitude.
5402 * - "G": Gauss (1.0E-4 T).
5403 * - "pixel": pixel.
5404 * - "pix": pixel.
5405 * - "barn": barn (1.0E-28 m**2).
5406 * - "D": Debye (1.0E-29/3 C.m).
5407 *
5408 * In addition, any other unknown unit symbol may be used (but of course
5409 * no mapping will be possible between unknown units).
5410 *
5411 * Unit symbols may be preceded with a numerical constant (for
5412 * instance "1000 m") or a standard multiplier symbol (for instance "km")
5413 * to represent some multiple of the unit. The following standard
5414 * multipliers are recognised:
5415 *
5416 * - "d": deci (1.0E-1)
5417 * - "c": centi (1.0E-2)
5418 * - "m": milli (1.0E-3)
5419 * - "u": micro (1.0E-6)
5420 * - "n": nano (1.0E-9)
5421 * - "p": pico (1.0E-12)
5422 * - "f": femto (1.0E-15)
5423 * - "a": atto (1.0E-18)
5424 * - "z": zepto (1.0E-21)
5425 * - "y": yocto (1.0E-24)
5426 * - "da": deca (1.0E1)
5427 * - "h": hecto (1.0E2)
5428 * - "k": kilo (1.0E3)
5429 * - "M": mega (1.0E6)
5430 * - "G": giga (1.0E9)
5431 * - "T": tera (1.0E12)
5432 * - "P": peta (1.0E15)
5433 * - "E": exa (1.0E18)
5434 * - "Z": zetta (1.0E21)
5435 * - "Y": yotta (1.0E24)
5436
5437 * Notes:
5438 * - NULL values are returned without error if the supplied units are
5439 * incompatible (for instance, if the input and output units are "kg"
5440 * and "m" ).
5441 * - NULL values are returned if this function is invoked with the
5442 * global error status set or if it should fail for any reason.
5443 *-
5444 */
5445
5446 /* Local Variables: */
5447 AstMapping *result;
5448 UnitNode **units;
5449 UnitNode *in_tree;
5450 UnitNode *intemp;
5451 UnitNode *inv;
5452 UnitNode *labtree;
5453 UnitNode *newtest;
5454 UnitNode *out_tree;
5455 UnitNode *outtemp;
5456 UnitNode *src;
5457 UnitNode *testtree;
5458 UnitNode *tmp;
5459 UnitNode *totaltree;
5460 UnitNode *totlabtree;
5461 const char *c;
5462 const char *exp;
5463 int i;
5464 int nc;
5465 int nunits;
5466 int ipass;
5467
5468 /* Initialise */
5469 result = NULL;
5470 if( in_lab ) *out_lab = NULL;
5471
5472 /* Check the global error status. */
5473 if ( !astOK ) return result;
5474
5475 /* A quick check for a common simple case: if the two strings are
5476 identical, return a UnitMap.*/
5477 if( !strcmp( in, out ) ) {
5478 if( in_lab ) *out_lab = astStore( NULL, in_lab, strlen( in_lab ) + 1 );
5479 return (AstMapping *) astUnitMap( 1, "", status );
5480 }
5481
5482 /* More initialisation. */
5483 in_tree = NULL;
5484 out_tree = NULL;
5485 units = NULL;
5486
5487 /* Parse the input units string, producing a tree of UnitNodes which
5488 represents the input units. A pointer to the UnitNode at the head of
5489 the tree is returned if succesfull. Report a context message if this
5490 fails. The returned tree contains branch nodes which correspond to
5491 operators or functions, and leaf nodes which represent constant values
5492 or named basic units (m, s, g, K, etc). Each branch node has one or more
5493 "arguments" (i.e. child nodes) which are operated on or combined by
5494 the branch node in some way to produce the nodes "value". This value
5495 is then used as an argument for the node's parent node (if any). If
5496 the string supplied by the user refers to any known derived units (e.g. "N",
5497 Newton) then each such unit is represented in the returned tree by a
5498 complete sub-tree in which the head node corresponds to the derived
5499 unit (e.g. "N") and the leaf nodes correspond to the basic units needed
5500 to define the derived unit ( for instance, "m", "s" and "g" - metres,
5501 seconds and grammes), or numerical constants. Thus every leaf node in the
5502 returned tree will be a basic unit (i.e. a unit which is not defined in
5503 terms of other units), or a numerical constant. */
5504 in_tree = CreateTree( in, 1, 1, status );
5505 if( !astOK ) astError( AST__BADUN, "astUnitMapper: Error parsing input "
5506 "units string '%s'.", status, in );
5507
5508 /* Do the same for the output units. */
5509 if( astOK ) {
5510 out_tree = CreateTree( out, 1, 1, status );
5511 if( !astOK ) astError( AST__BADUN, "astUnitMapper: Error parsing output "
5512 "units string '%s'.", status, out );
5513 }
5514
5515 /* If a blank string is supplied for both input and output units, then
5516 assume a UnitMap is the appropriate Mapping. */
5517 if( !in_tree && !out_tree && astOK ) {
5518 result = (AstMapping *) astUnitMap( 1, "", status );
5519 if( in_lab ) *out_lab = astStore( NULL, in_lab, strlen( in_lab ) + 1 );
5520
5521 /* Otherwise, if we have both input and output trees... */
5522 } else if( in_tree && out_tree && astOK ) {
5523
5524 /* Locate all the basic units used within either of these two trees. An
5525 array is formed in which each element is a pointer to a UnitNode
5526 contained within one of the trees created above. Each basic unit
5527 referred to in either tree will have a single entry in this array
5528 (even if the unit is referred to more than once). */
5529 units = NULL;
5530 nunits = 0;
5531 LocateUnits( in_tree, &units, &nunits, status );
5532 LocateUnits( out_tree, &units, &nunits, status );
5533
5534 /* Due to the simple nature of the simplification process in SimplifyTree,
5535 the following alogorithm sometimes fails to find a Mapping form input
5536 to output units, but can find a Mapping from output to input units.
5537 In this latter case, we can get the required Mapping from input to
5538 output simply by inverting the Mapign from output to input. So try
5539 first with the units in the original order. If this fails to find a
5540 Mapping, try again with the units swapped, and note that the final
5541 Mapping should be inverted before being used. */
5542 for( ipass = 0; ipass < 2; ipass++ ){
5543 if( ipass == 1 ) {
5544 tmp = in_tree;
5545 in_tree = out_tree;
5546 out_tree = tmp;
5547 }
5548
5549 /* We are going to create a new tree of UnitNodes in which the head node
5550 corresponds to the requested output units, and which has a single
5551 non-constant leaf node corresponding to the input units. Initialise a
5552 pointer to this new tree to indicate that it has not yet been created. */
5553 testtree = NULL;
5554
5555 /* Loop round each basic unit used in the definition of either the input
5556 or the output units (i.e. the elements of the array created above by
5557 "LocateUnits"). The unit selected by this loop is referred to as the
5558 "current" unit. On each pass through this loop, we create a tree which
5559 is a candidate for the final required tree (the "test tree" pointed to
5560 by the testtree pointer initialised above). In order for a mapping to
5561 be possible between input and output units, the test tree created on
5562 each pass through this loop must be equivalent to the test tree for the
5563 previous pass (in other words, all the test trees must be equivalent).
5564 We break out of the loop (and return a NULL Mapping) as soon as we find
5565 a test tree which differs from the previous test tree. */
5566 for( i = 0; i < nunits; i++ ) {
5567
5568 /* Create copies of the trees describing the input and output units, in which
5569 all units other than the current unit are set to a constant value of 1.
5570 This is done by replacing OP_LDVAR nodes (i.e. nodes which "load" the
5571 value of a named basic unit) by OP_LDCON nodes (i.e. nodes which load
5572 a specified constant value) in the tree copy. */
5573 intemp = FixUnits( in_tree, units[ i ], status );
5574 outtemp = FixUnits( out_tree, units[ i ], status );
5575
5576 /* Simplify these trees. An important side-effect of this simplification
5577 is that trees are "standardised" which allows them to be compared for
5578 equivalence. A single mathematical expression can often be represented
5579 in many different ways (for instance "A/B" is equivalent to "(B**(-1))*A").
5580 Standardisation is a process of forcing all equivalent representations
5581 into a single "standard" form. Without standardisation, trees representing
5582 the above two expressions would not be considered to be equivalent
5583 since thy would contain different nodes and have different structures.
5584 As a consequence of this standardisation, the "simplification" performed
5585 by SimplifyTree can sometimes actually make the tree more complicated
5586 (in terms of the number of nodes in the tree). */
5587 SimplifyTree( &intemp, 1, status );
5588 SimplifyTree( &outtemp, 1, status );
5589
5590 /* If either of the simplified trees does not depend on the current unit,
5591 then the node at the head of the simplified tree will have a constant
5592 value (because all the units other than the current unit have been fixed
5593 to a constant value of 1.0 above by FixUnits, leaving only the current
5594 unit to vary in value). If both simplified trees are constants, then
5595 neither tree depends on the current basic unit (i.e. references to the
5596 current basic unit cancel out within each string expression - for
5597 instance if converting from "m.s.Hz" to "km" and the current unit
5598 is "s", then the "s.Hz" term will cause the "s" units to cancel out). In
5599 this case ignore this basic unit and pass on to the next. */
5600 if( outtemp->con != AST__BAD && intemp->con != AST__BAD ) {
5601
5602 /* If just one simplified tree is constant, then the two units cannot
5603 match since one depends on the current basic unit and the other does
5604 not. Free any test tree from previous passes and break out of the loop. */
5605 } else if( outtemp->con != AST__BAD || intemp->con != AST__BAD ) {
5606 intemp = FreeTree( intemp, status );
5607 outtemp = FreeTree( outtemp, status );
5608 testtree = FreeTree( testtree, status );
5609 break;
5610
5611 /* If neither simplified tree is constant, both depend on the current
5612 basic unit and so we can continue to see if their dependencies are
5613 equivalent. */
5614 } else {
5615
5616 /* We are going to create a new tree which is the inverse of the above
5617 simplified "intemp" tree. That is, the new tree will have a head node
5618 corresponding to the current unit, and a single non-constant leaf node
5619 corresponding to the input units. Create an OP_LDVAR node which can be
5620 used as the leaf node for this inverted tree. If the input tree is
5621 inverted successfully, this root node becomes part of the inverted tree,
5622 and so does not need to be freed explicitly (it will be freed when the
5623 inverted tree is freed). */
5624 src = NewNode( NULL, OP_LDVAR, status );
5625 if( astOK ) src->name = astStore( NULL, "input_units", 12 );
5626
5627 /* Now produce the inverted input tree. If the tree cannot be inverted, a
5628 null pointer is returned. Check for this. Otherwise a pointer to the
5629 UnitNode at the head of the inverted tree is returned. */
5630 inv = InvertTree( intemp, src, status );
5631 if( inv ) {
5632
5633 /* Concatenate this tree (which goes from "input units" to "current unit")
5634 with the simplified output tree (which goes from "current unit" to
5635 "output units"), to get a new tree which goes from input units to output
5636 units. */
5637 totaltree = ConcatTree( inv, outtemp, status );
5638
5639 /* Simplify this tree. */
5640 SimplifyTree( &totaltree, 1, status );
5641
5642 /* Compare this simplified tree with the tree produced for the previous
5643 unit (if any). If they differ, we cannot map between the supplied
5644 units so annul the test tree and break out of the loop. If this is the
5645 first unit to be tested, use the total tree as the test tree for the
5646 next unit. */
5647 if( testtree ) {
5648 if( CmpTree( totaltree, testtree, 0, status ) ) testtree = FreeTree( testtree, status );
5649 totaltree = FreeTree( totaltree, status );
5650 if( !testtree ) break;
5651 } else {
5652 testtree = totaltree;
5653 }
5654 }
5655
5656 /* If the input tree was inverted, free the inverted tree. */
5657 if( inv ) {
5658 inv = FreeTree( inv, status );
5659
5660 /* If the input tree could not be inverted, we cannot convert between input
5661 and output units. Free the node which was created to be the root of the
5662 inverted tree (and which has consequently not been incorporated into the
5663 inverted tree), free any testtree and break out of the loop. */
5664 } else {
5665 src = FreeTree( src, status );
5666 testtree = FreeTree( testtree, status );
5667 break;
5668 }
5669 }
5670
5671 /* Free the other trees. */
5672 intemp = FreeTree( intemp, status );
5673 outtemp = FreeTree( outtemp, status );
5674
5675 }
5676
5677 /* If all the basic units used by either of the supplied system of units
5678 produced the same test tree, leave the "swap in and out units" loop. */
5679 if( testtree ) break;
5680
5681 }
5682
5683 /* If the input and output units have been swapped, swap them back to
5684 their original order, and invert the test tree (if there is one). */
5685 if( ipass > 0 ) {
5686 tmp = in_tree;
5687 in_tree = out_tree;
5688 out_tree = tmp;
5689 if( testtree ) {
5690 src = NewNode( NULL, OP_LDVAR, status );
5691 if( astOK ) src->name = astStore( NULL, "input_units", 12 );
5692 newtest = InvertTree( testtree, src, status );
5693 FreeTree( testtree, status );
5694 testtree = newtest;
5695 if( !newtest ) src = FreeTree( src, status );
5696 }
5697 }
5698
5699 /* If all the basic units used by either of the supplied system of units
5700 produced the same test tree, create a Mapping which is equivalent to the
5701 test tree and return it. */
5702 if( testtree ) {
5703 result = MakeMapping( testtree, status );
5704
5705 /* We now go on to produce the output axis label from the supplied input
5706 axis label. Get a tree of UnitNodes which describes the supplied label
5707 associated with the input axis. The tree will have single OP_LDVAR node
5708 corresponding to the basic label (i.e. the label without any single
5709 argument functions or exponentiation operators applied). */
5710 if( in_lab && astOK ) {
5711
5712 /* Get a pointer to the first non-blank character, and store the number of
5713 characters to examine (this excludes any trailing white space). */
5714 exp = in_lab;
5715 while( isspace( *exp ) ) exp++;
5716 c = exp + strlen( exp ) - 1;
5717 while( c >= exp && isspace( *c ) ) c--;
5718 nc = c - exp + 1;
5719
5720 /* Create the tree. */
5721 labtree = MakeLabelTree( exp, nc, status );
5722 if( astOK ) {
5723
5724 /* Concatenate this tree (which goes from "basic label" to "input label")
5725 with the test tree found above (which goes from "input units" to "output
5726 units"), to get a tree which goes from basic label to output label. */
5727 totlabtree = ConcatTree( labtree, testtree, status );
5728
5729 /* Simplify this tree. */
5730 SimplifyTree( &totlabtree, 1, status );
5731
5732 /* Create the output label from this tree. */
5733 *out_lab = (char *) MakeExp( totlabtree, 0, 1, status );
5734
5735 /* Free the trees. */
5736 totlabtree = FreeTree( totlabtree, status );
5737 labtree = FreeTree( labtree, status );
5738
5739 /* Report a context error if the input label could not be parsed. */
5740 } else {
5741 astError( AST__BADUN, "astUnitMapper: Error parsing axis "
5742 "label '%s'.", status, in_lab );
5743 }
5744 }
5745
5746 /* Free the units tree. */
5747 testtree = FreeTree( testtree, status );
5748
5749 }
5750 }
5751
5752 /* Free resources. */
5753 in_tree = FreeTree( in_tree, status );
5754 out_tree = FreeTree( out_tree, status );
5755 units = astFree( units );
5756
5757 /* If an error has occurred, annul the returned Mapping. */
5758 if( !astOK ) {
5759 result = astAnnul( result );
5760 if( in_lab ) *out_lab = astFree( *out_lab );
5761 }
5762
5763 /* Return the result. */
5764 return result;
5765
5766 }
5767
astUnitNormaliser_(const char * in,int * status)5768 const char *astUnitNormaliser_( const char *in, int *status ){
5769 /*
5770 *+
5771 * Name:
5772 * astUnitNormalizer
5773
5774 * Purpose:
5775 * Normalise a unit string into FITS-WCS format.
5776
5777 * Type:
5778 * Protected function.
5779
5780 * Synopsis:
5781 * #include "unit.h"
5782 * const char *astUnitNormaliser( const char *in )
5783
5784 * Class Membership:
5785 * Unit member function.
5786
5787 * Description:
5788 * This function returns a standard FITS-WCS form of the supplied unit
5789 * string.
5790
5791 * Parameters:
5792 * in
5793 * A string representation of the units, for instance "km/h".
5794
5795 * Returned Value:
5796 * A pointer to a dynamically allocated string holding the normalized
5797 * unit string. It should be freed using astFree when no longer needed.
5798
5799 * Notes:
5800 * - An error will be reported if the units string cannot be parsed.
5801 * - NULL is returned if this function is invoked with the
5802 * global error status set or if it should fail for any reason.
5803 *-
5804 */
5805
5806 /* Local Variables: */
5807 UnitNode *in_tree;
5808 double dval;
5809 const char *result;
5810
5811 /* Initialise */
5812 result = NULL;
5813
5814 /* Check the global error status. */
5815 if ( !astOK ) return result;
5816
5817 /* Parse the input units string, producing a tree of UnitNodes which
5818 represents the input units. A pointer to the UnitNode at the head of
5819 the tree is returned if succesfull. Report a context message if this
5820 fails. */
5821 in_tree = CreateTree( in, 0, 1, status );
5822 if( in_tree ) {
5823
5824 /* Simplify the units expression, only doing genuine simplifications. */
5825 SimplifyTree( &in_tree, 1, status );
5826
5827 /* Invert literal constant unit multipliers. This is because a constant of
5828 say 1000 for a unit of "m" means "multiply the value in metres by 1000",
5829 but a unit string of "1000 m" means "value in units of 1000 m" (i.e.
5830 *divide* the value in metres by 1000). */
5831 InvertConstants( &in_tree, status );
5832
5833 /* Convert the tree into string form. */
5834 result = MakeExp( in_tree, 2, 1, status );
5835
5836 /* If the result is a constant value, return a blank string. */
5837 if( 1 == astSscanf( result, "%lg", &dval ) ) {
5838 *((char *) result) = 0;
5839 }
5840
5841 /* Free the tree. */
5842 in_tree = FreeTree( in_tree, status );
5843
5844 } else {
5845 astError( AST__BADUN, "astUnitNormaliser: Error parsing input "
5846 "units string '%s'.", status, in );
5847 }
5848
5849 /* Return the result */
5850 return result;
5851 }
5852
Ustrcmp(const char * a,const char * b,int * status)5853 static int Ustrcmp( const char *a, const char *b, int *status ){
5854 /*
5855 * Name:
5856 * Ustrcmp
5857
5858 * Purpose:
5859 * A case blind version of strcmp.
5860
5861 * Type:
5862 * Private function.
5863
5864 * Synopsis:
5865 * #include "unit.h"
5866 * int Ustrcmp( const char *a, const char *b, int *status )
5867
5868 * Class Membership:
5869 * Unit member function.
5870
5871 * Description:
5872 * Returns 0 if there are no differences between the two strings, and 1
5873 * otherwise. Comparisons are case blind.
5874
5875 * Parameters:
5876 * a
5877 * Pointer to first string.
5878 * b
5879 * Pointer to second string.
5880 * status
5881 * Pointer to the inherited status variable.
5882
5883 * Returned Value:
5884 * Zero if the strings match, otherwise one.
5885
5886 * Notes:
5887 * - This function does not consider the sign of the difference between
5888 * the two strings, whereas "strcmp" does.
5889 * - This function attempts to execute even if an error has occurred.
5890
5891 */
5892
5893 /* Local Variables: */
5894 const char *aa; /* Pointer to next "a" character */
5895 const char *bb; /* Pointer to next "b" character */
5896 int ret; /* Returned value */
5897
5898 /* Initialise the returned value to indicate that the strings match. */
5899 ret = 0;
5900
5901 /* Initialise pointers to the start of each string. */
5902 aa = a;
5903 bb = b;
5904
5905 /* Loop round each character. */
5906 while( 1 ){
5907
5908 /* We leave the loop if either of the strings has been exhausted. */
5909 if( !(*aa ) || !(*bb) ){
5910
5911 /* If one of the strings has not been exhausted, indicate that the
5912 strings are different. */
5913 if( *aa || *bb ) ret = 1;
5914
5915 /* Break out of the loop. */
5916 break;
5917
5918 /* If neither string has been exhausted, convert the next characters to
5919 upper case and compare them, incrementing the pointers to the next
5920 characters at the same time. If they are different, break out of the
5921 loop. */
5922 } else {
5923
5924 if( toupper( (int) *(aa++) ) != toupper( (int) *(bb++) ) ){
5925 ret = 1;
5926 break;
5927 }
5928
5929 }
5930
5931 }
5932
5933 /* Return the result. */
5934 return ret;
5935
5936 }
5937
Ustrncmp(const char * a,const char * b,size_t n,int * status)5938 static int Ustrncmp( const char *a, const char *b, size_t n, int *status ){
5939 /*
5940 * Name:
5941 * Ustrncmp
5942
5943 * Purpose:
5944 * A case blind version of strncmp.
5945
5946 * Type:
5947 * Private function.
5948
5949 * Synopsis:
5950 * #include "unit.h"
5951 * int Ustrncmp( const char *a, const char *b, size_t n, int *status )
5952
5953 * Class Membership:
5954 * Unit member function.
5955
5956 * Description:
5957 * Returns 0 if there are no differences between the first "n"
5958 * characters of the two strings, and 1 otherwise. Comparisons are
5959 * case blind.
5960
5961 * Parameters:
5962 * a
5963 * Pointer to first string.
5964 * b
5965 * Pointer to second string.
5966 * n
5967 * The maximum number of characters to compare.
5968 * status
5969 * Pointer to the inherited status variable.
5970
5971 * Returned Value:
5972 * Zero if the strings match, otherwise one.
5973
5974 * Notes:
5975 * - This function does not consider the sign of the difference between
5976 * the two strings, whereas "strncmp" does.
5977 * - This function attempts to execute even if an error has occurred.
5978
5979 */
5980
5981 /* Local Variables: */
5982 const char *aa; /* Pointer to next "a" character */
5983 const char *bb; /* Pointer to next "b" character */
5984 int i; /* Character index */
5985 int ret; /* Returned value */
5986
5987
5988 /* Initialise the returned value to indicate that the strings match. */
5989 ret = 0;
5990
5991 /* Check pointer have been supplied. */
5992 if( !a || !b ) return ret;
5993
5994 /* Initialise pointers to the start of each string. */
5995 aa = a;
5996 bb = b;
5997
5998 /* Compare up to "n" characters. */
5999 for( i = 0; i < (int) n; i++ ){
6000
6001 /* We leave the loop if either of the strings has been exhausted. */
6002 if( !(*aa ) || !(*bb) ){
6003
6004 /* If one of the strings has not been exhausted, indicate that the
6005 strings are different. */
6006 if( *aa || *bb ) ret = 1;
6007
6008 /* Break out of the loop. */
6009 break;
6010
6011 /* If neither string has been exhausted, convert the next characters to
6012 upper case and compare them, incrementing the pointers to the next
6013 characters at the same time. If they are different, break out of the
6014 loop. */
6015 } else {
6016
6017 if( toupper( (int) *(aa++) ) != toupper( (int) *(bb++) ) ){
6018 ret = 1;
6019 break;
6020 }
6021
6022 }
6023
6024 }
6025
6026 /* Return the result. */
6027 return ret;
6028
6029 }
6030
6031
6032
6033
6034
6035
6036
6037
6038
6039
6040
6041
6042 /* The rest of this file contains functions which are of use for debugging
6043 this module. They are usually commented out.
6044
6045 static const char *DisplayTree( UnitNode *node, int ind ) {
6046 int i;
6047 char buf[200];
6048 const char *result;
6049 char *a;
6050 const char *arg[ 2 ];
6051 int rl;
6052 int slen;
6053 const opsym[ 100 ];
6054
6055 result = "";
6056
6057 for( i = 0; i < ind; i++ ) buf[ i ] = ' ';
6058 buf[ ind ] = 0;
6059
6060 if( !node ) {
6061 printf( "%s <null>\n", buf );
6062 } else {
6063
6064 printf( "%s Code: '%s' (%d)\n", buf, OpName( node->opcode ), node->opcode );
6065 printf( "%s Narg: %d\n", buf, node->narg );
6066 printf( "%s Constant: %g\n", buf, node->con );
6067 printf( "%s Name: %s\n", buf, node->name?node->name:"" );
6068 printf( "%s Unit: %s\n", buf, node->unit?node->unit->sym:"" );
6069 printf( "%s Mult: %s\n", buf, node->mult?node->mult->sym:"" );
6070
6071 OpSym( node, opsym );
6072 slen = strlen( opsym );
6073 rl = slen;
6074
6075 if( node->narg == 0 ) {
6076 result = astMalloc( rl + 1 );
6077 if( astOK ) strcpy( (char *) result, opsym );
6078
6079 } else if( node->narg == 1 ) {
6080 rl += 2;
6081 printf( "%s Arg 0:\n", buf );
6082 arg[ 0 ] = DisplayTree( (node->arg)[ 0 ], ind + 2 );
6083 rl += strlen( arg[ 0 ] );
6084
6085 result = astMalloc( rl + 1 );
6086 if( astOK ) {
6087 a = (char *) result;
6088 strcpy( a, opsym );
6089 a += slen;
6090 *(a++) = '(';
6091 strcpy( a, arg[0] );
6092 a += strlen( arg[ 0 ] );
6093 *(a++) = ')';
6094 }
6095
6096 } else {
6097 rl += 4;
6098 for( i = 0; i < node->narg; i++ ) {
6099 printf( "%s Arg %d:\n", buf, i );
6100 arg[ i ] = DisplayTree( (node->arg)[ i ], ind + 2 );
6101 rl += strlen( arg[ i ] );
6102 }
6103
6104 result = astMalloc( rl + 1 );
6105 if( astOK ) {
6106 a = (char *) result;
6107 *(a++) = '(';
6108 strcpy( a, arg[0] );
6109 a += strlen( arg[ 0 ] );
6110 *(a++) = ')';
6111 strcpy( a, opsym );
6112 a += slen;
6113 *(a++) = '(';
6114 strcpy( a, arg[1] );
6115 a += strlen( arg[ 1 ] );
6116 *(a++) = ')';
6117 }
6118 }
6119 }
6120
6121 if( !astOK ) {
6122 astFree( (void *) result );
6123 result = "";
6124 }
6125
6126 return result;
6127 }
6128
6129 static const char *OpName( Oper op ) {
6130 const char *name;
6131
6132 if( op == OP_LDCON ) {
6133 name = "LDCON";
6134 } else if( op == OP_LDVAR ) {
6135 name = "LDVAR";
6136 } else if( op == OP_LOG ) {
6137 name = "LOG";
6138 } else if( op == OP_LN ) {
6139 name = "LN";
6140 } else if( op == OP_EXP ) {
6141 name = "EXP";
6142 } else if( op == OP_SQRT ) {
6143 name = "SQRT";
6144 } else if( op == OP_POW ) {
6145 name = "POW";
6146 } else if( op == OP_DIV ) {
6147 name = "DIV";
6148 } else if( op == OP_MULT ) {
6149 name = "MULT";
6150 } else if( op == OP_LDPI ) {
6151 name = "LDPI";
6152 } else if( op == OP_LDE ) {
6153 name = "LDE";
6154 } else if( op == OP_NULL ) {
6155 name = "NULL";
6156 } else {
6157 name = "<unknown op code>";
6158 }
6159
6160 return name;
6161 }
6162
6163 static void OpSym( UnitNode *node, char *buff ) {
6164 const char *sym = NULL;
6165
6166 if( node->con != AST__BAD ) {
6167 sprintf( buff, "%g", node->con );
6168
6169 } else if( node->opcode == OP_LDVAR ) {
6170 sym = node->name;
6171
6172 } else if( node->opcode == OP_LOG ) {
6173 sym = "log";
6174
6175 } else if( node->opcode == OP_LN ) {
6176 sym = "ln";
6177
6178 } else if( node->opcode == OP_EXP ) {
6179 sym = "exp";
6180
6181 } else if( node->opcode == OP_SQRT ) {
6182 sym = "sqrt";
6183
6184 } else if( node->opcode == OP_POW ) {
6185 sym = "**";
6186
6187 } else if( node->opcode == OP_DIV ) {
6188 sym = "/";
6189
6190 } else if( node->opcode == OP_MULT ) {
6191 sym = "*";
6192
6193 } else if( node->opcode == OP_NULL ) {
6194 sym = "NULL";
6195
6196 } else {
6197 sym = "<unknown op code>";
6198 }
6199
6200 if( sym ) strcpy( buff, sym );
6201 }
6202
6203 static const char *TreeExp( UnitNode *node ) {
6204 char buff[ 100 ];
6205 char buff2[ 100 ];
6206
6207 if( node->narg == 0 ) {
6208 OpSym( node, buff );
6209
6210 } else if( node->narg == 1 ) {
6211 OpSym( node, buff2 );
6212 sprintf( buff, "%s(%s)", buff2, TreeExp( node->arg[ 0 ] ) );
6213
6214 } else if( node->narg == 2 ) {
6215 OpSym( node, buff2 );
6216 sprintf( buff, "(%s)%s(%s)", TreeExp( node->arg[ 0 ] ), buff2,
6217 TreeExp( node->arg[ 1 ] ) );
6218 }
6219
6220 return astStore( NULL, buff, strlen( buff ) + 1 );
6221 }
6222
6223 */
6224
6225
6226
6227
6228