diff options
author | William Joye <wjoye@cfa.harvard.edu> | 2018-09-21 17:04:03 (GMT) |
---|---|---|
committer | William Joye <wjoye@cfa.harvard.edu> | 2018-09-21 17:04:03 (GMT) |
commit | 28518ab5eb4726fe91ee0c916b2322f8c47eb3ad (patch) | |
tree | 2653f0ce7e9c70973f27fe6c0c38ec31e21ce5cf /ast/mathmap.c | |
parent | bd7d67f66c53df36bb50e3423bfc91eae8618201 (diff) | |
download | blt-28518ab5eb4726fe91ee0c916b2322f8c47eb3ad.zip blt-28518ab5eb4726fe91ee0c916b2322f8c47eb3ad.tar.gz blt-28518ab5eb4726fe91ee0c916b2322f8c47eb3ad.tar.bz2 |
update ast 8.6.3
Diffstat (limited to 'ast/mathmap.c')
-rw-r--r-- | ast/mathmap.c | 7421 |
1 files changed, 7421 insertions, 0 deletions
diff --git a/ast/mathmap.c b/ast/mathmap.c new file mode 100644 index 0000000..f3a21f9 --- /dev/null +++ b/ast/mathmap.c @@ -0,0 +1,7421 @@ +/* +*class++ +* Name: +* MathMap + +* Purpose: +* Transform coordinates using mathematical expressions. + +* Constructor Function: +c astMathMap +f AST_MATHMAP + +* Description: +c A MathMap is a Mapping which allows you to specify a set of forward +c and/or inverse transformation functions using arithmetic operations +c and mathematical functions similar to those available in C. The +c MathMap interprets these functions at run-time, whenever its forward +c or inverse transformation is required. Because the functions are not +c compiled in the normal sense (unlike an IntraMap), they may be used to +c describe coordinate transformations in a transportable manner. A +c MathMap therefore provides a flexible way of defining new types of +c Mapping whose descriptions may be stored as part of a dataset and +c interpreted by other programs. +f A MathMap is a Mapping which allows you to specify a set of forward +f and/or inverse transformation functions using arithmetic operations +f and mathematical functions similar to those available in Fortran. The +f MathMap interprets these functions at run-time, whenever its forward +f or inverse transformation is required. Because the functions are not +f compiled in the normal sense (unlike an IntraMap), they may be used to +f describe coordinate transformations in a transportable manner. A +f MathMap therefore provides a flexible way of defining new types of +f Mapping whose descriptions may be stored as part of a dataset and +f interpreted by other programs. + +* Inheritance: +* The MathMap class inherits from the Mapping class. + +* Attributes: +* In addition to those attributes common to all Mappings, every +* MathMap also has the following attributes: +* - Seed: Random number seed +* - SimpFI: Forward-inverse MathMap pairs simplify? +* - SimpIF: Inverse-forward MathMap pairs simplify? + +* Functions: +c The MathMap class does not define any new functions beyond those +f The MathMap class does not define any new routines beyond those +* which are applicable to all Mappings. + +* Copyright: +* Copyright (C) 1997-2006 Council for the Central Laboratory of the +* Research Councils + +* Licence: +* This program is free software: you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation, either +* version 3 of the License, or (at your option) any later +* version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU Lesser General Public License for more details. +* +* You should have received a copy of the GNU Lesser General +* License along with this program. If not, see +* <http://www.gnu.org/licenses/>. + +* Authors: +* RFWS: R.F. Warren-Smith (Starlink) + +* History: +* 3-SEP-1999 (RFWS): +* Original version. +* 8-JAN-2003 (DSB): +* Changed private InitVtab method to protected astInitMathMapVtab +* method. +* 14-FEB-2006 (DSB): +* Override astGetObjSize. +* 14-MAR-2006 (DSB): +* - Add QIF function. +* - Override astEqual method. +* 20-NOV-2006 (DSB): +* Re-implement the Equal method to avoid use of astSimplify. +* 30-AUG-2012 (DSB): +* Fix bug in undocumented Gaussian noise function. +*class-- +*/ + +/* Module Macros. */ +/* ============== */ +/* Set the name of the class we are implementing. This indicates to + the header files that define class interfaces that they should make + "protected" symbols available. */ +#define astCLASS MathMap + +/* Allocate pointer array. */ +/* ----------------------- */ +/* This macro allocates an array of pointers. If successful, each element + of the array is initialised to NULL. */ +#define MALLOC_POINTER_ARRAY(array_name,array_type,array_size) \ +\ +/* Allocate the array. */ \ + (array_name) = astMalloc( sizeof(array_type) * (size_t) (array_size) ); \ + if ( astOK ) { \ +\ +/* If successful, loop to initialise each element. */ \ + int array_index_; \ + for ( array_index_ = 0; array_index_ < (array_size); array_index_++ ) { \ + (array_name)[ array_index_ ] = NULL; \ + } \ + } + +/* Free pointer array. */ +/* ------------------- */ +/* This macro frees a dynamically allocated array of pointers, each of + whose elements may point at a further dynamically allocated array + (which is also to be freed). It also allows for the possibility of any + of the pointers being NULL. */ +#define FREE_POINTER_ARRAY(array_name,array_size) \ +\ +/* Check that the main array pointer is not NULL. */ \ + if ( (array_name) ) { \ +\ +/* If OK, loop to free each of the sub-arrays. */ \ + int array_index_; \ + for ( array_index_ = 0; array_index_ < (array_size); array_index_++ ) { \ +\ +/* Check that each sub-array pointer is not NULL before freeing it. */ \ + if ( (array_name)[ array_index_ ] ) { \ + (array_name)[ array_index_ ] = \ + astFree( (array_name)[ array_index_ ] ); \ + } \ + } \ +\ +/* Free the main pointer array. */ \ + (array_name) = astFree( (array_name) ); \ + } + +/* SizeOf pointer array. */ +/* --------------------- */ +/* This macro increments "result" by the number of bytes allocated for an + array of pointers, each of whose elements may point at a further + dynamically allocated array (which is also to be included). It also + allows for the possibility of any of the pointers being NULL. */ +#define SIZEOF_POINTER_ARRAY(array_name,array_size) \ +\ +/* Check that the main array pointer is not NULL. */ \ + if ( (array_name) ) { \ +\ +/* If OK, loop to measure each of the sub-arrays. */ \ + int array_index_; \ + for ( array_index_ = 0; array_index_ < (array_size); array_index_++ ) { \ +\ +/* Check that each sub-array pointer is not NULL before measuring it. */ \ + if ( (array_name)[ array_index_ ] ) { \ + result += astTSizeOf( (array_name)[ array_index_ ] ); \ + } \ + } \ +\ +/* Include the main pointer array. */ \ + result += astTSizeOf( (array_name) ); \ + } + +/* Header files. */ +/* ============= */ +/* Interface definitions. */ +/* ---------------------- */ +#include "channel.h" /* I/O channels */ + +#include "globals.h" /* Thread-safe global data access */ +#include "error.h" /* Error reporting facilities */ +#include "mapping.h" /* Coordinate mappings (parent class) */ +#include "cmpmap.h" /* Compound Mappings */ +#include "mathmap.h" /* Interface definition for this class */ +#include "memory.h" /* Memory allocation facilities */ +#include "globals.h" /* Thread-safe global data access */ +#include "object.h" /* Base Object class */ +#include "pointset.h" /* Sets of points */ +#include "unitmap.h" /* Unit Mapping */ + +/* Error code definitions. */ +/* ----------------------- */ +#include "ast_err.h" /* AST error codes */ + +/* C header files. */ +/* --------------- */ +#include <ctype.h> +#include <errno.h> +#include <limits.h> +#include <math.h> +#include <stddef.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <time.h> + +/* Module Variables. */ +/* ================= */ +/* This type is made obscure since it is publicly accessible (but not + useful). Provide shorthand for use within this module. */ +typedef AstMathMapRandContext_ Rcontext; + + + +/* Address of this static variable is used as a unique identifier for + member of this class. */ +static int class_check; + +/* Pointers to parent class methods which are extended by this class. */ +static int (* parent_getobjsize)( AstObject *, int * ); +static AstPointSet *(* parent_transform)( AstMapping *, AstPointSet *, int, AstPointSet *, int * ); +static const char *(* parent_getattrib)( AstObject *, const char *, int * ); +static int (* parent_testattrib)( AstObject *, const char *, int * ); +static void (* parent_clearattrib)( AstObject *, const char *, int * ); +static void (* parent_setattrib)( AstObject *, const char *, int * ); + +/* This declaration enumerates the operation codes recognised by the + EvaluateFunction function which evaluates arithmetic expressions. */ +typedef enum { + +/* User-supplied constants and variables. */ + OP_LDCON, /* Load constant */ + OP_LDVAR, /* Load variable */ + +/* System constants. */ + OP_LDBAD, /* Load bad value (AST__BAD) */ + OP_LDDIG, /* Load # decimal digits (AST__DBL_DIG) */ + OP_LDEPS, /* Load relative precision (DBL_EPSILON) */ + OP_LDMAX, /* Load largest value (DBL_MAX) */ + OP_LDMAX10E, /* Max. decimal exponent (DBL_MAX_10_EXP) */ + OP_LDMAXE, /* Load maximum exponent (DBL_MAX_EXP) */ + OP_LDMDIG, /* Load # mantissa digits (DBL_MANT_DIG) */ + OP_LDMIN, /* Load smallest value (DBL_MIN) */ + OP_LDMIN10E, /* Min. decimal exponent (DBL_MIN_10_EXP) */ + OP_LDMINE, /* Load minimum exponent (DBL_MIN_EXP) */ + OP_LDRAD, /* Load floating radix (FLT_RADIX) */ + OP_LDRND, /* Load rounding mode (FLT_ROUNDS) */ + +/* Mathematical constants. */ + OP_LDE, /* Load e (base of natural logarithms) */ + OP_LDPI, /* Load pi */ + +/* Functions with one argument. */ + OP_ABS, /* Absolute value (sign removal) */ + OP_ACOS, /* Inverse cosine (radians) */ + OP_ACOSD, /* Inverse cosine (degrees) */ + OP_ACOSH, /* Inverse hyperbolic cosine */ + OP_ACOTH, /* Inverse hyperbolic cotangent */ + OP_ACSCH, /* Inverse hyperbolic cosecant */ + OP_ASECH, /* Inverse hyperbolic secant */ + OP_ASIN, /* Inverse sine (radians) */ + OP_ASIND, /* Inverse sine (degrees) */ + OP_ASINH, /* Inverse hyperbolic sine */ + OP_ATAN, /* Inverse tangent (radians) */ + OP_ATAND, /* Inverse tangent (degrees) */ + OP_ATANH, /* Inverse hyperbolic tangent */ + OP_CEIL, /* C ceil function (round up) */ + OP_COS, /* Cosine (radians) */ + OP_COSD, /* Cosine (degrees) */ + OP_COSH, /* Hyperbolic cosine */ + OP_COTH, /* Hyperbolic cotangent */ + OP_CSCH, /* Hyperbolic cosecant */ + OP_EXP, /* Exponential function */ + OP_FLOOR, /* C floor function (round down) */ + OP_INT, /* Integer value (round towards zero) */ + OP_ISBAD, /* Test for bad value */ + OP_LOG, /* Natural logarithm */ + OP_LOG10, /* Base 10 logarithm */ + OP_NINT, /* Fortran NINT function (round to nearest) */ + OP_POISS, /* Poisson random number */ + OP_SECH, /* Hyperbolic secant */ + OP_SIN, /* Sine (radians) */ + OP_SINC, /* Sinc function [= sin(x)/x] */ + OP_SIND, /* Sine (degrees) */ + OP_SINH, /* Hyperbolic sine */ + OP_SQR, /* Square */ + OP_SQRT, /* Square root */ + OP_TAN, /* Tangent (radians) */ + OP_TAND, /* Tangent (degrees) */ + OP_TANH, /* Hyperbolic tangent */ + +/* Functions with two arguments. */ + OP_ATAN2, /* Inverse tangent (2 arguments, radians) */ + OP_ATAN2D, /* Inverse tangent (2 arguments, degrees) */ + OP_DIM, /* Fortran DIM (positive difference) fn. */ + OP_GAUSS, /* Gaussian random number */ + OP_MOD, /* Modulus function */ + OP_POW, /* Raise to power */ + OP_RAND, /* Uniformly distributed random number */ + OP_SIGN, /* Transfer of sign function */ + +/* Functions with three arguments. */ + OP_QIF, /* C "question mark" operator "a?b:c" */ + +/* Functions with variable numbers of arguments. */ + OP_MAX, /* Maximum of 2 or more values */ + OP_MIN, /* Minimum of 2 or more values */ + +/* Unary arithmetic operators. */ + OP_NEG, /* Negate (change sign) */ + +/* Unary boolean operators. */ + OP_NOT, /* Boolean NOT */ + +/* Binary arithmetic operators. */ + OP_ADD, /* Add */ + OP_DIV, /* Divide */ + OP_MUL, /* Multiply */ + OP_SUB, /* Subtract */ + +/* Bit-shift operators. */ + OP_SHFTL, /* Shift bits left */ + OP_SHFTR, /* Shift bits right */ + +/* Relational operators. */ + OP_EQ, /* Relational equal */ + OP_GE, /* Greater than or equal */ + OP_GT, /* Greater than */ + OP_LE, /* Less than or equal */ + OP_LT, /* Less than */ + OP_NE, /* Not equal */ + +/* Bit-wise operators. */ + OP_BITAND, /* Bit-wise AND */ + OP_BITOR, /* Bit-wise OR */ + OP_BITXOR, /* Bit-wise exclusive OR */ + +/* Binary boolean operators. */ + OP_AND, /* Boolean AND */ + OP_EQV, /* Fortran logical .EQV. operation */ + OP_OR, /* Boolean OR */ + OP_XOR, /* Boolean exclusive OR */ + +/* Null operation. */ + OP_NULL /* Null operation */ +} Oper; + +/* This structure holds a description of each symbol which may appear + in an expression. */ +typedef struct { + const char *text; /* Symbol text as it appears in expressions */ + const int size; /* Size of symbol text */ + const int operleft; /* An operator when seen from the left? */ + const int operright; /* An operator when seen from the right? */ + const int unarynext; /* May be followed by a unary +/- ? */ + const int unaryoper; /* Is a unary +/- ? */ + const int leftpriority; /* Priority when seen from the left */ + const int rightpriority; /* Priority when seen from the right */ + const int parincrement; /* Change in parenthesis level */ + const int stackincrement; /* Change in evaluation stack size */ + const int nargs; /* Number of function arguments */ + const Oper opcode; /* Resulting operation code */ +} Symbol; + +/* This initialises an array of Symbol structures to hold data on all + the supported symbols. The order is not important, but symbols are + arranged here in approximate order of descending evaluation + priority. The end of the array is indicated by an element with a NULL + "text" component. */ +static const Symbol symbol[] = { + +/* User-supplied constants and variables. */ + { "" , 0, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDCON }, + { "" , 0, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDVAR }, + +/* System constants. */ + { "<bad>" , 5, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDBAD }, + { "<dig>" , 5, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDDIG }, + { "<epsilon>" , 9, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDEPS }, + { "<mant_dig>" , 10, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDMDIG }, + { "<max>" , 5, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDMAX }, + { "<max_10_exp>", 12, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDMAX10E }, + { "<max_exp>" , 9, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDMAXE }, + { "<min>" , 5, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDMIN }, + { "<min_10_exp>", 12, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDMIN10E }, + { "<min_exp>" , 9, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDMINE }, + { "<radix>" , 7, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDRAD }, + { "<rounds>" , 8, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDRND }, + +/* Mathematical constants. */ + { "<e>" , 3, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDE }, + { "<pi>" , 4, 0, 0, 0, 0, 19, 19, 0, 1, 0, OP_LDPI }, + +/* Functions with one argument. */ + { "abs(" , 4, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ABS }, + { "acos(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ACOS }, + { "acosd(" , 6, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ACOSD }, + { "acosh(" , 6, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ACOSH }, + { "acoth(" , 6, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ACOTH }, + { "acsch(" , 6, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ACSCH }, + { "aint(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_INT }, + { "asech(" , 6, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ASECH }, + { "asin(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ASIN }, + { "asind(" , 6, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ASIND }, + { "asinh(" , 6, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ASINH }, + { "atan(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ATAN }, + { "atand(" , 6, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ATAND }, + { "atanh(" , 6, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ATANH }, + { "ceil(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_CEIL }, + { "cos(" , 4, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_COS }, + { "cosd(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_COSD }, + { "cosh(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_COSH }, + { "coth(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_COTH }, + { "csch(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_CSCH }, + { "exp(" , 4, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_EXP }, + { "fabs(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ABS }, + { "floor(" , 6, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_FLOOR }, + { "int(" , 4, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_INT }, + { "isbad(" , 6, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_ISBAD }, + { "log(" , 4, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_LOG }, + { "log10(" , 6, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_LOG10 }, + { "nint(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_NINT }, + { "poisson(" , 8, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_POISS }, + { "sech(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_SECH }, + { "sin(" , 4, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_SIN }, + { "sinc(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_SINC }, + { "sind(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_SIND }, + { "sinh(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_SINH }, + { "sqr(" , 4, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_SQR }, + { "sqrt(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_SQRT }, + { "tan(" , 4, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_TAN }, + { "tand(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_TAND }, + { "tanh(" , 5, 0, 1, 1, 0, 19, 1, 1, 0, 1, OP_TANH }, + +/* Functions with two arguments. */ + { "atan2(" , 6, 0, 1, 1, 0, 19, 1, 1, -1, 2, OP_ATAN2 }, + { "atan2d(" , 7, 0, 1, 1, 0, 19, 1, 1, -1, 2, OP_ATAN2D }, + { "dim(" , 4, 0, 1, 1, 0, 19, 1, 1, -1, 2, OP_DIM }, + { "fmod(" , 5, 0, 1, 1, 0, 19, 1, 1, -1, 2, OP_MOD }, + { "gauss(" , 6, 0, 1, 1, 0, 19, 1, 1, -1, 2, OP_GAUSS }, + { "mod(" , 4, 0, 1, 1, 0, 19, 1, 1, -1, 2, OP_MOD }, + { "pow(" , 4, 0, 1, 1, 0, 19, 1, 1, -1, 2, OP_POW }, + { "rand(" , 5, 0, 1, 1, 0, 19, 1, 1, -1, 2, OP_RAND }, + { "sign(" , 5, 0, 1, 1, 0, 19, 1, 1, -1, 2, OP_SIGN }, + +/* Functions with two arguments. */ + { "qif(" , 4, 0, 1, 1, 0, 19, 1, 1, -2, 3, OP_QIF }, + +/* Functions with variable numbers of arguments. */ + { "max(" , 4, 0, 1, 1, 0, 19, 1, 1, -1, -2, OP_MAX }, + { "min(" , 4, 0, 1, 1, 0, 19, 1, 1, -1, -2, OP_MIN }, + +/* Parenthesised expressions. */ + { ")" , 1, 1, 0, 0, 0, 2, 19, -1, 0, 0, OP_NULL }, + { "(" , 1, 0, 1, 1, 0, 19, 1, 1, 0, 0, OP_NULL }, + +/* Unary arithmetic operators. */ + { "+" , 1, 0, 1, 1, 1, 17, 16, 0, 0, 0, OP_NULL }, + { "-" , 1, 0, 1, 1, 1, 17, 16, 0, 0, 0, OP_NEG }, + +/* Unary boolean operators. */ + { "!" , 1, 0, 1, 1, 0, 17, 16, 0, 0, 0, OP_NOT }, + { ".not." , 5, 0, 1, 1, 0, 17, 16, 0, 0, 0, OP_NOT }, + +/* Binary arithmetic operators. */ + { "**" , 2, 1, 1, 1, 0, 18, 15, 0, -1, 0, OP_POW }, + { "*" , 1, 1, 1, 1, 0, 14, 14, 0, -1, 0, OP_MUL }, + { "/" , 1, 1, 1, 1, 0, 14, 14, 0, -1, 0, OP_DIV }, + { "+" , 1, 1, 1, 1, 0, 13, 13, 0, -1, 0, OP_ADD }, + { "-" , 1, 1, 1, 1, 0, 13, 13, 0, -1, 0, OP_SUB }, + +/* Bit-shift operators. */ + { "<<" , 2, 1, 1, 1, 0, 12, 12, 0, -1, 0, OP_SHFTL }, + { ">>" , 2, 1, 1, 1, 0, 12, 12, 0, -1, 0, OP_SHFTR }, + +/* Relational operators. */ + { "<" , 1, 1, 1, 1, 0, 11, 11, 0, -1, 0, OP_LT }, + { ".lt." , 4, 1, 1, 1, 0, 11, 11, 0, -1, 0, OP_LT }, + { "<=" , 2, 1, 1, 1, 0, 11, 11, 0, -1, 0, OP_LE }, + { ".le." , 4, 1, 1, 1, 0, 11, 11, 0, -1, 0, OP_LE }, + { ">" , 1, 1, 1, 1, 0, 11, 11, 0, -1, 0, OP_GT }, + { ".gt." , 4, 1, 1, 1, 0, 11, 11, 0, -1, 0, OP_GT }, + { ">=" , 2, 1, 1, 1, 0, 11, 11, 0, -1, 0, OP_GE }, + { ".ge." , 4, 1, 1, 1, 0, 11, 11, 0, -1, 0, OP_GE }, + { "==" , 2, 1, 1, 1, 0, 10, 10, 0, -1, 0, OP_EQ }, + { ".eq." , 4, 1, 1, 1, 0, 10, 10, 0, -1, 0, OP_EQ }, + { "!=" , 2, 1, 1, 1, 0, 10, 10, 0, -1, 0, OP_NE }, + { ".ne." , 4, 1, 1, 1, 0, 10, 10, 0, -1, 0, OP_NE }, + +/* Bit-wise operators. */ + { "&" , 1, 1, 1, 1, 0, 9, 9, 0, -1, 0, OP_BITAND }, + { "^" , 1, 1, 1, 1, 0, 8, 8, 0, -1, 0, OP_BITXOR }, + { "|" , 1, 1, 1, 1, 0, 7, 7, 0, -1, 0, OP_BITOR }, + +/* Binary boolean operators. */ + { "&&" , 2, 1, 1, 1, 0, 6, 6, 0, -1, 0, OP_AND }, + { ".and." , 5, 1, 1, 1, 0, 6, 6, 0, -1, 0, OP_AND }, + { "^^" , 2, 1, 1, 1, 0, 5, 5, 0, -1, 0, OP_XOR }, + { "||" , 2, 1, 1, 1, 0, 4, 4, 0, -1, 0, OP_OR }, + { ".or." , 4, 1, 1, 1, 0, 4, 4, 0, -1, 0, OP_OR }, + { ".eqv." , 5, 1, 1, 1, 0, 3, 3, 0, -1, 0, OP_EQV }, + { ".neqv." , 6, 1, 1, 1, 0, 3, 3, 0, -1, 0, OP_XOR }, + { ".xor." , 5, 1, 1, 1, 0, 3, 3, 0, -1, 0, OP_XOR }, + +/* Separators. */ + { "," , 1, 1, 1, 1, 0, 2, 2, 0, 0, 0, OP_NULL }, + +/* End of symbol data. */ + { NULL , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, OP_NULL } +}; + +/* These variables identify indices in the above array which hold + special symbols used explicitly in the code. */ +static const int symbol_ldcon = 0; /* Load a constant */ +static const int symbol_ldvar = 1; /* Load a variable */ + +/* Define macros for accessing each item of thread specific global data. */ +#ifdef THREAD_SAFE + +/* Define how to initialise thread-specific globals. */ +#define GLOBAL_inits \ + globals->Class_Init = 0; \ + globals->GetAttrib_Buff[ 0 ] = 0; + +/* Create the function that initialises global data for this module. */ +astMAKE_INITGLOBALS(MathMap) + +/* Define macros for accessing each item of thread specific global data. */ +#define class_init astGLOBAL(MathMap,Class_Init) +#define class_vtab astGLOBAL(MathMap,Class_Vtab) +#define getattrib_buff astGLOBAL(MathMap,GetAttrib_Buff) + + + +static pthread_mutex_t mutex2 = PTHREAD_MUTEX_INITIALIZER; +#define LOCK_MUTEX2 pthread_mutex_lock( &mutex2 ); +#define UNLOCK_MUTEX2 pthread_mutex_unlock( &mutex2 ); + +static pthread_mutex_t mutex3 = PTHREAD_MUTEX_INITIALIZER; +#define LOCK_MUTEX3 pthread_mutex_lock( &mutex3 ); +#define UNLOCK_MUTEX3 pthread_mutex_unlock( &mutex3 ); + +static pthread_mutex_t mutex4 = PTHREAD_MUTEX_INITIALIZER; +#define LOCK_MUTEX4 pthread_mutex_lock( &mutex4 ); +#define UNLOCK_MUTEX4 pthread_mutex_unlock( &mutex4 ); + +static pthread_mutex_t mutex5 = PTHREAD_MUTEX_INITIALIZER; +#define LOCK_MUTEX5 pthread_mutex_lock( &mutex5 ); +#define UNLOCK_MUTEX5 pthread_mutex_unlock( &mutex5 ); + +static pthread_mutex_t mutex6 = PTHREAD_MUTEX_INITIALIZER; +#define LOCK_MUTEX6 pthread_mutex_lock( &mutex6 ); +#define UNLOCK_MUTEX6 pthread_mutex_unlock( &mutex6 ); + +static pthread_mutex_t mutex7 = PTHREAD_MUTEX_INITIALIZER; +#define LOCK_MUTEX7 pthread_mutex_lock( &mutex7 ); +#define UNLOCK_MUTEX7 pthread_mutex_unlock( &mutex7 ); + +/* If thread safety is not needed, declare and initialise globals at static + variables. */ +#else + +static char getattrib_buff[ 51 ]; + + +/* Define the class virtual function table and its initialisation flag + as static variables. */ +static AstMathMapVtab class_vtab; /* Virtual function table */ +static int class_init = 0; /* Virtual function table initialised? */ + +#define LOCK_MUTEX2 +#define UNLOCK_MUTEX2 + +#define LOCK_MUTEX3 +#define UNLOCK_MUTEX3 + +#define LOCK_MUTEX4 +#define UNLOCK_MUTEX4 + +#define LOCK_MUTEX5 +#define UNLOCK_MUTEX5 + +#define LOCK_MUTEX6 +#define UNLOCK_MUTEX6 + +#define LOCK_MUTEX7 +#define UNLOCK_MUTEX7 + +#endif + + +/* External Interface Function Prototypes. */ +/* ======================================= */ +/* The following functions have public prototypes only (i.e. no + protected prototypes), so we must provide local prototypes for use + within this module. */ +AstMathMap *astMathMapId_( int, int, int, const char *[], int, const char *[], const char *, ... ); + +/* Prototypes for Private Member Functions. */ +/* ======================================== */ +static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * ); +static int GetObjSize( AstObject *, int * ); +static const char *GetAttrib( AstObject *, const char *, int * ); +static double Gauss( Rcontext *, int * ); +static double LogGamma( double, int * ); +static double Poisson( Rcontext *, double, int * ); +static double Rand( Rcontext *, int * ); +static int DefaultSeed( const Rcontext *, int * ); +static int Equal( AstObject *, AstObject *, int * ); +static int GetSeed( AstMathMap *, int * ); +static int GetSimpFI( AstMathMap *, int * ); +static int GetSimpIF( AstMathMap *, int * ); +static int MapMerge( AstMapping *, int, int, int *, AstMapping ***, int **, int * ); +static int TestAttrib( AstObject *, const char *, int * ); +static int TestSeed( AstMathMap *, int * ); +static int TestSimpFI( AstMathMap *, int * ); +static int TestSimpIF( AstMathMap *, int * ); +static void CleanFunctions( int, const char *[], char ***, int * ); +static void ClearAttrib( AstObject *, const char *, int * ); +static void ClearSeed( AstMathMap *, int * ); +static void ClearSimpFI( AstMathMap *, int * ); +static void ClearSimpIF( AstMathMap *, int * ); +static void CompileExpression( const char *, const char *, const char *, int, const char *[], int **, double **, int *, int * ); +static void CompileMapping( const char *, const char *, int, int, int, const char *[], int, const char *[], int ***, int ***, double ***, double ***, int *, int *, int * ); +static void Copy( const AstObject *, AstObject *, int * ); +static void Delete( AstObject *, int * ); +static void Dump( AstObject *, AstChannel *, int * ); +static void EvaluateFunction( Rcontext *, int, const double **, const int *, const double *, int, double *, int * ); +static void EvaluationSort( const double [], int, int [], int **, int *, int * ); +static void ExtractExpressions( const char *, const char *, int, const char *[], int, char ***, int * ); +static void ExtractVariables( const char *, const char *, int, const char *[], int, int, int, int, int, char ***, int * ); +static void ParseConstant( const char *, const char *, const char *, int, int *, double *, int * ); +static void ParseName( const char *, int, int *, int * ); +static void ParseVariable( const char *, const char *, const char *, int, int, const char *[], int *, int *, int * ); +static void SetAttrib( AstObject *, const char *, int * ); +static void SetSeed( AstMathMap *, int, int * ); +static void SetSimpFI( AstMathMap *, int, int * ); +static void SetSimpIF( AstMathMap *, int, int * ); +static void ValidateSymbol( const char *, const char *, const char *, int, int, int *, int **, int **, int *, double **, int * ); + +/* Member functions. */ +/* ================= */ +static void CleanFunctions( int nfun, const char *fun[], char ***clean, int *status ) { +/* +* Name: +* CleanFunctions + +* Purpose: +* Make a clean copy of a set of functions. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* void CleanFunctions( int nfun, const char *fun[], char ***clean, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* This function copies an array of strings, eliminating any white space +* characters and converting to lower case. It is intended for cleaning +* up arrays of function definitions prior to compilation. The returned +* copy is stored in dynamically allocated memory. + +* Parameters: +* nfun +* The number of functions to be cleaned. +* fun +* Pointer to an array, with "nfun" elements, of pointers to null +* terminated strings which contain each of the functions. +* clean +* Address in which to return a pointer to an array (with "nfun" +* elements) of pointers to null terminated strings containing the +* cleaned functions (i.e. this returns an array of strings). +* +* Both the returned array of pointers, and the strings to which they +* point, will be dynamically allocated and should be freed by the +* caller (using astFree) when no longer required. +* status +* Pointer to the inherited status variable. + +* Notes: +* - A NULL value will be returned for "*clean" if this function is +* invoked with the global error status set, or if it should fail for +* any reason. +*/ + +/* Local Variables: */ + char c; /* Character from function string */ + int i; /* Loop counter for characters */ + int ifun; /* Loop counter for functions */ + int nc; /* Count of non-blank characters */ + +/* Initialise. */ + *clean = NULL; + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Allocate and initialise an array to hold the returned pointers. */ + MALLOC_POINTER_ARRAY( *clean, char *, nfun ) + +/* Loop through all the input functions. */ + if ( astOK ) { + for ( ifun = 0; ifun < nfun; ifun++ ) { + +/* Count the number of non-blank characters in each function string. */ + nc = 0; + for ( i = 0; ( c = fun[ ifun ][ i ] ); i++ ) nc += !isspace( c ); + +/* Allocate a string long enough to hold the function with all the + white space removed, storing its pointer in the array allocated + earlier. Check for errors. */ + ( *clean )[ ifun ] = astMalloc( sizeof( char ) * + (size_t) ( nc + 1 ) ); + if ( !astOK ) break; + +/* Loop to copy the non-blank function characters into the new + string. */ + nc = 0; + for ( i = 0; ( c = fun[ ifun ][ i ] ); i++ ) { + if ( !isspace( c ) ) ( *clean )[ ifun ][ nc++ ] = tolower( c ); + } + +/* Null-terminate the result. */ + ( *clean )[ ifun ][ nc ] = '\0'; + } + +/* If an error occurred, then free the main pointer array together + with any strings that have been allocated, resetting the output + value. */ + if ( !astOK ) { + FREE_POINTER_ARRAY( *clean, nfun ) + } + } +} + +static void ClearAttrib( AstObject *this_object, const char *attrib, int *status ) { +/* +* Name: +* ClearAttrib + +* Purpose: +* Clear an attribute value for a MathMap. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* void ClearAttrib( AstObject *this, const char *attrib, int *status ) + +* Class Membership: +* MathMap member function (over-rides the astClearAttrib protected +* method inherited from the Mapping class). + +* Description: +* This function clears the value of a specified attribute for a +* MathMap, so that the default value will subsequently be used. + +* Parameters: +* this +* Pointer to the MathMap. +* attrib +* Pointer to a null terminated string specifying the attribute +* name. This should be in lower case with no surrounding white +* space. +* status +* Pointer to the inherited status variable. +*/ + +/* Local Variables: */ + AstMathMap *this; /* Pointer to the MathMap structure */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Obtain a pointer to the MathMap structure. */ + this = (AstMathMap *) this_object; + +/* Check the attribute name and clear the appropriate attribute. */ + +/* Seed. */ +/* ----- */ + if ( !strcmp( attrib, "seed" ) ) { + astClearSeed( this ); + +/* SimpFI. */ +/* ------- */ + } else if ( !strcmp( attrib, "simpfi" ) ) { + astClearSimpFI( this ); + +/* SimpIF. */ +/* ------- */ + } else if ( !strcmp( attrib, "simpif" ) ) { + astClearSimpIF( this ); + +/* If the attribute is not recognised, pass it on to the parent method + for further interpretation. */ + } else { + (*parent_clearattrib)( this_object, attrib, status ); + } +} + +static void CompileExpression( const char *method, const char *class, + const char *exprs, int nvar, const char *var[], + int **code, double **con, int *stacksize, int *status ) { +/* +* Name: +* CompileExpression + +* Purpose: +* Compile a mathematical expression. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* void CompileExpression( const char *method, const char *class, +* const char *exprs, int nvar, const char *var[], +* int **code, double **con, int *stacksize ) + +* Class Membership: +* MathMap member function. + +* Description: +* This function checks and compiles a mathematical expression. It +* produces a sequence of operation codes (opcodes) and a set of +* numerical constants which may subsequently be used to evaluate the +* expression on a push-down stack. + +* Parameters: +* method +* Pointer to a constant null-terminated character string +* containing the name of the method that invoked this function. +* This method name is used solely for constructing error messages. +* class +* Pointer to a constant null-terminated character string containing the +* class name of the Object being processed. This name is used solely +* for constructing error messages. +* exprs +* Pointer to a null-terminated string containing the expression +* to be compiled. This is case sensitive and should contain no white +* space. +* nvar +* The number of variable names defined for use in the expression. +* var +* An array of pointers (with "nvar" elements) to null-terminated +* strings. Each of these should contain a variable name which may +* appear in the expression. These strings are case sensitive and +* should contain no white space. +* code +* Address of a pointer which will be set to point at a dynamically +* allocated array of int containing the set of opcodes (cast to int) +* produced by this function. The first element of this array will +* contain a count of the number of opcodes which follow. +* +* The allocated space must be freed by the caller (using astFree) when +* no longer required. +* con +* Address of a pointer which will be set to point at a dynamically +* allocated array of double containing the set of constants +* produced by this function (this may be NULL if no constants are +* produced). +* +* The allocated space must be freed by the caller (using astFree) when +* no longer required. +* stacksize +* Pointer to an int in which to return the size of the push-down stack +* required to evaluate the expression using the returned opcodes and +* constants. + +* Algorithm: +* The function passes through the input expression searching for +* symbols. It looks for standard symbols (arithmetic operators, +* parentheses, function calls and delimiters) in the next part of the +* expression to be parsed, using identification information stored in +* the static "symbol" array. It ignores certain symbols, according to +* whether they appear to be operators or operands. The choice depends on +* what the previous symbol was; for instance, two operators may not +* occur in succession. Unary +/- operators are also ignored in +* situations where they are not permitted. +* +* If a standard symbol is found, it is passed to the ValidateSymbol +* function, which keeps track of the current level of parenthesis in the +* expression and of the number of arguments supplied to any (possibly +* nested) function calls. This function then accepts or rejects the +* symbol according to whether it is valid within the current context. An +* error is reported if it is rejected. +* +* If the part of the expression currently being parsed did not contain a +* standard symbol, an attempt is made to parse it first as a constant, +* then as a variable name. If either of these succeeds, an appropriate +* symbol number is added to the list of symbols identified so far, and a +* value is added to the list of constants - this is either the value of +* the constant itself, or the identification number of the variable. If +* the expression cannot be parsed, an error is reported. +* +* When the entire expression has been analysed as a sequence of symbols +* (and associated constants), the EvaluationSort function is +* invoked. This sorts the symbols into evaluation order, which is the +* order in which the associated operations must be performed on a +* push-down arithmetic stack to evaluate the expression. This routine +* also substitutes operation codes (defined in the "Oper" enum) for the +* symbol numbers and calculates the size of evaluation stack which will +* be required. + +* Notes: +* - A value of NULL will be returned for the "*code" and "*con" pointers +* and a value of zero will be returned for the "*stacksize" value if this +* function is invoked with the global error status set, or if it should +* fail for any reason. +*/ + +/* Local Variables: */ + double c; /* Value of parsed constant */ + int *argcount; /* Array of argument count information */ + int *opensym; /* Array of opening parenthesis information */ + int *symlist; /* Array of symbol indices */ + int found; /* Standard symbol identified? */ + int iend; /* Ending index in the expression string */ + int istart; /* Staring index in the expression string */ + int isym; /* Loop counter for symbols */ + int ivar; /* Index of variable name */ + int lpar; /* Parenthesis level */ + int ncon; /* Number of constants generated */ + int nsym; /* Number of symbols identified */ + int opernext; /* Next symbol an operator (from left)? */ + int size; /* Size of symbol matched */ + int sym; /* Index of symbol in static "symbol" array */ + int unarynext; /* Next symbol may be unary +/- ? */ + +/* Initialise. */ + *code = NULL; + *con = NULL; + *stacksize = 0; + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Further initialisation. */ + argcount = NULL; + lpar = 0; + ncon = 0; + nsym = 0; + opensym = NULL; + symlist = NULL; + sym = 0; + ivar = 0; + +/* The first symbol to be encountered must not look like an operator + from the left. It may be a unary + or - operator. */ + opernext = 0; + unarynext = 1; + +/* Search through the expression to classify each symbol which appears + in it. Stop when there are no more input characters or an error is + detected. */ + istart = 0; + for ( istart = 0; astOK && exprs[ istart ]; istart = iend + 1 ) { + +/* Compare each of the symbols in the symbol data with the next + section of the expression, looking for the longest symbol text which + will match. Stop if a NULL "text" value is found, which acts as the + end flag. */ + found = 0; + size = 0; + for ( isym = 0; symbol[ isym ].text; isym++ ) { + +/* Only consider symbols which have text associated with them and + which look like operators or operands from the left, according to the + setting of the "opernext" flag. Thus, if an operator or operand is + missing from the input expression, the next symbol will not be + identified, because it will be of the wrong type. Also exclude unary + +/- operators if they are out of context. */ + if ( symbol[ isym ].size && + ( symbol[ isym ].operleft == opernext ) && + ( !symbol[ isym ].unaryoper || unarynext ) ) { + +/* Test if the text of the symbol matches the expression at the + current position. If so, note that a match has been found. */ + if ( !strncmp( exprs + istart, symbol[ isym ].text, + (size_t) symbol[ isym ].size ) ) { + found = 1; + +/* If this symbol matches more characters than any previous symbol, + then store the symbol's index and note its size. */ + if ( symbol[ isym ].size > size ) { + sym = isym; + size = symbol[ isym ].size; + +/* Calculate the index of the last symbol character in the expression + string. */ + iend = istart + size - 1; + } + } + } + } + +/* If the symbol was identified as one of the standard symbols, then + validate it, updating the parenthesis level and argument count + information at the same time. */ + if ( found ) { + ValidateSymbol( method, class, exprs, iend, sym, &lpar, &argcount, + &opensym, &ncon, con, status ); + +/* If it was not one of the standard symbols, then check if the next + symbol was expected to be an operator. If so, then there is a missing + operator, so report an error. */ + } else { + if ( opernext ) { + astError( AST__MIOPR, + "%s(%s): Missing or invalid operator in the expression " + "\"%.*s\".", status, + method, class, istart + 1, exprs ); + +/* If the next symbol was expected to be an operand, then it may be a + constant, so try to parse it as one. */ + } else { + ParseConstant( method, class, exprs, istart, &iend, &c, status ); + if ( astOK ) { + +/* If successful, set the symbol number to "symbol_ldcon" (load + constant) and extend the "*con" array to accommodate a new + constant. Check for errors. */ + if ( iend >= istart ) { + sym = symbol_ldcon; + *con = astGrow( *con, ncon + 1, sizeof( double ) ); + if ( astOK ) { + +/* Append the constant to the "*con" array. */ + ( *con )[ ncon++ ] = c; + } + +/* If the symbol did not parse as a constant, then it may be a + variable name, so try to parse it as one. */ + } else { + ParseVariable( method, class, exprs, istart, nvar, var, + &ivar, &iend, status ); + if ( astOK ) { + +/* If successful, set the symbol to "symbol_ldvar" (load variable) and + extend the "*con" array to accommodate a new constant. Check for + errors. */ + if ( ivar != -1 ) { + sym = symbol_ldvar; + *con = astGrow( *con, ncon + 1, sizeof( double ) ); + if ( astOK ) { + +/* Append the variable identification number as a constant to the + "*con" array. */ + ( *con )[ ncon++ ] = (double) ivar; + } + +/* If the expression did not parse as a variable name, then there is a + missing operand in the expression, so report an error. */ + } else { + astError( AST__MIOPA, + "%s(%s): Missing or invalid operand in the " + "expression \"%.*s\".", status, + method, class, istart + 1, exprs ); + } + } + } + } + } + } + +/* If there has been no error, then the next symbol in the input + expression has been identified and is valid. */ + if ( astOK ) { + +/* Decide whether the next symbol should look like an operator or an + operand from the left. This is determined by the nature of the symbol + just identified (seen from the right) - two operands or two operators + cannot be adjacent. */ + opernext = !symbol[ sym ].operright; + +/* Also decide whether the next symbol may be a unary +/- operator, + according to the "unarynext" symbol data entry for the symbol just + identified. */ + unarynext = symbol[ sym ].unarynext; + +/* Extend the "symlist" array to accommodate the symbol just + identified. Check for errors. */ + symlist = astGrow( symlist, nsym + 1, sizeof( int ) ); + if ( astOK ) { + +/* Append the symbol's index to the end of this list. */ + symlist[ nsym++ ] = sym; + } + } + } + +/* If there has been no error, check the final context after + identifying all the symbols... */ + if ( astOK ) { + +/* If an operand is still expected, then there is an unsatisfied + operator on the end of the expression, so report an error. */ + if ( !opernext ) { + astError( AST__MIOPA, + "%s(%s): Missing or invalid operand in the expression " + "\"%s\".", status, + method, class, exprs ); + +/* If the final parenthesis level is positive, then there is a missing + right parenthesis, so report an error. */ + } else if ( lpar > 0 ) { + astError( AST__MRPAR, + "%s(%s): Missing right parenthesis in the expression " + "\"%s\".", status, + method, class, exprs ); + } + } + +/* Sort the symbols into evaluation order to produce output opcodes. */ + EvaluationSort( *con, nsym, symlist, code, stacksize, status ); + +/* Free any memory used as workspace. */ + if ( argcount ) argcount = astFree( argcount ); + if ( opensym ) opensym = astFree( opensym ); + if ( symlist ) symlist = astFree( symlist ); + +/* If OK, re-allocate the "*con" array to have the correct size (since + astGrow may have over-allocated space). */ + if ( astOK && *con ) { + *con = astRealloc( *con, sizeof( double ) * (size_t) ncon ); + } + +/* If an error occurred, free any allocated memory and reset the + output values. */ + if ( !astOK ) { + *code = astFree( *code ); + *con = astFree( *con ); + *stacksize = 0; + } +} + +static void CompileMapping( const char *method, const char *class, + int nin, int nout, + int nfwd, const char *fwdfun[], + int ninv, const char *invfun[], + int ***fwdcode, int ***invcode, + double ***fwdcon, double ***invcon, + int *fwdstack, int *invstack, int *status ) { +/* +* Name: +* CompileMapping + +* Purpose: +* Compile the transformation functions for a MathMap. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* void CompileMapping( const char *method, const char *class, +* int nin, int nout, +* int nfwd, const char *fwdfun[], +* int ninv, const char *invfun[], +* int ***fwdcode, int ***invcode, +* double ***fwdcon, double ***invcon, +* int *fwdstack, int *invstack, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* This function checks and compiles the transformation functions required +* to create a MathMap. It produces sequences of operation codes (opcodes) +* and numerical constants which may subsequently be used to evaluate the +* functions on a push-down stack. + +* Parameters: +* method +* Pointer to a constant null-terminated character string +* containing the name of the method that invoked this function. +* This method name is used solely for constructing error messages. +* class +* Pointer to a constant null-terminated character string containing the +* class name of the Object being processed. This name is used solely +* for constructing error messages. +* nin +* Number of input variables for the MathMap. +* nout +* Number of output variables for the MathMap. +* nfwd +* The number of forward transformation functions being supplied. +* This must be at least equal to "nout". +* fwdfun +* Pointer to an array, with "nfwd" elements, of pointers to null +* terminated strings which contain each of the forward transformation +* functions. These must be in lower case and should contain no white +* space. +* ninv +* The number of inverse transformation functions being supplied. +* This must be at least equal to "nin". +* invfun +* Pointer to an array, with "ninv" elements, of pointers to null +* terminated strings which contain each of the inverse transformation +* functions. These must be in lower case and should contain no white +* space. +* fwdcode +* Address in which to return a pointer to an array (with "nfwd" +* elements) of pointers to arrays of int containing the set of opcodes +* (cast to int) for each forward transformation function. The number +* of opcodes produced for each function is given by the first element +* of the opcode array. +* +* Both the returned array of pointers, and the arrays of int to which +* they point, will be stored in dynamically allocated memory and should +* be freed by the caller (using astFree) when no longer required. +* +* If the right hand sides (including the "=" sign) of all the supplied +* functions are absent, then this indicates an undefined transformation +* and the returned pointer value will be NULL. An error results if +* an "=" sign is present but no expression follows it. +* invcode +* Address in which to return a pointer to an array (with "ninv" +* elements) of pointers to arrays of int containing the set of opcodes +* (cast to int) for each inverse transformation function. The number +* of opcodes produced for each function is given by the first element +* of the opcode array. +* +* Both the returned array of pointers, and the arrays of int to which +* they point, will be stored in dynamically allocated memory and should +* be freed by the caller (using astFree) when no longer required. +* +* If the right hand sides (including the "=" sign) of all the supplied +* functions are absent, then this indicates an undefined transformation +* and the returned pointer value will be NULL. An error results if +* an "=" sign is present but no expression follows it. +* fwdcon +* Address in which to return a pointer to an array (with "nfwd" +* elements) of pointers to arrays of double containing the set of +* constants for each forward transformation function. +* +* Both the returned array of pointers, and the arrays of double to which +* they point, will be stored in dynamically allocated memory and should +* be freed by the caller (using astFree) when no longer required. Note +* that any of the pointers to the arrays of double may be NULL if no +* constants are associated with a particular function. +* +* If the forward transformation is undefined, then the returned pointer +* value will be NULL. +* invcon +* Address in which to return a pointer to an array (with "ninv" +* elements) of pointers to arrays of double containing the set of +* constants for each inverse transformation function. +* +* Both the returned array of pointers, and the arrays of double to which +* they point, will be stored in dynamically allocated memory and should +* be freed by the caller (using astFree) when no longer required. Note +* that any of the pointers to the arrays of double may be NULL if no +* constants are associated with a particular function. +* +* If the inverse transformation is undefined, then the returned pointer +* value will be NULL. +* fwdstack +* Pointer to an int in which to return the size of the push-down stack +* required to evaluate the forward transformation functions. +* invstack +* Pointer to an int in which to return the size of the push-down stack +* required to evaluate the inverse transformation functions. +* status +* Pointer to the inherited status variable. + +* Notes: +* - A value of NULL will be returned for the "*fwdcode", "*invcode", +* "*fwdcon" and "*invcon" pointers and a value of zero will be returned +* for the "*fwdstack" and "*invstack" values if this function is invoked +* with the global error status set, or if it should fail for any reason. +*/ + +/* Local Variables: */ + char **exprs; /* Pointer to array of expressions */ + char **var; /* Pointer to array of variable names */ + const char **strings; /* Pointer to temporary array of strings */ + int ifun; /* Loop counter for functions */ + int nvar; /* Number of variables to extract */ + int stacksize; /* Required stack size */ + +/* Initialise. */ + *fwdcode = NULL; + *invcode = NULL; + *fwdcon = NULL; + *invcon = NULL; + *fwdstack = 0; + *invstack = 0; + nvar = 0; + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Further initialisation. */ + exprs = NULL; + var = NULL; + +/* Compile the forward transformation. */ +/* ----------------------------------- */ +/* Allocate space for an array of pointers to the functions from which + we will extract variable names. */ + strings = astMalloc( sizeof( char * ) * (size_t) ( nin + nfwd ) ); + +/* Fill the first elements of this array with pointers to the inverse + transformation functions ("nin" in number) which yield the final input + values. These will have the names of the input variables on their left + hand sides. */ + if ( astOK ) { + nvar = 0; + for ( ifun = ninv - nin; ifun < ninv; ifun++ ) { + strings[ nvar++ ] = invfun[ ifun ]; + } + +/* Fill the remaining elements of the array with pointers to the + forward transformation functions. These will have the names of any + intermediate variables plus the final output variables on their left + hand sides. */ + for ( ifun = 0; ifun < nfwd; ifun++ ) strings[ nvar++ ] = fwdfun[ ifun ]; + +/* Extract the variable names from the left hand sides of these + functions and check them for validity and absence of duplication. */ + ExtractVariables( method, class, nvar, strings, nin, nout, nfwd, ninv, 1, + &var, status ); + } + +/* Free the temporary array of string pointers. */ + strings = astFree( strings ); + +/* Extract the expressions from the right hand sides of the forward + transformation functions. */ + ExtractExpressions( method, class, nfwd, fwdfun, 1, &exprs, status ); + +/* If OK, and the forward transformation is defined, then allocate and + initialise space for an array of pointers to the opcodes for each + expression and, similarly, for the constants for each expression. */ + if ( astOK && exprs ) { + MALLOC_POINTER_ARRAY( *fwdcode, int *, nfwd ) + MALLOC_POINTER_ARRAY( *fwdcon, double *, nfwd ) + +/* If OK, loop to compile each of the expressions, storing pointers to + the resulting opcodes and constants in the arrays allocated above. On + each loop, we make progressively more of the variable names in "var" + visible to the compilation function. This ensures that each expression + can only use variables which have been defined earlier. */ + if ( astOK ) { + for ( ifun = 0; ifun < nfwd; ifun++ ) { + CompileExpression( method, class, exprs[ ifun ], + nin + ifun, (const char **) var, + &( *fwdcode )[ ifun ], &( *fwdcon )[ ifun ], + &stacksize, status ); + +/* If an error occurs, then report contextual information and quit. */ + if ( !astOK ) { + astError( astStatus, + "Error in forward transformation function %d.", status, + ifun + 1 ); + break; + } + +/* If OK, calculate the maximum evaluation stack size required by any + of the expressions. */ + *fwdstack = ( *fwdstack > stacksize ) ? *fwdstack : stacksize; + } + } + } + +/* Free the memory containing the extracted expressions and variables. */ + FREE_POINTER_ARRAY( exprs, nfwd ) + FREE_POINTER_ARRAY( var, nvar ) + +/* Compile the inverse transformation. */ +/* ----------------------------------- */ +/* Allocate space for an array of pointers to the functions from which + we will extract variable names. */ + strings = astMalloc( sizeof( char * ) * (size_t) ( nout + ninv ) ); + +/* Fill the first elements of this array with pointers to the forward + transformation functions ("nout" in number) which yield the final + output values. These will have the names of the output variables on + their left hand sides. */ + if ( astOK ) { + nvar = 0; + for ( ifun = nfwd - nout; ifun < nfwd; ifun++ ) { + strings[ nvar++ ] = fwdfun[ ifun ]; + } + +/* Fill the remaining elements of the array with pointers to the + inverse transformation functions. These will have the names of any + intermediate variables plus the final input variables on their left + hand sides. */ + for ( ifun = 0; ifun < ninv; ifun++ ) strings[ nvar++ ] = invfun[ ifun ]; + +/* Extract the variable names from the left hand sides of these + functions and check them for validity and absence of duplication. */ + ExtractVariables( method, class, nvar, strings, nin, nout, nfwd, ninv, 0, + &var, status ); + } + +/* Free the temporary array of string pointers. */ + strings = astFree( strings ); + +/* Extract the expressions from the right hand sides of the inverse + transformation functions. */ + ExtractExpressions( method, class, ninv, invfun, 0, &exprs, status ); + +/* If OK, and the forward transformation is defined, then allocate and + initialise space for an array of pointers to the opcodes for each + expression and, similarly, for the constants for each expression. */ + if ( astOK && exprs ) { + MALLOC_POINTER_ARRAY( *invcode, int *, ninv ) + MALLOC_POINTER_ARRAY( *invcon, double *, ninv ) + +/* If OK, loop to compile each of the expressions, storing pointers to + the resulting opcodes and constants in the arrays allocated above. On + each loop, we make progressively more of the variable names in "var" + visible to the compilation function. This ensures that each expression + can only use variables which have been defined earlier. */ + if ( astOK ) { + for ( ifun = 0; ifun < ninv; ifun++ ) { + CompileExpression( method, class, exprs[ ifun ], + nout + ifun, (const char **) var, + &( *invcode )[ ifun ], &( *invcon )[ ifun ], + &stacksize, status ); + +/* If an error occurs, then report contextual information and quit. */ + if ( !astOK ) { + astError( astStatus, + "Error in inverse transformation function %d.", status, + ifun + 1 ); + break; + } + +/* If OK, calculate the maximum evaluation stack size required by any + of the expressions. */ + *invstack = ( *invstack > stacksize ) ? *invstack : stacksize; + } + } + } + +/* Free the memory containing the extracted expressions and variables. */ + FREE_POINTER_ARRAY( exprs, ninv ) + FREE_POINTER_ARRAY( var, nvar ) + +/* If an error occurred, then free all remaining allocated memory and + reset the output values. */ + if ( !astOK ) { + FREE_POINTER_ARRAY( *fwdcode, nfwd ) + FREE_POINTER_ARRAY( *invcode, ninv ) + FREE_POINTER_ARRAY( *fwdcon, nfwd ) + FREE_POINTER_ARRAY( *invcon, ninv ) + *fwdstack = 0; + *invstack = 0; + } +} + +static int DefaultSeed( const Rcontext *context, int *status ) { +/* +* Name: +* DefaultSeed + +* Purpose: +* Generate an unpredictable seed for a random number generator. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* int DefaultSeed( Rcontext *context, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* On each invocation this function returns an integer value which is +* highly unpredictable. This value may be used as a default seed for the +* random number generator associated with a MathMap, so that it +* generates a different sequence on each occasion. + +* Parameters: +* context +* Pointer to the random number generator context associated with +* the MathMap. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* The unpredictable integer. + +* Notes: +* - This function does not perform error checking and will execute even +* if the global error status is set. +*/ + +/* Local Constants: */ + const int nwarm = 5; /* Number of warm-up iterations */ + const long int a = 8121L; /* Constants for random number generator... */ + const long int c = 28411L; + const long int m = 134456L; + +/* Local Variables; */ + int iwarm; /* Loop counter for warm-up iterations */ + static long init = 0; /* Local initialisation performed? */ + static long int rand; /* Local random integer */ + unsigned long int bits; /* Bit pattern for producing result */ + +/* On the first invocation, initialise a local random number generator + to a value derived by combining bit patterns obtained from the system + clock and the processor time used. The result needs to be positive and + lie in the range 0 to "m-1". */ + LOCK_MUTEX5 + if ( !init ) { + rand = (long int) ( ( (unsigned long int) time( NULL ) ^ + (unsigned long int) clock() ) % + (unsigned long int) m ); + +/* These values will typically only change in their least significant + bits between programs run successively, but by using the bit pattern + as a seed, we ensure that these differences are rapidly propagated to + other bits. To hasten this process, we "warm up" the local generator + with a few iterations. This is a quick and dirty generator using + constants from Press et al. (Numerical recipes). */ + for ( iwarm = 0; iwarm < nwarm; iwarm++ ) { + rand = ( rand * a + c ) % m; + } + +/* Note that this initialisation has been performed. */ + init = 1; + } + UNLOCK_MUTEX5 + +/* Generate a new bit pattern from the system time. Apart from the + first invocation, this will be a different time to that used above. */ + bits = (unsigned long int) time( NULL ); + +/* Mask in a pattern derived from the CPU time used. */ + bits ^= (unsigned long int) clock(); + +/* The system time may change quite slowly (e.g. every second), so + also mask in the address of the random number generator context + supplied. This makes the seed depend on which MathMap is in use. */ + bits ^= (unsigned long int) context; + +/* Now mask in the last random integer produced by the random number + generator whose context has been supplied. This makes the seed depend + on the MathMap's past use of random numbers. */ + bits ^= (unsigned long int) context->random_int; + +/* Finally, in order to produce different seeds when this function is + invoked twice in rapid succession on the same object (with no + intermediate processing), we also mask in a pseudo-random value + generated here. Generate the next local random integer. */ + rand = ( rand * a + c ) % m; + +/* We then scale this value to give an integer in the range 0 to + ULONG_MAX and mask the corresponding bit pattern into our seed. */ + bits ^= (unsigned long int) ( ( (double) rand / (double) ( m - 1UL ) ) * + ( ( (double) ULONG_MAX + 1.0 ) * + ( 1.0 - DBL_EPSILON ) ) ); + +/* Return the integer value of the seed (which may involve discarding + some unwanted bits). */ + return (int) bits; +} + +static int Equal( AstObject *this_object, AstObject *that_object, int *status ) { +/* +* Name: +* Equal + +* Purpose: +* Test if two MathMaps are equivalent. + +* Type: +* Private function. + +* Synopsis: +* #include "mapping.h" +* int Equal( AstObject *this, AstObject *that, int *status ) + +* Class Membership: +* MathMap member function (over-rides the astEqual protected +* method inherited from the Object class). + +* Description: +* This function returns a boolean result (0 or 1) to indicate whether +* two MathMaps are equivalent. + +* Parameters: +* this +* Pointer to the first Object (a MathMap). +* that +* Pointer to the second Object. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* One if the MathMaps are equivalent, zero otherwise. + +* Notes: +* - The two MathMaps are considered equivalent if the combination of +* the first in series with the inverse of the second simplifies to a +* UnitMap. +* - A value of zero will be returned if this function is invoked +* with the global status set, or if it should fail for any reason. +*/ + +/* Local Variables: */ + AstMathMap *that; /* Pointer to the second MathMap structure */ + AstMathMap *this; /* Pointer to the first MathMap structure */ + double **that_con; /* Lists of constants from "that" */ + double **this_con; /* Lists of constants from "this" */ + int **that_code; /* Lists of opcodes from "that" */ + int **this_code; /* Lists of opcodes from "this" */ + int code; /* Opcode value */ + int icode; /* Opcode index */ + int icon; /* Constant index */ + int ifun; /* Function index */ + int ncode; /* No. of opcodes for current "this" function */ + int ncode_that; /* No. of opcodes for current "that" function */ + int nin; /* Number of inputs */ + int nout; /* Number of outputs */ + int pass; /* Check fwd or inv */ + int result; /* Result value to return */ + int that_nfun; /* Number of functions from "that" */ + int this_nfun; /* Number of functions from "this" */ + +/* Initialise. */ + result = 0; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Obtain pointers to the two MathMap structures. */ + this = (AstMathMap *) this_object; + that = (AstMathMap *) that_object; + +/* Check the second object is a MathMap. We know the first is a + MathMap since we have arrived at this implementation of the virtual + function. */ + if( astIsAMathMap( that ) ) { + +/* Check they have the same number of inputs and outputs */ + nin = astGetNin( this ); + nout = astGetNout( this ); + if( astGetNout( that ) == nout && astGetNin( that ) == nin ) { + +/* Assume equality. */ + result = 1; + +/* The first pass through this next loop compares forward functions, and + the second pass compares inverse functions. */ + for( pass = 0; pass < 2 && result; pass++ ) { + +/* On the first pass, get pointers to the lists of opcodes and constants for + the effective forward transformations (taking into account the value + of the Invert attribute), together with the number of such functions. */ + if( pass == 0 ) { + if( !astGetInvert( this ) ) { + this_code = this->fwdcode; + this_con = this->fwdcon; + this_nfun = this->nfwd; + } else { + this_code = this->invcode; + this_con = this->invcon; + this_nfun = this->ninv; + } + + if( !astGetInvert( that ) ) { + that_code = that->fwdcode; + that_con = that->fwdcon; + that_nfun = that->nfwd; + } else { + that_code = that->invcode; + that_con = that->invcon; + that_nfun = that->ninv; + } + +/* On the second pass, get pointers to the lists of opcodes and constants for + the effective inverse transformations, together with the number of such + functions. */ + } else { + + if( astGetInvert( this ) ) { + this_code = this->fwdcode; + this_con = this->fwdcon; + this_nfun = this->nfwd; + } else { + this_code = this->invcode; + this_con = this->invcon; + this_nfun = this->ninv; + } + + if( astGetInvert( that ) ) { + that_code = that->fwdcode; + that_con = that->fwdcon; + that_nfun = that->nfwd; + } else { + that_code = that->invcode; + that_con = that->invcon; + that_nfun = that->ninv; + } + } + +/* Check that "this" and "that" have the same number of functions */ + if( that_nfun != this_nfun ) result = 0; + +/* Loop round each function. */ + for( ifun = 0; ifun < this_nfun && result; ifun++ ) { + +/* The first element in the opcode array is the number of subsequent + opcodes. Obtain and compare these counts. */ + ncode = this_code ? this_code[ ifun ][ 0 ] : 0; + ncode_that = that_code ? that_code[ ifun ][ 0 ] : 0; + if( ncode != ncode_that ) result = 0; + +/* Compare the following opcodes. Some opcodes consume constants from the + list of constants associated with the MathMap. Compare the constants + for such opcodes. */ + icon = 0; + for( icode = 0; icode < ncode && result; icode++ ){ + code = this_code[ ifun ][ icode ]; + if( that_code[ ifun ][ icode ] != code ) { + result = 0; + + } else if( code == OP_LDCON || + code == OP_LDVAR || + code == OP_MAX || + code == OP_MIN ) { + + if( this_con[ ifun ][ icon ] != + that_con[ ifun ][ icon ] ) { + result = 0; + } else { + icon++; + } + } + } + } + } + } + } + +/* If an error occurred, clear the result value. */ + if ( !astOK ) result = 0; + +/* Return the result, */ + return result; +} + +static void EvaluateFunction( Rcontext *rcontext, int npoint, + const double **ptr_in, const int *code, + const double *con, int stacksize, double *out, int *status ) { +/* +* Name: +* EvaluateFunction + +* Purpose: +* Evaluate a compiled function. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* void EvaluateFunction( Rcontext *rcontext, int npoint, +* const double **ptr_in, const int *code, +* const double *con, int stacksize, double *out, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* This function implements a "virtual machine" which executes operations +* on an arithmetic stack in order to evaluate transformation functions. +* Each operation is specified by an input operation code (opcode) and +* results in the execution of a vector operation on a stack. The final +* result, after executing all the supplied opcodes, is returned as a +* vector. +* +* This function detects arithmetic errors (such as overflow and division +* by zero) and propagates any "bad" coordinate values, including those +* present in the input, to the output. + +* Parameters: +* npoint +* The number of points to be transformd (i.e. the size of the vector +* of values on which operations are to be performed). +* ptr_in +* Pointer to an array of pointers to arrays of double (with "npoint" +* elements). These arrays should contain the input coordinate values, +* such that coordinate number "coord" for point number "point" can be +* found in "ptr_in[coord][point]". +* code +* Pointer to an array of int containing the set of opcodes (cast to int) +* for the operations to be performed. The first element of this array +* should contain a count of the number of opcodes which follow. +* con +* Pointer to an array of double containing the set of constants required +* to evaluate the function (this may be NULL if no constants are +* required). +* stacksize +* The size of the stack required to evaluate the expression using the +* opcodes and constants supplied. This value should be calculated during +* expression compilation. +* out +* Pointer to an array of double (with "npoint" elements) in which to +* return the vector of result values. +* status +* Pointer to the inherited status variable. +*/ + +/* Local Constants: */ + const int bits = /* Number of bits in an unsigned long */ + sizeof( unsigned long ) * CHAR_BIT; + const double eps = /* Smallest number subtractable from 2.0 */ + 2.0 * DBL_EPSILON; + const double scale = /* 2.0 raised to the power "bits" */ + ldexp( 1.0, bits ); + const double scale1 = /* 2.0 raised to the power "bits-1" */ + scale * 0.5; + const double rscale = /* Reciprocal scale factor */ + 1.0 / scale; + const double rscale1 = /* Reciprocal initial scale factor */ + 1.0 / scale1; + const int nblock = /* Number of blocks of bits to process */ + ( sizeof( double ) + sizeof( unsigned long ) - 1 ) / + sizeof( unsigned long ); + const unsigned long signbit = /* Mask for extracting sign bit */ + 1UL << ( bits - 1 ); + +/* Local Variables: */ + double **stack; /* Array of pointers to stack elements */ + double *work; /* Pointer to stack workspace */ + double *xv1; /* Pointer to first argument vector */ + double *xv2; /* Pointer to second argument vector */ + double *xv3; /* Pointer to third argument vector */ + double *xv; /* Pointer to sole argument vector */ + double *y; /* Pointer to result */ + double *yv; /* Pointer to result vector */ + double abs1; /* Absolute value (temporary variable) */ + double abs2; /* Absolute value (temporary variable) */ + double frac1; /* First (maybe normalised) fraction */ + double frac2; /* Second (maybe normalised) fraction */ + double frac; /* Sole normalised fraction */ + double newexp; /* New power of 2 exponent value */ + double ran; /* Random number */ + double result; /* Function result value */ + double unscale; /* Factor for removing scaling */ + double value; /* Value to be assigned to stack vector */ + double x1; /* First argument value */ + double x2; /* Second argument value */ + double x3; /* Third argument value */ + double x; /* Sole argument value */ + int expon1; /* First power of 2 exponent */ + int expon2; /* Second power of 2 exponent */ + int expon; /* Sole power of 2 exponent */ + int iarg; /* Loop counter for arguments */ + int iblock; /* Loop counter for blocks of bits */ + int icode; /* Opcode value */ + int icon; /* Counter for number of constants used */ + int istk; /* Loop counter for stack elements */ + int ivar; /* Input variable number */ + int narg; /* Number of function arguments */ + int ncode; /* Number of opcodes to process */ + int point; /* Loop counter for stack vector elements */ + int sign; /* Argument is non-negative? */ + int tos; /* Top of stack index */ + static double d2r; /* Degrees to radians conversion factor */ + static double log2; /* Natural logarithm of 2.0 */ + static double pi; /* Value of PI */ + static double r2d; /* Radians to degrees conversion factor */ + static double rsafe_sq; /* Reciprocal of "safe_sq" */ + static double safe_sq; /* Huge value that can safely be squared */ + static int init = 0; /* Initialisation performed? */ + unsigned long b1; /* Block of bits from first argument */ + unsigned long b2; /* Block of bits from second argument */ + unsigned long b; /* Block of bits for result */ + unsigned long neg; /* Result is negative? (sign bit) */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* If this is the first invocation of this function, then initialise + constant values. */ + LOCK_MUTEX2 + if ( !init ) { + +/* Trigonometrical conversion factors. */ + pi = acos( -1.0 ); + r2d = 180.0 / pi; + d2r = pi / 180.0; + +/* Natural logarithm of 2.0. */ + log2 = log( 2.0 ); + +/* This value must be safe to square without producing overflow, yet + large enough that adding or subtracting 1.0 from the square makes no + difference. We also need its reciprocal. */ + safe_sq = 0.9 * sqrt( DBL_MAX ); + rsafe_sq = 1.0 / safe_sq; + +/* Note that initialisation has been performed. */ + init = 1; + } + UNLOCK_MUTEX2 + +/* Allocate space for an array of pointers to elements of the + workspace stack (each stack element being an array of double). */ + stack = astMalloc( sizeof( double * ) * (size_t) stacksize ); + +/* Allocate space for the stack itself. */ + work = astMalloc( sizeof( double ) * + (size_t) ( npoint * ( stacksize - 1 ) ) ); + +/* If OK, then initialise the stack pointer array to identify the + start of each vector on the stack. The first element points at the + output array (in which the result will be accumulated), while other + elements point at successive vectors within the workspace allocated + above. */ + if ( astOK ) { + stack[ 0 ] = out; + for ( istk = 1; istk < stacksize; istk++ ) { + stack[ istk ] = work + ( istk - 1 ) * npoint; + } + +/* Define stack operations. */ +/* ======================== */ +/* We now define a set of macros for performing vector operations on + elements of the stack. Each is in the form of a "case" block for + execution in response to the appropriate operation code (opcode). */ + +/* Zero-argument operation. */ +/* ------------------------ */ +/* This macro performs a zero-argument operation, which results in the + insertion of a new vector on to the stack. */ +#define ARG_0(oper,setup,function) \ +\ +/* Test for the required opcode value. */ \ + case oper: \ +\ +/* Perform any required initialisation. */ \ + {setup;} \ +\ +/* Increment the top of stack index and obtain a pointer to the new stack \ + element (vector). */ \ + yv = stack[ ++tos ]; \ +\ +/* Loop to access each vector element, obtaining a pointer to it. */ \ + for ( point = 0; point < npoint; point++ ) { \ + y = yv + point; \ +\ +/* Perform the processing, which results in assignment to this element. */ \ + {function;} \ + } \ +\ +/* Break out of the "case" block. */ \ + break; + +/* One-argument operation. */ +/* ----------------------- */ +/* This macro performs a one-argument operation, which processes the + top stack element without changing the stack size. */ +#define ARG_1(oper,function) \ +\ +/* Test for the required opcode value. */ \ + case oper: \ +\ +/* Obtain a pointer to the top stack element (vector). */ \ + xv = stack[ tos ]; \ +\ +/* Loop to access each vector element, obtaining its value and \ + checking that it is not bad. */ \ + for ( point = 0; point < npoint; point++ ) { \ + if ( ( x = xv[ point ] ) != AST__BAD ) { \ +\ +/* Also obtain a pointer to the element. */ \ + y = xv + point; \ +\ +/* Perform the processing, which uses the element's value and then \ + assigns the result to this element. */ \ + {function;} \ + } \ + } \ +\ +/* Break out of the "case" block. */ \ + break; + +/* One-argument boolean operation. */ +/* ------------------------------- */ +/* This macro is similar in function to ARG_1 above, except that no + checks are made for bad argument values. It is intended for use with + boolean functions where bad values are handled explicitly. */ +#define ARG_1B(oper,function) \ +\ +/* Test for the required opcode value. */ \ + case oper: \ +\ +/* Obtain a pointer to the top stack element (vector). */ \ + xv = stack[ tos ]; \ +\ +/* Loop to access each vector element, obtaining the argument value \ + and a pointer to the element. */ \ + for ( point = 0; point < npoint; point++ ) { \ + x = xv[ point ]; \ + y = xv + point; \ +\ +/* Perform the processing, which uses the element's value and then \ + assigns the result to this element. */ \ + {function;} \ + } \ +\ +/* Break out of the "case" block. */ \ + break; + +/* Two-argument operation. */ +/* ----------------------- */ +/* This macro performs a two-argument operation, which processes the + top two stack elements and produces a single result, resulting in the + stack size decreasing by one. In this case, we first define a macro + without the "case" block statements present. */ +#define DO_ARG_2(function) \ +\ +/* Obtain pointers to the top two stack elements (vectors), decreasing \ + the top of stack index by one. */ \ + xv2 = stack[ tos-- ]; \ + xv1 = stack[ tos ]; \ +\ +/* Loop to access each vector element, obtaining the value of the \ + first argument and checking that it is not bad. */ \ + for ( point = 0; point < npoint; point++ ) { \ + if ( ( x1 = xv1[ point ] ) != AST__BAD ) { \ +\ +/* Also obtain a pointer to the element which is to receive the \ + result. */ \ + y = xv1 + point; \ +\ +/* Obtain the value of the second argument, again checking that it is \ + not bad. */ \ + if ( ( x2 = xv2[ point ] ) != AST__BAD ) { \ +\ +/* Perform the processing, which uses the two argument values and then \ + assigns the result to the appropriate top of stack element. */ \ + {function;} \ +\ +/* If the second argument was bad, so is the result. */ \ + } else { \ + *y = AST__BAD; \ + } \ + } \ + } + +/* This macro simply wraps the one above up in a "case" block. */ +#define ARG_2(oper,function) \ + case oper: \ + DO_ARG_2(function) \ + break; + +/* Two-argument boolean operation. */ +/* ------------------------------- */ +/* This macro is similar in function to ARG_2 above, except that no + checks are made for bad argument values. It is intended for use with + boolean functions where bad values are handled explicitly. */ +#define ARG_2B(oper,function) \ +\ +/* Test for the required opcode value. */ \ + case oper: \ +\ +/* Obtain pointers to the top two stack elements (vectors), decreasing \ + the top of stack index by one. */ \ + xv2 = stack[ tos-- ]; \ + xv1 = stack[ tos ]; \ +\ +/* Loop to access each vector element, obtaining the value of both \ + arguments and a pointer to the element which is to receive the \ + result. */ \ + for ( point = 0; point < npoint; point++ ) { \ + x1 = xv1[ point ]; \ + x2 = xv2[ point ]; \ + y = xv1 + point; \ +\ +/* Perform the processing, which uses the two argument values and then \ + assigns the result to the appropriate top of stack element. */ \ + {function;} \ + } \ +\ +/* Break out of the "case" block. */ \ + break; + +/* Three-argument boolean operation. */ +/* --------------------------------- */ +/* This macro is similar in function to ARG_2B above, except that it + takes three values of the stack and puts one back. It performs no + checks for bad values. */ +#define ARG_3B(oper,function) \ +\ +/* Test for the required opcode value. */ \ + case oper: \ +\ +/* Obtain pointers to the top three stack elements (vectors), decreasing \ + the top of stack index by two. */ \ + xv3 = stack[ tos-- ]; \ + xv2 = stack[ tos-- ]; \ + xv1 = stack[ tos ]; \ +\ +/* Loop to access each vector element, obtaining the value of all 3 \ + arguments and a pointer to the element which is to receive the \ + result. */ \ + for ( point = 0; point < npoint; point++ ) { \ + x1 = xv1[ point ]; \ + x2 = xv2[ point ]; \ + x3 = xv3[ point ]; \ + y = xv1 + point; \ +\ +/* Perform the processing, which uses the three argument values and then \ + assigns the result to the appropriate top of stack element. */ \ + {function;} \ + } \ +\ +/* Break out of the "case" block. */ \ + break; + +/* Define arithmetic operations. */ +/* ============================= */ +/* We now define macros for performing some of the arithmetic + operations we will require in a "safe" way - i.e. trapping numerical + problems such as overflow and invalid arguments and translating them + into the AST__BAD value. */ + +/* Absolute value. */ +/* --------------- */ +/* This is just shorthand. */ +#define ABS(x) ( ( (x) >= 0.0 ) ? (x) : -(x) ) + +/* Integer part. */ +/* ------------- */ +/* This implements rounding towards zero without involving conversion + to an integer (which could overflow). */ +#define INT(x) ( ( (x) >= 0.0 ) ? floor( (x) ) : ceil( (x) ) ) + +/* Trap maths overflow. */ +/* -------------------- */ +/* This macro calls a C maths library function and checks for overflow + in the result. */ +#define CATCH_MATHS_OVERFLOW(function) \ + ( \ +\ +/* Clear the "errno" value. */ \ + errno = 0, \ +\ +/* Evaluate the function. */ \ + result = (function), \ +\ +/* Check if "errno" and the returned result indicate overflow and \ + return the appropriate result. */ \ + ( ( errno == ERANGE ) && ( ABS( result ) == HUGE_VAL ) ) ? AST__BAD : \ + result \ + ) + +/* Trap maths errors. */ +/* ------------------ */ +/* This macro is similar to the one above, except that it also checks + for domain errors (i.e. invalid argument values). */ +#define CATCH_MATHS_ERROR(function) \ + ( \ +\ +/* Clear the "errno" value. */ \ + errno = 0, \ +\ +/* Evaluate the function. */ \ + result = (function), \ +\ +/* Check if "errno" and the returned result indicate a domain error or \ + overflow and return the appropriate result. */ \ + ( ( errno == EDOM ) || \ + ( ( errno == ERANGE ) && ( ABS( result ) == HUGE_VAL ) ) ) ? \ + AST__BAD : result \ + ) + +/* Tri-state boolean OR. */ +/* --------------------- */ +/* This evaluates a boolean OR using tri-state logic. For example, + "a||b" may evaluate to 1 if "a" is bad but "b" is non-zero, so that + the normal rules of bad value propagation do not apply. */ +#define TRISTATE_OR(x1,x2) \ +\ +/* Test if the first argument is bad. */ \ + ( (x1) == AST__BAD ) ? ( \ +\ +/* If so, test the second argument. */ \ + ( ( (x2) == 0.0 ) || ( (x2) == AST__BAD ) ) ? AST__BAD : 1.0 \ + ) : ( \ +\ +/* Test if the second argument is bad. */ \ + ( (x2) == AST__BAD ) ? ( \ +\ +/* If so, test the first argument. */ \ + ( (x1) == 0.0 ) ? AST__BAD : 1.0 \ +\ +/* If neither argument is bad, use the normal OR operator. */ \ + ) : ( \ + ( (x1) != 0.0 ) || ( (x2) != 0.0 ) \ + ) \ + ) + +/* Tri-state boolean AND. */ +/* ---------------------- */ +/* This evaluates a boolean AND using tri-state logic. */ +#define TRISTATE_AND(x1,x2) \ +\ +/* Test if the first argument is bad. */ \ + ( (x1) == AST__BAD ) ? ( \ +\ +/* If so, test the second argument. */ \ + ( (x2) != 0.0 ) ? AST__BAD : 0.0 \ + ) : ( \ +\ +/* Test if the second argument is bad. */ \ + ( (x2) == AST__BAD ) ? ( \ +\ +/* If so, test the first argument. */ \ + ( (x1) != 0.0 ) ? AST__BAD : 0.0 \ +\ +/* If neither argument is bad, use the normal AND operator. */ \ + ) : ( \ + ( (x1) != 0.0 ) && ( (x2) != 0.0 ) \ + ) \ + ) + +/* Safe addition. */ +/* -------------- */ +/* This macro performs addition while avoiding possible overflow. */ +#define SAFE_ADD(x1,x2) ( \ +\ +/* Test if the first argument is non-negative. */ \ + ( (x1) >= 0.0 ) ? ( \ +\ +/* If so, then we can perform addition if the second argument is \ + non-positive. Otherwise, we must calculate the most positive safe \ + second argument value that can be added and test for this (the test \ + itself is safe against overflow). */ \ + ( ( (x2) <= 0.0 ) || ( ( (DBL_MAX) - (x1) ) >= (x2) ) ) ? ( \ +\ +/* Perform addition if it is safe, otherwise return AST__BAD. */ \ + (x1) + (x2) \ + ) : ( \ + AST__BAD \ + ) \ +\ +/* If the first argument is negative, then we can perform addition if \ + the second argument is non-negative. Otherwise, we must calculate the \ + most negative second argument value that can be added and test for \ + this (the test itself is safe against overflow). */ \ + ) : ( \ + ( ( (x2) >= 0.0 ) || ( ( (DBL_MAX) + (x1) ) >= -(x2) ) ) ? ( \ +\ +/* Perform addition if it is safe, otherwise return AST__BAD. */ \ + (x1) + (x2) \ + ) : ( \ + AST__BAD \ + ) \ + ) \ +) + +/* Safe subtraction. */ +/* ----------------- */ +/* This macro performs subtraction while avoiding possible overflow. */ +#define SAFE_SUB(x1,x2) ( \ +\ +/* Test if the first argument is non-negative. */ \ + ( (x1) >= 0.0 ) ? ( \ +\ +/* If so, then we can perform subtraction if the second argument is \ + also non-negative. Otherwise, we must calculate the most negative safe \ + second argument value that can be subtracted and test for this (the \ + test itself is safe against overflow). */ \ + ( ( (x2) >= 0.0 ) || ( ( (DBL_MAX) - (x1) ) >= -(x2) ) ) ? ( \ +\ +/* Perform subtraction if it is safe, otherwise return AST__BAD. */ \ + (x1) - (x2) \ + ) : ( \ + AST__BAD \ + ) \ +\ +/* If the first argument is negative, then we can perform subtraction \ + if the second argument is non-positive. Otherwise, we must calculate \ + the most positive second argument value that can be subtracted and \ + test for this (the test itself is safe against overflow). */ \ + ) : ( \ + ( ( (x2) <= 0.0 ) || ( ( (DBL_MAX) + (x1) ) >= (x2) ) ) ? ( \ +\ +/* Perform subtraction if it is safe, otherwise return AST__BAD. */ \ + (x1) - (x2) \ + ) : ( \ + AST__BAD \ + ) \ + ) \ +) + +/* Safe multiplication. */ +/* -------------------- */ +/* This macro performs multiplication while avoiding possible overflow. */ +#define SAFE_MUL(x1,x2) ( \ +\ +/* Multiplication is safe if the absolute value of either argument is \ + unity or less. Otherwise, we must use the first argument to calculate \ + the maximum absolute value that the second argument may have and test \ + for this (the test itself is safe against overflow). */ \ + ( ( ( abs1 = ABS( (x1) ) ) <= 1.0 ) || \ + ( ( abs2 = ABS( (x2) ) ) <= 1.0 ) || \ + ( ( (DBL_MAX) / abs1 ) >= abs2 ) ) ? ( \ +\ +/* Perform multiplication if it is safe, otherwise return AST__BAD. */ \ + (x1) * (x2) \ + ) : ( \ + AST__BAD \ + ) \ +) + +/* Safe division. */ +/* -------------- */ +/* This macro performs division while avoiding possible overflow. */ +#define SAFE_DIV(x1,x2) ( \ +\ +/* Division is unsafe if the second argument is zero. Otherwise, it is \ + safe if the abolute value of the second argument is unity or \ + more. Otherwise, we must use the second argument to calculate the \ + maximum absolute value that the first argument may have and test for \ + this (the test itself is safe against overflow). */ \ + ( ( (x2) != 0.0 ) && \ + ( ( ( abs2 = ABS( (x2) ) ) >= 1.0 ) || \ + ( ( (DBL_MAX) * abs2 ) >= ABS( (x1) ) ) ) ) ? ( \ +\ +/* Perform division if it is safe, otherwise return AST__BAD. */ \ + (x1) / (x2) \ + ) : ( \ + AST__BAD \ + ) \ +) + +/* Bit-shift operation. */ +/* -------------------- */ +/* This macro shifts the bits in a double value a specified number of + places to the left, which simply corresponds to multiplying by the + appropriate power of two. */ +#define SHIFT_BITS(x1,x2) ( \ +\ +/* Decompose the value into a normalised fraction and a power of 2. */ \ + frac = frexp( (x1), &expon ), \ +\ +/* Calculate the new power of 2 which should apply after the shift, \ + rounding towards zero to give an integer value. */ \ + newexp = INT( (x2) ) + (double) expon, \ +\ +/* If the new exponent is too negative to convert to an integer, then \ + the result must underflow to zero. */ \ + ( newexp < (double) -INT_MAX ) ? ( \ + 0.0 \ +\ +/* Otherwise, if it is too positive to convert to an integer, then the \ + result must overflow, unless the normalised fraction is zero. */ \ + ) : ( ( newexp > (double) INT_MAX ) ? ( \ + ( frac == 0.0 ) ? 0.0 : AST__BAD \ +\ +/* Otherwise, convert the new exponent to an integer and apply \ + it. Trap any overflow which may still occur. */ \ + ) : ( \ + CATCH_MATHS_OVERFLOW( ldexp( frac, (int) newexp ) ) \ + ) ) \ +) + +/* Two-argument bit-wise boolean operation. */ +/* ---------------------------------------- */ +/* This macro expands to code which performs a bit-wise boolean + operation on a pair of arguments and assigns the result to the + variable "result". It operates on floating point (double) values, + which are regarded as if they are fixed-point binary numbers with + negative values expressed in twos-complement notation. This means that + it delivers the same results for integer values as the normal + (integer) C bit-wise operations. However, it will also operate on the + fraction bits of floating point numbers. It also offers greater + precision (the first 53 or so significant bits of the result being + preserved for typical IEEE floating point implementations). */ +#define BIT_OPER(oper,x1,x2) \ +\ +/* Convert each argument to a normalised fraction in the range \ + [0.5,1.0) and a power of two exponent, removing any sign \ + information. */ \ + frac1 = frexp( ABS( (x1) ), &expon1 ); \ + frac2 = frexp( ABS( (x2) ), &expon2 ); \ +\ +/* Set "expon" to be the larger of the two exponents. If the two \ + exponents are not equal, divide the fraction with the smaller exponent \ + by 2 to the power of the exponent difference. This gives both \ + fractions the same effective exponent (although one of them may no \ + longer be normalised). Note that overflow is avoided because all \ + numbers remain less than 1.0, but underflow may occur. */ \ + expon = expon1; \ + if ( expon2 > expon1 ) { \ + expon = expon2; \ + frac1 = ldexp( frac1, expon1 - expon ); \ + } else if ( expon1 > expon2 ) { \ + frac2 = ldexp( frac2, expon2 - expon ); \ + } \ +\ +/* If either of the original arguments is negative, we now subtract \ + the corresponding fraction from 2.0. If we think of the fraction as \ + represented in fixed-point binary notation, this corresponds to \ + converting negative numbers into the twos-complement form normally used \ + for integers (the sign bit being the bit with value 1) instead \ + of having a separate sign bit as for floating point numbers. \ +\ + Note that one of the fractions may have underflowed during the \ + scaling above. In that case (if the original argument was negative), \ + we must subtract the value "eps" (= 2.0 * DBL_EPSILON) from 2.0 \ + instead, so that we produce the largest number less than 2.0. In \ + twos-complement notation this represents the smallest possible \ + negative number and corresponds to extending the sign bit of the \ + original number up into more significant bits. This causes all bits to \ + be set as we require (rather than all being clear if the underflow \ + is simply ignored). */ \ + if ( (x1) < 0.0 ) frac1 = 2.0 - ( ( frac1 > eps ) ? frac1 : eps ); \ + if ( (x2) < 0.0 ) frac2 = 2.0 - ( ( frac2 > eps ) ? frac2 : eps ); \ +\ +/* We now extract the bits from the fraction values into integer \ + variables so that we may perform bit-wise operations on them. However, \ + since a double may be longer than any available integer, we may \ + have to handle several successive blocks of bits individually. */ \ +\ +/* Extract the first block of bits by scaling by the required power of \ + 2 to shift the required bits to the left of the binary point. Then \ + extract the integer part. Note that this initial shift is one bit less \ + than the number of bits in an unsigned long, because we have \ + introduced an extra sign bit. */ \ + frac1 *= scale1; \ + frac2 *= scale1; \ + b1 = (unsigned long) frac1; \ + b2 = (unsigned long) frac2; \ +\ +/* Perform the required bit-wise operation on the extracted blocks of \ + bits. */ \ + b = b1 oper b2; \ +\ +/* Extract the sign bit from this initial result. This determines \ + whether the final result bit pattern should represent a negative \ + floating point number. */ \ + neg = b & signbit; \ +\ +/* Initialise the floating point result by setting it to the integer \ + result multipled by the reciprocal of the scale factor used to shift \ + the bits above. This returns the result bits to their correct \ + significance. */ \ + unscale = rscale1; \ + result = (double) b * unscale; \ +\ +/* We now loop to extract and process further blocks of bits (if \ + present). The number of blocks is determined by the relative lengths \ + of a double and an unsigned long. In practice, some bits of the double \ + will be used by its exponent, so the last block may be incomplete and \ + will simply be padded with zeros. */ \ + for ( iblock = 1; iblock < nblock; iblock++ ) { \ +\ +/* Subtract the integer part (which has already been processed) from \ + each fraction, to leave the bits which remain to be processed. Then \ + multiply by a scale factor to shift the next set of bits to the left \ + of the binary point. This time, we use as many bits as will fit into \ + an unsigned long. */ \ + frac1 = ( frac1 - (double) b1 ) * scale; \ + frac2 = ( frac2 - (double) b2 ) * scale; \ +\ +/* Extract the integer part, which contains the required bits. */ \ + b1 = (unsigned long) frac1; \ + b2 = (unsigned long) frac2; \ +\ +/* Perform the required bit-wise operation on the extracted blocks of \ + bits. */ \ + b = b1 oper b2; \ +\ +/* Update the result floating point value by adding the new integer \ + result multiplied by a scale factor to return the bits to their \ + original significance. */ \ + unscale *= rscale; \ + result += (double) b * unscale; \ + } \ +\ +/* If the (normalised fraction) result represents a negative number, \ + then subtract 2.0 from it (equivalent to subtracting it from 2 and \ + negating the result). This converts back to using a separate sign bit \ + instead of twos-complement notation. */ \ + if ( neg ) result -= 2.0; \ +\ +/* Scale by the required power of 2 to remove the initial \ + normalisation applied and assign the result to the "result" \ + variable. */ \ + result = ldexp( result, expon ) + +/* Gaussian random number. */ +/* ----------------------- */ +/* This macro expands to code which assigns a pseudo-random value to + the "result" variable. The value is drawn from a Gaussian distribution + with mean "x1" and standard deviation "ABS(x2)". */ +#define GAUSS(x1,x2) \ +\ +/* Loop until a satisfactory result is obtained. */ \ + do { \ +\ +/* Obtain a value drawn from a standard Gaussian distribution. */ \ + ran = Gauss( rcontext, status ); \ +\ +/* Multiply by "ABS(x2)", trapping possible overflow. */ \ + result = ABS( (x2) ); \ + result = SAFE_MUL( ran, result ); \ +\ +/* If OK, add "x1", again trapping possible overflow. */ \ + if ( result != AST__BAD ) result = SAFE_ADD( result, (x1) ); \ +\ +/* Continue generating values until one is found which does not cause \ + overflow. */ \ + } while ( result == AST__BAD ); + +/* Implement the stack-based arithmetic. */ +/* ===================================== */ +/* Initialise the top of stack index and constant counter. */ + tos = -1; + icon = 0; + +/* Determine the number of opcodes to be processed and loop to process + them, executing the appropriate "case" block for each one. */ + ncode = code[ 0 ]; + for ( icode = 1; icode <= ncode; icode++ ) { + switch ( (Oper) code[ icode ] ) { + +/* Ignore any null opcodes (which shouldn't occur). */ + case OP_NULL: break; + +/* Otherwise, perform the required vector operation on the stack... */ + +/* User-supplied constants and variables. */ +/* -------------------------------------- */ +/* Loading a constant involves incrementing the constant count and + assigning the next constant's value to the top of stack element. */ + ARG_0( OP_LDCON, value = con[ icon++ ], *y = value ) + +/* Loading a variable involves obtaining the variable's index by + consuming a constant (as above), and then copying the variable's + values into the top of stack element. */ + ARG_0( OP_LDVAR, ivar = (int) ( con[ icon++ ] + 0.5 ), + *y = ptr_in[ ivar ][ point ] ) + +/* System constants. */ +/* ----------------- */ +/* Loading a "bad" value simply means assigning AST__BAD to the top of + stack element. */ + ARG_0( OP_LDBAD, ;, *y = AST__BAD ) + +/* The following load constants associated with the (double) floating + point representation into the top of stack element. */ + ARG_0( OP_LDDIG, ;, *y = (double) AST__DBL_DIG ) + ARG_0( OP_LDEPS, ;, *y = DBL_EPSILON ) + ARG_0( OP_LDMAX, ;, *y = DBL_MAX ) + ARG_0( OP_LDMAX10E, ;, *y = (double) DBL_MAX_10_EXP ) + ARG_0( OP_LDMAXE, ;, *y = (double) DBL_MAX_EXP ) + ARG_0( OP_LDMDIG, ;, *y = (double) DBL_MANT_DIG ) + ARG_0( OP_LDMIN, ;, *y = DBL_MIN ) + ARG_0( OP_LDMIN10E, ;, *y = (double) DBL_MIN_10_EXP ) + ARG_0( OP_LDMINE, ;, *y = (double) DBL_MIN_EXP ) + ARG_0( OP_LDRAD, ;, *y = (double) FLT_RADIX ) + ARG_0( OP_LDRND, ;, *y = (double) FLT_ROUNDS ) + +/* Mathematical constants. */ +/* ----------------------- */ +/* The following load mathematical constants into the top of stack + element. */ + ARG_0( OP_LDE, value = exp( 1.0 ), *y = value ) + ARG_0( OP_LDPI, ;, *y = pi ) + +/* Functions with one argument. */ +/* ---------------------------- */ +/* The following simply evaluate a function of the top of stack + element and assign the result to the same element. */ + ARG_1( OP_ABS, *y = ABS( x ) ) + ARG_1( OP_ACOS, *y = ( ABS( x ) <= 1.0 ) ? + acos( x ) : AST__BAD ) + ARG_1( OP_ACOSD, *y = ( ABS( x ) <= 1.0 ) ? + acos( x ) * r2d : AST__BAD ) + ARG_1( OP_ACOSH, *y = ( x < 1.0 ) ? AST__BAD : + ( ( x > safe_sq ) ? log( x ) + log2 : + log( x + sqrt( x * x - 1.0 ) ) ) ) + ARG_1( OP_ACOTH, *y = ( ABS( x ) <= 1.0 ) ? AST__BAD : + 0.5 * ( log( ( x + 1.0 ) / + ( x - 1.0 ) ) ) ) + ARG_1( OP_ACSCH, *y = ( ( x == 0.0 ) ? AST__BAD : + ( sign = ( x >= 0.0 ), x = ABS( x ), + ( sign ? 1.0 : -1.0 ) * + ( ( x < rsafe_sq ) ? log2 - log( x ) : + ( x = 1.0 / x, + log( x + sqrt( x * x + 1.0 ) ) ) ) ) ) ) + ARG_1( OP_ASECH, *y = ( ( x <= 0 ) || ( x > 1.0 ) ) ? AST__BAD : + ( ( x < rsafe_sq ) ? log2 - log( x ) : + ( x = 1.0 / x, + log( x + sqrt( x * x - 1.0 ) ) ) ) ) + ARG_1( OP_ASIN, *y = ( ABS( x ) <= 1.0 ) ? + asin( x ) : AST__BAD ) + ARG_1( OP_ASIND, *y = ( ABS( x ) <= 1.0 ) ? + asin( x ) * r2d : AST__BAD ) + ARG_1( OP_ASINH, *y = ( sign = ( x >= 0.0 ), x = ABS( x ), + ( sign ? 1.0 : -1.0 ) * + ( ( x > safe_sq ) ? log( x ) + log2 : + log( x + sqrt( x * x + 1.0 ) ) ) ) ) + ARG_1( OP_ATAN, *y = atan( x ) ) + ARG_1( OP_ATAND, *y = atan( x ) * r2d ) + ARG_1( OP_ATANH, *y = ( ABS( x ) >= 1.0 ) ? AST__BAD : + 0.5 * ( log( ( 1.0 + x ) / + ( 1.0 - x ) ) ) ) + ARG_1( OP_CEIL, *y = ceil( x ) ) + ARG_1( OP_COS, *y = cos( x ) ) + ARG_1( OP_COSD, *y = cos( x * d2r ) ) + ARG_1( OP_COSH, *y = CATCH_MATHS_OVERFLOW( cosh( x ) ) ) + ARG_1( OP_COTH, *y = ( x = tanh( x ), SAFE_DIV( 1.0, x ) ) ) + ARG_1( OP_CSCH, *y = ( x = CATCH_MATHS_OVERFLOW( sinh( x ) ), + ( x == AST__BAD ) ? + 0.0 : SAFE_DIV( 1.0, x ) ) ) + ARG_1( OP_EXP, *y = CATCH_MATHS_OVERFLOW( exp( x ) ) ) + ARG_1( OP_FLOOR, *y = floor( x ) ) + ARG_1( OP_INT, *y = INT( x ) ) + ARG_1B( OP_ISBAD, *y = ( x == AST__BAD ) ) + ARG_1( OP_LOG, *y = ( x > 0.0 ) ? log( x ) : AST__BAD ) + ARG_1( OP_LOG10, *y = ( x > 0.0 ) ? log10( x ) : AST__BAD ) + ARG_1( OP_NINT, *y = ( x >= 0 ) ? + floor( x + 0.5 ) : ceil( x - 0.5 ) ) + ARG_1( OP_POISS, *y = Poisson( rcontext, x, status ) ) + ARG_1( OP_SECH, *y = ( x = CATCH_MATHS_OVERFLOW( cosh( x ) ), + ( x == AST__BAD ) ? 0.0 : 1.0 / x ) ) + ARG_1( OP_SIN, *y = sin( x ) ) + ARG_1( OP_SINC, *y = ( x == 0.0 ) ? 1.0 : sin( x ) / x ) + ARG_1( OP_SIND, *y = sin( x * d2r ) ) + ARG_1( OP_SINH, *y = CATCH_MATHS_OVERFLOW( sinh( x ) ) ) + ARG_1( OP_SQR, *y = SAFE_MUL( x, x ) ) + ARG_1( OP_SQRT, *y = ( x >= 0.0 ) ? sqrt( x ) : AST__BAD ) + ARG_1( OP_TAN, *y = CATCH_MATHS_OVERFLOW( tan( x ) ) ) + ARG_1( OP_TAND, *y = tan( x * d2r ) ) + ARG_1( OP_TANH, *y = tanh( x ) ) + +/* Functions with two arguments. */ +/* ----------------------------- */ +/* These evaluate a function of the top two entries on the stack. */ + ARG_2( OP_ATAN2, *y = atan2( x1, x2 ) ) + ARG_2( OP_ATAN2D, *y = atan2( x1, x2 ) * r2d ) + ARG_2( OP_DIM, *y = ( x1 > x2 ) ? x1 - x2 : 0.0 ) + ARG_2( OP_GAUSS, GAUSS( x1, x2 ); *y = result ) + ARG_2( OP_MOD, *y = ( x2 != 0.0 ) ? + fmod( x1, x2 ) : AST__BAD ) + ARG_2( OP_POW, *y = CATCH_MATHS_ERROR( pow( x1, x2 ) ) ) + ARG_2( OP_RAND, ran = Rand( rcontext, status ); + *y = x1 * ran + x2 * ( 1.0 - ran ); ) + ARG_2( OP_SIGN, *y = ( ( x1 >= 0.0 ) == ( x2 >= 0.0 ) ) ? + x1 : -x1 ) + +/* Functions with three arguments. */ +/* ------------------------------- */ +/* These evaluate a function of the top three entries on the stack. */ + ARG_3B( OP_QIF, *y = ( ( x1 ) ? ( x2 ) : ( x3 ) ) ) + + +/* Functions with variable numbers of arguments. */ +/* --------------------------------------------- */ +/* These operations take a variable number of arguments, the actual + number being determined by consuming a constant. We then loop to + perform a 2-argument operation on the stack (as above) the required + number of times. */ + case OP_MAX: + narg = (int) ( con[ icon++ ] + 0.5 ); + for ( iarg = 0; iarg < ( narg - 1 ); iarg++ ) { + DO_ARG_2( *y = ( x1 >= x2 ) ? x1 : x2 ) + } + break; + case OP_MIN: + narg = (int) ( con[ icon++ ] + 0.5 ); + for ( iarg = 0; iarg < ( narg - 1 ); iarg++ ) { + DO_ARG_2( *y = ( x1 <= x2 ) ? x1 : x2 ) + } + break; + +/* Unary arithmetic operators. */ +/* --------------------------- */ + ARG_1( OP_NEG, *y = -x ) + +/* Unary boolean operators. */ +/* ------------------------ */ + ARG_1( OP_NOT, *y = ( x == 0.0 ) ) + +/* Binary arithmetic operators. */ +/* ---------------------------- */ + ARG_2( OP_ADD, *y = SAFE_ADD( x1, x2 ) ) + ARG_2( OP_SUB, *y = SAFE_SUB( x1, x2 ) ) + ARG_2( OP_MUL, *y = SAFE_MUL( x1, x2 ) ) + ARG_2( OP_DIV , *y = SAFE_DIV( x1, x2 ) ) + +/* Bit-shift operators. */ +/* -------------------- */ + ARG_2( OP_SHFTL, *y = SHIFT_BITS( x1, x2 ) ) + ARG_2( OP_SHFTR, *y = SHIFT_BITS( x1, -x2 ) ) + +/* Relational operators. */ +/* --------------------- */ + ARG_2( OP_EQ, *y = ( x1 == x2 ) ) + ARG_2( OP_GE, *y = ( x1 >= x2 ) ) + ARG_2( OP_GT, *y = ( x1 > x2 ) ) + ARG_2( OP_LE, *y = ( x1 <= x2 ) ) + ARG_2( OP_LT, *y = ( x1 < x2 ) ) + ARG_2( OP_NE, *y = ( x1 != x2 ) ) + +/* Bit-wise operators. */ +/* ------------------- */ + ARG_2( OP_BITOR, BIT_OPER( |, x1, x2 ); *y = result ) + ARG_2( OP_BITXOR, BIT_OPER( ^, x1, x2 ); *y = result ) + ARG_2( OP_BITAND, BIT_OPER( &, x1, x2 ); *y = result ) + +/* Binary boolean operators. */ +/* ------------------------- */ + ARG_2B( OP_AND, *y = TRISTATE_AND( x1, x2 ) ) + ARG_2( OP_EQV, *y = ( ( x1 != 0.0 ) == ( x2 != 0.0 ) ) ) + ARG_2B( OP_OR, *y = TRISTATE_OR( x1, x2 ) ) + ARG_2( OP_XOR, *y = ( ( x1 != 0.0 ) != ( x2 != 0.0 ) ) ) + } + } + } + +/* When all opcodes have been processed, the result of the function + evaluation will reside in the lowest stack entry - i.e. the output + array. */ + +/* Free the workspace arrays. */ + work = astFree( work ); + stack = astFree( stack ); + +/* Undefine macros local to this function. */ +#undef ARG_0 +#undef ARG_1 +#undef ARG_1B +#undef DO_ARG_2 +#undef ARG_2 +#undef ARG_2B +#undef ABS +#undef INT +#undef CATCH_MATHS_OVERFLOW +#undef CATCH_MATHS_ERROR +#undef TRISTATE_OR +#undef TRISTATE_AND +#undef SAFE_ADD +#undef SAFE_SUB +#undef SAFE_MUL +#undef SAFE_DIV +#undef SHIFT_BITS +#undef BIT_OPER +#undef GAUSS +} + +static void EvaluationSort( const double con[], int nsym, int symlist[], + int **code, int *stacksize, int *status ) { +/* +* Name: +* EvaluationSort + +* Purpose: +* Perform an evaluation-order sort on parsed expression symbols. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* void EvaluationSort( const double con[], int nsym, int symlist[], +* int **code, int *stacksize, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* This function sorts a sequence of numbers representing symbols +* identified in an expression. The symbols (i.e. the expression syntax) +* must have been fully validated beforehand, as no validation is +* performed here. +* +* The symbols are sorted into the order in which corresponding +* operations must be performed on a push-down arithmetic stack in order +* to evaluate the expression. Operation codes (opcodes), as defined in +* the "Oper" enum, are then substituted for the symbol numbers. + +* Parameters: +* con +* Pointer to an array of double containing the set of constants +* generated while parsing the expression (these are required in order +* to determine the number of arguments associated with functions which +* take a variable number of arguments). +* nsym +* The number of symbols identified while parsing the expression. +* symlist +* Pointer to an array of int, with "nsym" elements. On entry, this +* should contain the indices in the static "symbol" array of the +* symbols identified while parsing the expression. On exit, the +* contents are undefined. +* code +* Address of a pointer which will be set to point at a dynamically +* allocated array of int containing the set of opcodes (cast to int) +* produced by this function. The first element of this array will +* contain a count of the number of opcodes which follow. +* +* The allocated space must be freed by the caller (using astFree) when +* no longer required. +* stacksize +* Pointer to an int in which to return the size of the push-down stack +* required to evaluate the expression using the returned opcodes. +* status +* Pointer to the inherited status variable. + +* Notes: +* - A value of NULL will be returned for the "*code" pointer and a value +* of zero will be returned for the "*stacksize" value if this function is +* invoked with the global error status set, or if it should fail for any +* reason. +*/ + +/* Local Variables: */ + int flush; /* Flush parenthesised symbol sequence? */ + int icon; /* Input constant counter */ + int isym; /* Input symbol counter */ + int ncode; /* Number of opcodes generated */ + int nstack; /* Evaluation stack size */ + int push; /* Push a new symbol on to stack? */ + int sym; /* Variable for symbol number */ + int tos; /* Top of sort stack index */ + +/* Initialise */ + *code = NULL; + *stacksize = 0; + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Further initialisation. */ + flush = 0; + icon = 0; + isym = 0; + ncode = 0; + nstack = 0; + tos = -1; + +/* Loop to generate output opcodes until the sort stack is empty and + there are no further symbols to process, or an error is detected. */ + while ( astOK && ( ( tos > -1 ) || ( isym < nsym ) ) ) { + +/* Decide whether to push a symbol on to the sort stack (which + "diverts" it so that higher-priority symbols can be output), or to pop + the top symbol off the sort stack and send it to the output + stream... */ + +/* We must push a symbol on to the sort stack if the stack is + currently empty. */ + if ( tos == -1 ) { + push = 1; + +/* We must pop the top symbol off the sort stack if there are no more + input symbols to process. */ + } else if ( isym >= nsym ) { + push = 0; + +/* If the sort stack is being flushed to complete the evaluation of a + parenthesised expression, then the top symbol (which will be the + opening parenthesis or function call) must be popped. This is only + done once, so reset the "flush" flag before the next loop. */ + } else if ( flush ) { + push = 0; + flush = 0; + +/* In all other circumstances, we must push a symbol on to the sort + stack if its evaluation priority (seen from the left) is higher than + that of the current top of stack symbol (seen from the right). This + means it will eventually be sent to the output stream ahead of the + current top of stack symbol. */ + } else { + push = ( symbol[ symlist[ isym ] ].leftpriority > + symbol[ symlist[ tos ] ].rightpriority ); + } + +/* If a symbol is being pushed on to the sort stack, then get the next + input symbol which is to be used. */ + if ( push ) { + sym = symlist[ isym++ ]; + +/* If the symbol decreases the parenthesis level (a closing + parenthesis), then all the sort stack entries down to the symbol which + opened the current level of parenthesis (the matching opening + parenthesis or function call) will already have been sent to the + output stream as a consequence of the evaluation priority defined for + a closing parenthesis in the symbol data. The opening parenthesis (or + function call) must next be flushed from the sort stack, so set the + "flush" flag which is interpreted on the next loop. Ignore the current + symbol, which cancels with the opening parenthesis on the stack. */ + if ( symbol[ sym ].parincrement < 0 ) { + flush = 1; + +/* All other symbols are pushed on to the sort stack. The stack + occupies that region of the "symlist" array from which the input + symbol numbers have already been extracted. */ + } else { + symlist[ ++tos ] = sym; + } + +/* If a symbol is being popped from the top of the sort stack, then + the top of stack entry is transferred to the output stream. Obtain the + symbol number from the stack. Increment the local constant counter if + the associated operation will use a constant. */ + } else { + sym = symlist[ tos-- ]; + icon += ( ( sym == symbol_ldvar ) || ( sym == symbol_ldcon ) ); + +/* If the output symbol does not represent a "null" operation, + increase the size of the output opcode array to accommodate it, + checking for errors. Note that we allocate one extra array element + (the first) which will eventually hold a count of all the opcodes + generated. */ + if ( symbol[ sym ].opcode != OP_NULL ) { + *code = astGrow( *code, ncode + 2, sizeof( int ) ); + if ( astOK ) { + +/* Append the new opcode to the end of this array. */ + ( *code )[ ++ncode ] = (int) symbol[ sym ].opcode; + +/* Increment/decrement the counter representing the stack size + required for evaluation of the expression. If the symbol is a + function with a variable number of arguments (indicated by a negative + "nargs" entry in the symbol data table), then the change in stack size + must be determined from the argument number stored in the constant + table. */ + if ( symbol[ sym ].nargs >= 0 ) { + nstack += symbol[ sym ].stackincrement; + } else { + nstack -= (int) ( con[ icon++ ] + 0.5 ) - 1; + } + +/* Note the maximum size of the stack. */ + *stacksize = ( nstack > *stacksize ) ? nstack : *stacksize; + } + } + } + } + +/* If no "*code" array has been allocated, then allocate one simply to + store the number of opcodes generated, i.e. zero (this shouldn't + normally happen as this represents an invalid expression). */ + if ( !*code ) *code = astMalloc( sizeof( int ) ); + +/* If no error has occurred, store the count of opcodes generated in + the first element of the "*code" array and re-allocate the array to + its final size (since astGrow may have over-allocated space). */ + if ( astOK ) { + ( *code )[ 0 ] = ncode; + *code = astRealloc( *code, sizeof( int ) * (size_t) ( ncode + 1 ) ); + } + +/* If an error occurred, free any memory that was allocated and reset + the output values. */ + if ( !astOK ) { + *code = astFree( *code ); + *stacksize = 0; + } +} + +static void ExtractExpressions( const char *method, const char *class, + int nfun, const char *fun[], int forward, + char ***exprs, int *status ) { +/* +* Name: +* ExtractExpressions + +* Purpose: +* Extract and validate expressions. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* void ExtractExpressions( const char *method, const char *class, +* int nfun, const char *fun[], int forward, +* char ***exprs, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* This function extracts expressions from the right hand sides of a set +* of functions. These expressions are then validated to check that they +* are either all present, or all absent (absence indicating an undefined +* transformation). An error is reported if anything is found to be +* wrong. +* +* Note that the syntax of the expressions is not checked by this function +* (i.e. they are not compiled). + +* Parameters: +* method +* Pointer to a constant null-terminated character string +* containing the name of the method that invoked this function. +* This method name is used solely for constructing error messages. +* class +* Pointer to a constant null-terminated character string containing the +* class name of the Object being processed. This name is used solely +* for constructing error messages. +* nfun +* The number of functions to be analysed. +* fun +* Pointer to an array, with "nfun" elements, of pointers to null +* terminated strings which contain each of the functions. These +* strings should contain no white space. +* forward +* A non-zero value indicates the the MathMap's forward transformation +* functions are being processed, while a zero value indicates processing +* of the inverse transformation functions. This value is used solely for +* constructing error messages. +* exprs +* Address in which to return a pointer to an array (with "nfun" +* elements) of pointers to null terminated strings containing the +* extracted expressions (i.e. this returns an array of strings). +* +* Both the returned array of pointers, and the strings to which they +* point, will be stored in dynamically allocated memory and should +* be freed by the caller (using astFree) when no longer required. +* +* If the right hand sides (including the "=" sign) of all the supplied +* functions are absent, then this indicates an undefined transformation +* and the returned pointer value will be NULL. An error results if +* an "=" sign is present but no expression follows it. +* status +* Pointer to the inherited status variable. + +* Notes: +* - A NULL value will be returned for "*exprs" if this function is +* invoked with the global error status set, or if it should fail for +* any reason. +*/ + +/* Local Variables: */ + char *ex; /* Pointer to start of expression string */ + int ifun; /* Loop counter for functions */ + int iud; /* Index of first undefined function */ + int nud; /* Number of undefined expressions */ + +/* Initialise. */ + *exprs = NULL; + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Further initialisation. */ + nud = 0; + iud = 0; + +/* Allocate and initialise memory for the returned array of pointers. */ + MALLOC_POINTER_ARRAY( *exprs, char *, nfun ) + +/* Loop to inspect each function in turn. */ + if ( astOK ) { + for ( ifun = 0; ifun < nfun; ifun++ ) { + +/* Search for the first "=" sign. */ + if ( ( ex = strchr( fun[ ifun ], '=' ) ) ) { + +/* If found, and there are more characters after the "=" sign, then + find the length of the expression which follows. Allocate a string to + hold this expression, storing its pointer in the array allocated + above. Check for errors. */ + if ( *++ex ) { + ( *exprs )[ ifun ] = astMalloc( strlen( ex ) + (size_t) 1 ); + if ( !astOK ) break; + +/* If OK, extract the expression string. */ + (void) strcpy( ( *exprs )[ ifun ], ex ); + +/* If an "=" sign was found but there are no characters following it, + then there is a missing right hand side to a function, so report an + error and quit. */ + } else { + astError( AST__NORHS, + "%s(%s): Missing right hand side in expression: " + "\"%s\".", status, + method, class, fun[ ifun ] ); + astError( astStatus, + "Error in %s transformation function %d.", status, + forward ? "forward" : "inverse", ifun + 1 ); + break; + } + +/* If no "=" sign was found, then the transformation may be undefined, + in which case each function should only contain a variable name. Count + the number of times this happens and record the index of the first + instance. */ + } else { + nud++; + if ( nud == 1 ) iud = ifun; + } + } + } + +/* Either all functions should have an "=" sign (in which case the + transformation is defined), or none of them should have (in which case + it is undefined). If some do and some don't, then report an error, + citing the first instance of a missing "=" sign. */ + if ( astOK && ( nud != 0 ) && ( nud != nfun ) ) { + astError( AST__NORHS, + "%s(%s): Missing right hand side in function: \"%s\".", status, + method, class, fun[ iud ] ); + astError( astStatus, + "Error in %s transformation function %d.", status, + forward ? "forward" : "inverse", iud + 1 ); + } + +/* If an error occurred, or all the expressions were absent, then free any + allocated memory and reset the output value. */ + if ( !astOK || nud ) { + FREE_POINTER_ARRAY( *exprs, nfun ) + } +} + +static void ExtractVariables( const char *method, const char *class, + int nfun, const char *fun[], + int nin, int nout, int nfwd, int ninv, + int forward, char ***var, int *status ) { +/* +* Name: +* ExtractVariables + +* Purpose: +* Extract and validate variable names. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* void ExtractVariables( const char *method, const char *class, +* int nfun, const char *fun[], +* int nin, int nout, int nfwd, int ninv, +* int forward, char ***var, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* This function extracts variable names from the left hand sides of a +* set of transformation functions belonging to a MathMap. These variable +* names are then validated to check for correct syntax and no +* duplication. An error is reported if anything is wrong with the +* variable names obtained. + +* Parameters: +* method +* Pointer to a constant null-terminated character string +* containing the name of the method that invoked this function. +* This method name is used solely for constructing error messages. +* class +* Pointer to a constant null-terminated character string containing the +* class name of the Object being processed. This name is used solely +* for constructing error messages. +* nfun +* The number of functions to be analysed. +* fun +* Pointer to an array, with "nfun" elements, of pointers to null +* terminated strings which contain each of the functions. These strings +* are case sensitive and should contain no white space. +* +* The first elements of this array should point to the functions that +* define the primary input/output variables (depending on direction). +* These should be followed by any functions which define intermediate +* variables (taken from the set of functions which transform in the +* opposite direction to the first ones). +* nin +* Number of input variables for the MathMap. +* nout +* Number of output variables for the MathMap. +* nfwd +* Number of forward transformation functions for the MathMap. +* ninv +* Number of inverse transformation functions for the MathMap. +* forward +* A non-zero value indicates the the MathMap's forward transformation +* functions are being processed, while a zero value indicates processing +* of the inverse transformation functions. This value, together with +* "nin", "nout", "nfwd" and "ninv" are used solely for constructing +* error messages. +* var +* Address in which to return a pointer to an array (with "nfun" +* elements) of pointers to null terminated strings containing the +* extracted variable names (i.e. this returns an array of strings). +* +* Both the returned array of pointers, and the strings to which they +* point, will be stored in dynamically allocated memory and should +* be freed by the caller (using astFree) when no longer required. +* status +* Pointer to the inherited status variable. + +* Notes: +* - A NULL value will be returned for "*var" if this function is +* invoked with the global error status set, or if it should fail for +* any reason. +*/ + +/* Local Variables: */ + char *duser1; /* Transformation direction for function */ + char *duser2; /* Transformation direction for function */ + char c; /* Extracted character */ + int i1; /* Loop counter for detecting duplicates */ + int i2; /* Loop counter for detecting duplicates */ + int i; /* Loop counter for characters */ + int iend; /* Last character index in parsed name */ + int ifun; /* Loop counter for functions */ + int iuser1; /* Function number as known to the user */ + int iuser2; /* Function number as known to the user */ + int nc; /* Character count */ + int nextra; /* Number of intermediate functions */ + int nprimary; /* Number of primary input/output variables */ + +/* Initialise. */ + *var = NULL; + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Obtain the number of primary input/output variables, depending on + the direction of the coordinate transformation. */ + nprimary = ( forward ? nin : nout ); + +/* Deterine the number of extra (intermediate) functions that come + before these primary ones. These affect the numbering of + transformation functions as known to the user, and must be accounted + for when reporting error messages. */ + nextra = ( forward ? ninv - nin : nfwd - nout ); + +/* Allocate and initialise memory for the returned array of pointers. */ + MALLOC_POINTER_ARRAY( *var, char *, nfun ) + +/* Loop to process each function in turn. */ + if ( astOK ) { + for ( ifun = 0; ifun < nfun; ifun++ ) { + +/* Count the number of characters appearing before the "=" sign (or in + the entire string if the "=" is absent). */ + for ( nc = 0; ( c = fun[ ifun ][ nc ] ); nc++ ) if ( c == '=' ) break; + +/* If no characters were counted, then report an appropriate error + message, depending on whether the function string was entirely + blank. */ + if ( !nc ) { + if ( c ) { + astError( AST__MISVN, + "%s(%s): No left hand side in expression: \"%s\".", status, + method, class, fun[ ifun ] ); + } else { + astError( AST__MISVN, + "%s: Transformation function contains no variable " + "name.", status, + method ); + } + break; + } + +/* If OK, allocate memory to hold the output string and check for + errors. */ + ( *var )[ ifun ] = astMalloc( sizeof( char ) * (size_t) ( nc + 1 ) ) ; + if ( !astOK ) break; + +/* If OK, copy the characters before the "=" sign to the new + string. */ + nc = 0; + for ( i = 0; ( c = fun[ ifun ][ i ] ); i++ ) { + if ( c == '=' ) break; + ( *var )[ ifun ][ nc++] = c; + } + +/* Null terminate the result. */ + ( *var )[ ifun ][ nc ] = '\0'; + +/* Try to parse the contents of the extracted string as a name. */ + ParseName( ( *var )[ ifun ], 0, &iend, status ); + +/* If unsuccessful, or if all the characters were not parsed, then we + have an invalid variable name, so report an error and quit. */ + if ( ( iend < 0 ) || ( *var )[ ifun ][ iend + 1 ] ) { + astError( AST__VARIN, + "%s(%s): Variable name is invalid: \"%s\".", status, + method, class, ( *var )[ ifun ] ); + break; + } + } + +/* If an error occurred above, then determine the function number, and + the direction of the transformation of which it forms part, as known + to the user. */ + if ( !astOK ) { + if ( ifun < nprimary ) { + iuser1 = ifun + 1 + nextra; + duser1 = ( forward ? "inverse" : "forward" ); + } else { + iuser1 = ifun + 1 - nprimary; + duser1 = ( forward ? "forward" : "inverse" ); + } + +/* Report a contextual error message. */ + astError( astStatus, + "Error in %s transformation function %d.", status, + duser1, iuser1 ); + } + } + +/* If there has been no error, loop to compare all the variable names + with each other to detect duplication. */ + if ( astOK ) { + for ( i1 = 1; i1 < nfun; i1++ ) { + for ( i2 = 0; i2 < i1; i2++ ) { + +/* If a duplicate variable name is found, report an error. */ + if ( !strcmp( ( *var )[ i1 ], ( *var )[ i2 ] ) ) { + astError( AST__DUVAR, + "%s(%s): Duplicate definition of variable name: " + "\"%s\".", status, + method, class, ( *var )[ i1 ] ); + +/* For each transformation function involved, determine the function + number and the direction of the transformation of which it forms part, + as known to the user. */ + if ( i1 < nprimary ) { + iuser1 = i1 + 1 + nextra; + duser1 = ( forward ? "inverse" : "forward" ); + } else { + iuser1 = i1 + 1 - nprimary; + duser1 = ( forward ? "forward" : "inverse" ); + } + if ( i2 < nprimary ) { + iuser2 = i2 + 1 + nextra; + duser2 = ( forward ? "inverse" : "forward" ); + } else { + iuser2 = i2 + 1 - nprimary; + duser2 = ( forward ? "forward" : "inverse" ); + } + +/* Report a contextual error message. */ + astError( astStatus, + "Conflict between %s function %d and %s function %d.", status, + duser1, iuser1, duser2, iuser2 ); + break; + } + } + if ( !astOK ) break; + } + } + +/* If an error occurred, free any allocated memory and reset the + output value. */ + if ( !astOK ) { + FREE_POINTER_ARRAY( *var, nfun ) + } +} + +static double Gauss( Rcontext *context, int *status ) { +/* +* Name: +* Gauss + +* Purpose: +* Produce a pseudo-random sample from a standard Gaussian distribution. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* double Gauss( Rcontext *context, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* On each invocation, this function returns a pseudo-random sample drawn +* from a standard Gaussian distribution with mean zero and standard +* deviation unity. The Box-Muller transformation method is used. + +* Parameters: +* context +* Pointer to an Rcontext structure which holds the random number +* generator's context between invocations. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* A sample from a standard Gaussian distribution. + +* Notes: +* - The sequence of numbers returned is determined by the "seed" +* value in the Rcontext structure supplied. +* - If the seed value is changed, the "active" flag must also be cleared +* so that this function can re-initiallise the Rcontext structure before +* generating the next pseudo-random number. The "active" flag should +* also be clear to force initialisation the first time an Rcontext +* structure is used. +* - This function does not perform error checking and does not generate +* errors. It will execute even if the global error status is set. +*/ + +/* Local Variables: */ + double rsq; /* Squared radius */ + double s; /* Scale factor */ + double x; /* First result value */ + static double y; /* Second result value */ + static int ysaved = 0; /* Previously-saved value available? */ + + LOCK_MUTEX7 + +/* If the random number generator context is not active, then it will + be (re)initialised on the first invocation of Rand (below). Ensure + that any previously-saved value within this function is first + discarded. */ + if ( !context->active ) ysaved = 0; + +/* If there is a previously-saved value available, then use it and + mark it as no longer available. */ + if ( ysaved ) { + x = y; + ysaved = 0; + +/* Otherwise, loop until a suitable new pair of values has been + obtained. */ + } else { + while ( 1 ) { + +/* Loop to obtain two random values uniformly distributed inside the + unit circle, while avoiding the origin (which maps to an infinite + result). */ + do { + x = 2.0 * Rand( context, status ) - 1.0; + y = 2.0 * Rand( context, status ) - 1.0; + rsq = x * x + y * y; + } while ( ( rsq >= 1.0 ) || ( rsq == 0.0 ) ); + +/* Perform the Box-Muller transformation, checking that this will not + produce overflow (which is extremely unlikely). If overflow would + occur, we simply repeat the above steps with a new pair of random + numbers. */ + s = -2.0 * log( rsq ); + if ( ( DBL_MAX * rsq ) >= s ) { + s = sqrt( s / rsq ); + +/* Scale the original random values to give a pair of results. One will be + returned and the second kept until next time. */ + x *= s; + y *= s; + break; + } + } + +/* Note that a saved value is available. */ + ysaved = 1; + } + + UNLOCK_MUTEX7 + +/* Return the current result. */ + return x; +} + +static int GetObjSize( AstObject *this_object, int *status ) { +/* +* Name: +* GetObjSize + +* Purpose: +* Return the in-memory size of an Object. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* int GetObjSize( AstObject *this, int *status ) + +* Class Membership: +* MathMap member function (over-rides the astGetObjSize protected +* method inherited from the parent class). + +* Description: +* This function returns the in-memory size of the supplied MathMap, +* in bytes. + +* Parameters: +* this +* Pointer to the MathMap. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* The Object size, in bytes. + +* Notes: +* - A value of zero will be returned if this function is invoked +* with the global status set, or if it should fail for any reason. +*/ + +/* Local Variables: */ + AstMathMap *this; /* Pointer to MathMap structure */ + int result; /* Result value to return */ + +/* Initialise. */ + result = 0; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Obtain a pointers to the MathMap structure. */ + this = (AstMathMap *) this_object; + +/* Invoke the GetObjSize method inherited from the parent class, and then + add on any components of the class structure defined by thsi class + which are stored in dynamically allocated memory. */ + result = (*parent_getobjsize)( this_object, status ); + + SIZEOF_POINTER_ARRAY( this->fwdfun, this->nfwd ) + SIZEOF_POINTER_ARRAY( this->invfun, this->ninv ) + SIZEOF_POINTER_ARRAY( this->fwdcode, this->nfwd ) + SIZEOF_POINTER_ARRAY( this->invcode, this->ninv ) + SIZEOF_POINTER_ARRAY( this->fwdcon, this->nfwd ) + SIZEOF_POINTER_ARRAY( this->invcon, this->ninv ) + +/* If an error occurred, clear the result value. */ + if ( !astOK ) result = 0; + +/* Return the result, */ + return result; +} + +static const char *GetAttrib( AstObject *this_object, const char *attrib, int *status ) { +/* +* Name: +* GetAttrib + +* Purpose: +* Get the value of a specified attribute for a MathMap. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* const char *GetAttrib( AstObject *this, const char *attrib, int *status ) + +* Class Membership: +* MathMap member function (over-rides the protected astGetAttrib +* method inherited from the Mapping class). + +* Description: +* This function returns a pointer to the value of a specified +* attribute for a MathMap, formatted as a character string. + +* Parameters: +* this +* Pointer to the MathMap. +* attrib +* Pointer to a null-terminated string containing the name of +* the attribute whose value is required. This name should be in +* lower case, with all white space removed. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* - Pointer to a null-terminated string containing the attribute +* value. + +* Notes: +* - The returned string pointer may point at memory allocated +* within the MathMap, or at static memory. The contents of the +* string may be over-written or the pointer may become invalid +* following a further invocation of the same function or any +* modification of the MathMap. A copy of the string should +* therefore be made if necessary. +* - A NULL pointer will be returned if this function is invoked +* with the global error status set, or if it should fail for any +* reason. +*/ + +/* Local Variables: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ + AstMathMap *this; /* Pointer to the MathMap structure */ + const char *result; /* Pointer value to return */ + int ival; /* Integer attribute value */ + +/* Initialise. */ + result = NULL; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(this_object); + +/* Obtain a pointer to the MathMap structure. */ + this = (AstMathMap *) this_object; + +/* Compare "attrib" with each recognised attribute name in turn, + obtaining the value of the required attribute. If necessary, write + the value into "getattrib_buff" as a null-terminated string in an appropriate + format. Set "result" to point at the result string. */ + +/* Seed. */ +/* ----- */ + if ( !strcmp( attrib, "seed" ) ) { + ival = astGetSeed( this ); + if ( astOK ) { + (void) sprintf( getattrib_buff, "%d", ival ); + result = getattrib_buff; + } + +/* SimpFI. */ +/* ------- */ + } else if ( !strcmp( attrib, "simpfi" ) ) { + ival = astGetSimpFI( this ); + if ( astOK ) { + (void) sprintf( getattrib_buff, "%d", ival ); + result = getattrib_buff; + } + +/* SimpIF. */ +/* ------- */ + } else if ( !strcmp( attrib, "simpif" ) ) { + ival = astGetSimpIF( this ); + if ( astOK ) { + (void) sprintf( getattrib_buff, "%d", ival ); + result = getattrib_buff; + } + +/* If the attribute name was not recognised, pass it on to the parent + method for further interpretation. */ + } else { + result = (*parent_getattrib)( this_object, attrib, status ); + } + +/* Return the result. */ + return result; + +} + +void astInitMathMapVtab_( AstMathMapVtab *vtab, const char *name, int *status ) { +/* +*+ +* Name: +* astInitMathMapVtab + +* Purpose: +* Initialise a virtual function table for a MathMap. + +* Type: +* Protected function. + +* Synopsis: +* #include "mathmap.h" +* void astInitMathMapVtab( AstMathMapVtab *vtab, const char *name ) + +* Class Membership: +* MathMap vtab initialiser. + +* Description: +* This function initialises the component of a virtual function +* table which is used by the MathMap class. + +* Parameters: +* vtab +* Pointer to the virtual function table. The components used by +* all ancestral classes will be initialised if they have not already +* been initialised. +* name +* Pointer to a constant null-terminated character string which contains +* the name of the class to which the virtual function table belongs (it +* is this pointer value that will subsequently be returned by the Object +* astClass function). +*- +*/ + +/* Local Variables: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ + AstMappingVtab *mapping; /* Pointer to Mapping component of Vtab */ + AstObjectVtab *object; /* Pointer to Object component of Vtab */ + +/* Check the local error status. */ + if ( !astOK ) return; + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(NULL); + +/* Initialize the component of the virtual function table used by the + parent class. */ + astInitMappingVtab( (AstMappingVtab *) vtab, name ); + +/* Store a unique "magic" value in the virtual function table. This + will be used (by astIsAMathMap) to determine if an object belongs + to this class. We can conveniently use the address of the (static) + class_check variable to generate this unique value. */ + vtab->id.check = &class_check; + vtab->id.parent = &(((AstMappingVtab *) vtab)->id); + +/* Initialise member function pointers. */ +/* ------------------------------------ */ +/* Store pointers to the member functions (implemented here) that + provide virtual methods for this class. */ + vtab->ClearSeed = ClearSeed; + vtab->ClearSimpFI = ClearSimpFI; + vtab->ClearSimpIF = ClearSimpIF; + vtab->GetSeed = GetSeed; + vtab->GetSimpFI = GetSimpFI; + vtab->GetSimpIF = GetSimpIF; + vtab->SetSeed = SetSeed; + vtab->SetSimpFI = SetSimpFI; + vtab->SetSimpIF = SetSimpIF; + vtab->TestSeed = TestSeed; + vtab->TestSimpFI = TestSimpFI; + vtab->TestSimpIF = TestSimpIF; + +/* Save the inherited pointers to methods that will be extended, and + replace them with pointers to the new member functions. */ + object = (AstObjectVtab *) vtab; + mapping = (AstMappingVtab *) vtab; + parent_getobjsize = object->GetObjSize; + object->GetObjSize = GetObjSize; + + parent_clearattrib = object->ClearAttrib; + object->ClearAttrib = ClearAttrib; + parent_getattrib = object->GetAttrib; + object->GetAttrib = GetAttrib; + parent_setattrib = object->SetAttrib; + object->SetAttrib = SetAttrib; + parent_testattrib = object->TestAttrib; + object->TestAttrib = TestAttrib; + + parent_transform = mapping->Transform; + mapping->Transform = Transform; + +/* Store replacement pointers for methods which will be over-ridden by + new member functions implemented here. */ + object->Equal = Equal; + mapping->MapMerge = MapMerge; + +/* Declare the copy constructor, destructor and class dump function. */ + astSetCopy( vtab, Copy ); + astSetDelete( vtab, Delete ); + astSetDump( vtab, Dump, "MathMap", + "Transformation using mathematical functions" ); + +/* If we have just initialised the vtab for the current class, indicate + that the vtab is now initialised, and store a pointer to the class + identifier in the base "object" level of the vtab. */ + if( vtab == &class_vtab ) { + class_init = 1; + astSetVtabClassIdentifier( vtab, &(vtab->id) ); + } +} + +static double LogGamma( double x, int *status ) { +/* +* Name: +* LogGamma + +* Purpose: +* Calculate the logarithm of the gamma function. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* double LogGamma( double x, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* This function returns the natural logarithm of the gamma function +* for real arguments x>0. It uses the approximation of Lanczos, with +* constants from Press et al. (Numerical Recipes), giving a maximum +* fractional error (on the gamma function) of less than 2e-10. + +* Parameters: +* x +* The function argument, which must be greater than zero. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* The natural logarithm of the gamma function with "x" as argument, +* or AST__BAD if "x" is not greater than zero. + +* Notes: +* - This function does not generate errors and does not perform error +* reporting. It will execute even if the global error status is set. +*/ + +/* Local Constants: */ + const double c0 = 1.000000000190015; /* Coefficients for series sum... */ + const double c1 = 76.18009172947146; + const double c2 = -86.50532032941677; + const double c3 = 24.01409824083091; + const double c4 = -1.231739572450155; + const double c5 = 0.1208650973866179e-2; + const double c6 = -0.5395239384953e-5; + const double g = 5.0; + +/* Local Variables: */ + double result; /* Result value to return */ + double sum; /* Series sum */ + double xx; /* Denominator for summing series */ + static double root_twopi; /* sqrt( 2.0 * pi ) */ + static int init = 0; /* Initialisation performed? */ + +/* If initialisation has not yet been performed, calculate the + constant required below. */ + LOCK_MUTEX3 + if ( !init ) { + root_twopi = sqrt( 2.0 * acos( -1.0 ) ); + +/* Note that initialisation has been performed. */ + init = 1; + } + UNLOCK_MUTEX3 + +/* Return a bad value if "x" is not greater than zero. */ + if ( x <= 0.0 ) { + result = AST__BAD; + +/* Otherwise, form the series sum. Since we only use 6 terms, the loop + that would normally be used has been completely unrolled here. */ + } else { + xx = x; + sum = c0; + sum += c1 / ++xx; + sum += c2 / ++xx; + sum += c3 / ++xx; + sum += c4 / ++xx; + sum += c5 / ++xx; + sum += c6 / ++xx; + +/* Calculate the result. */ + result = x + g + 0.5; + result -= ( x + 0.5 ) * log( result ); + result = log( root_twopi * sum / x ) - result; + } + +/* Return the result. */ + return result; +} + +static int MapMerge( AstMapping *this, int where, int series, int *nmap, + AstMapping ***map_list, int **invert_list, int *status ) { +/* +* Name: +* MapMerge + +* Purpose: +* Simplify a sequence of Mappings containing a MathMap. + +* Type: +* Private function. + +* Synopsis: +* #include "mapping.h" +* int MapMerge( AstMapping *this, int where, int series, int *nmap, +* AstMapping ***map_list, int **invert_list, int *status ) + +* Class Membership: +* MathMap method (over-rides the protected astMapMerge method +* inherited from the Mapping class). + +* Description: +* This function attempts to simplify a sequence of Mappings by +* merging a nominated MathMap in the sequence with its neighbours, +* so as to shorten the sequence if possible. +* +* In many cases, simplification will not be possible and the +* function will return -1 to indicate this, without further +* action. +* +* In most cases of interest, however, this function will either +* attempt to replace the nominated MathMap with one which it +* considers simpler, or to merge it with the Mappings which +* immediately precede it or follow it in the sequence (both will +* normally be considered). This is sufficient to ensure the +* eventual simplification of most Mapping sequences by repeated +* application of this function. +* +* In some cases, the function may attempt more elaborate +* simplification, involving any number of other Mappings in the +* sequence. It is not restricted in the type or scope of +* simplification it may perform, but will normally only attempt +* elaborate simplification in cases where a more straightforward +* approach is not adequate. + +* Parameters: +* this +* Pointer to the nominated MathMap which is to be merged with +* its neighbours. This should be a cloned copy of the MathMap +* pointer contained in the array element "(*map_list)[where]" +* (see below). This pointer will not be annulled, and the +* MathMap it identifies will not be modified by this function. +* where +* Index in the "*map_list" array (below) at which the pointer +* to the nominated MathMap resides. +* series +* A non-zero value indicates that the sequence of Mappings to +* be simplified will be applied in series (i.e. one after the +* other), whereas a zero value indicates that they will be +* applied in parallel (i.e. on successive sub-sets of the +* input/output coordinates). +* nmap +* Address of an int which counts the number of Mappings in the +* sequence. On entry this should be set to the initial number +* of Mappings. On exit it will be updated to record the number +* of Mappings remaining after simplification. +* map_list +* Address of a pointer to a dynamically allocated array of +* Mapping pointers (produced, for example, by the astMapList +* method) which identifies the sequence of Mappings. On entry, +* the initial sequence of Mappings to be simplified should be +* supplied. +* +* On exit, the contents of this array will be modified to +* reflect any simplification carried out. Any form of +* simplification may be performed. This may involve any of: (a) +* removing Mappings by annulling any of the pointers supplied, +* (b) replacing them with pointers to new Mappings, (c) +* inserting additional Mappings and (d) changing their order. +* +* The intention is to reduce the number of Mappings in the +* sequence, if possible, and any reduction will be reflected in +* the value of "*nmap" returned. However, simplifications which +* do not reduce the length of the sequence (but improve its +* execution time, for example) may also be performed, and the +* sequence might conceivably increase in length (but normally +* only in order to split up a Mapping into pieces that can be +* more easily merged with their neighbours on subsequent +* invocations of this function). +* +* If Mappings are removed from the sequence, any gaps that +* remain will be closed up, by moving subsequent Mapping +* pointers along in the array, so that vacated elements occur +* at the end. If the sequence increases in length, the array +* will be extended (and its pointer updated) if necessary to +* accommodate any new elements. +* +* Note that any (or all) of the Mapping pointers supplied in +* this array may be annulled by this function, but the Mappings +* to which they refer are not modified in any way (although +* they may, of course, be deleted if the annulled pointer is +* the final one). +* invert_list +* Address of a pointer to a dynamically allocated array which, +* on entry, should contain values to be assigned to the Invert +* attributes of the Mappings identified in the "*map_list" +* array before they are applied (this array might have been +* produced, for example, by the astMapList method). These +* values will be used by this function instead of the actual +* Invert attributes of the Mappings supplied, which are +* ignored. +* +* On exit, the contents of this array will be updated to +* correspond with the possibly modified contents of the +* "*map_list" array. If the Mapping sequence increases in +* length, the "*invert_list" array will be extended (and its +* pointer updated) if necessary to accommodate any new +* elements. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* If simplification was possible, the function returns the index +* in the "map_list" array of the first element which was +* modified. Otherwise, it returns -1 (and makes no changes to the +* arrays supplied). + +* Notes: +* - A value of -1 will be returned if this function is invoked +* with the global error status set, or if it should fail for any +* reason. +*/ + +/* Local Variables: */ + AstMapping *new; /* Pointer to replacement Mapping */ + AstMathMap *mathmap1; /* Pointer to first MathMap */ + AstMathMap *mathmap2; /* Pointer to second MathMap */ + char **fwd1; /* Pointer to first forward function array */ + char **fwd2; /* Pointer to second forward function array */ + char **inv1; /* Pointer to first inverse function array */ + char **inv2; /* Pointer to second inverse function array */ + int ifun; /* Loop counter for functions */ + int imap1; /* Index of first Mapping */ + int imap2; /* Index of second Mapping */ + int imap; /* Loop counter for Mappings */ + int invert1; /* Invert flag for first MathMap */ + int invert2; /* Invert flag for second MathMap */ + int nfwd1; /* No. forward functions for first MathMap */ + int nfwd2; /* No. forward functions for second MathMap */ + int nin1; /* Number input coords for first MathMap */ + int ninv1; /* No. inverse functions for first MathMap */ + int ninv2; /* No. inverse functions for second MathMap */ + int nout2; /* Number output coords for second MathMap */ + int result; /* Result value to return */ + int simplify; /* Mappings may simplify? */ + +/* Initialise the returned result. */ + result = -1; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Initialise variables to avoid "used of uninitialised variable" + messages from dumb compilers. */ + mathmap1 = NULL; + mathmap2 = NULL; + imap1 = 0; + imap2 = 0; + invert1 = 0; + invert2 = 0; + nfwd1 = 0; + nin1 = 0; + ninv1 = 0; + +/* MathMaps are only worth simplifying if they occur in series. */ + simplify = series; + +/* If simplification appears possible, then obtain the indices of the + nominated mapping and of the one which follows it. Check that a + mapping exists for the second index. */ + if ( simplify ) { + imap1 = where; + imap2 = imap1 + 1; + simplify = ( imap2 < *nmap ); + } + +/* If OK, check whether the class of both Mappings is "MathMap" (a + MathMap can only combine with another MathMap). */ + if ( simplify ) { + simplify = !strcmp( astGetClass( ( *map_list )[ imap1 ] ), "MathMap" ); + } + if ( astOK && simplify ) { + simplify = !strcmp( astGetClass( ( *map_list )[ imap2 ] ), "MathMap" ); + } + +/* If still OK, obtain pointers to the two MathMaps and the associated + invert flag values. */ + if ( astOK && simplify ) { + mathmap1 = (AstMathMap *) ( *map_list )[ imap1 ]; + mathmap2 = (AstMathMap *) ( *map_list )[ imap2 ]; + invert1 = ( *invert_list )[ imap1 ]; + invert2 = ( *invert_list )[ imap2 ]; + +/* Depending on the invert flag values, obtain the SimpFI or SimpIF + attribute value from each MathMap and check whether they are set so as + to permit simplification. */ + simplify = ( ( invert1 ? astGetSimpIF( mathmap1 ) : + astGetSimpFI( mathmap1 ) ) && + ( invert2 ? astGetSimpFI( mathmap2 ) : + astGetSimpIF( mathmap2 ) ) ); + } + +/* If still OK, obtain the effective numbers of input coordinates for + the first MathMap and output coordinates for the second. Take account + of the associated invert flags and the way the Invert attribute of + each MathMap is currently set. */ + if ( astOK && simplify ) { + nin1 = ( invert1 == astGetInvert( mathmap1 ) ) ? + astGetNin( mathmap1 ) : astGetNout( mathmap1 ); + nout2 = ( invert2 == astGetInvert( mathmap2 ) ) ? + astGetNout( mathmap2 ) : astGetNin( mathmap2 ); + +/* Simplification is only possible if these two numbers are equal + (otherwise the the two MathMaps cannot be identical). */ + simplify = ( nin1 == nout2 ); + } + +/* If still OK, obtain the effective number of forward transformation + functions for the first MathMap (allowing for the associated invert + flag). Similarly, obtain the effective number of inverse + transformation functions for the second MathMap. */ + if ( astOK && simplify ) { + nfwd1 = !invert1 ? mathmap1->nfwd : mathmap1->ninv; + ninv2 = !invert2 ? mathmap2->ninv : mathmap2->nfwd; + +/* Check whether these values are equal. The MathMaps cannot be + identical if they are not. */ + simplify = ( nfwd1 == ninv2 ); + } + +/* As above, obtain pointers to the array of effective forward + transformation functions for the first MathMap, and the effective + inverse transformation functions for the second MathMap. */ + if ( astOK && simplify ) { + fwd1 = !invert1 ? mathmap1->fwdfun : mathmap1->invfun; + inv2 = !invert2 ? mathmap2->invfun : mathmap2->fwdfun; + +/* Loop to check whether these two sets of functions are + identical. The MathMaps cannot be merged unless they are. */ + for ( ifun = 0; ifun < nfwd1; ifun++ ) { + simplify = !strcmp( fwd1[ ifun ], inv2[ ifun ] ); + if ( !simplify ) break; + } + } + +/* If OK, repeat the above process to compare the effective inverse + transformation functions of the first MathMap with the forward + functions of the second one. */ + if ( astOK && simplify ) { + ninv1 = !invert1 ? mathmap1->ninv : mathmap1->nfwd; + nfwd2 = !invert2 ? mathmap2->nfwd : mathmap2->ninv; + simplify = ( ninv1 == nfwd2 ); + } + if ( astOK && simplify ) { + inv1 = !invert1 ? mathmap1->invfun : mathmap1->fwdfun; + fwd2 = !invert2 ? mathmap2->fwdfun : mathmap2->invfun; + for ( ifun = 0; ifun < ninv1; ifun++ ) { + simplify = !strcmp( inv1[ ifun ], fwd2[ ifun ] ); + if ( !simplify ) break; + } + } + +/* If the two MathMaps can be merged, create a UnitMap as a + replacement. */ + if ( astOK && simplify ) { + new = (AstMapping *) astUnitMap( nin1, "", status ); + +/* If OK, annul the pointers to the original MathMaps. */ + if ( astOK ) { + ( *map_list )[ imap1 ] = astAnnul( ( *map_list )[ imap1 ] ); + ( *map_list )[ imap2 ] = astAnnul( ( *map_list )[ imap2 ] ); + +/* Insert the pointer to the replacement UnitMap and store the + associated invert flag. */ + ( *map_list )[ imap1 ] = new; + ( *invert_list )[ imap1 ] = 0; + +/* Loop to move the following Mapping pointers and invert flags down + in their arrays to close the gap. */ + for ( imap = imap2 + 1; imap < *nmap; imap++ ) { + ( *map_list )[ imap - 1 ] = ( *map_list )[ imap ]; + ( *invert_list )[ imap - 1 ] = ( *invert_list )[ imap ]; + } + +/* Clear the final entry in each array. */ + ( *map_list )[ *nmap - 1 ] = NULL; + ( *invert_list )[ *nmap - 1 ] = 0; + +/* Decrement the Mapping count and return the index of the first + modified element. */ + ( *nmap )--; + result = imap1; + } + } + +/* If an error occurred, clear the returned value. */ + if ( !astOK ) result = -1; + +/* Return the result. */ + return result; +} + +static void ParseConstant( const char *method, const char *class, + const char *exprs, int istart, int *iend, + double *con, int *status ) { +/* +* Name: +* ParseConstant + +* Purpose: +* Parse a constant. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* void ParseConstant( const char *method, const char *class, +* const char *exprs, int istart, int *iend, +* double *con, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* This routine parses an expression, looking for a constant starting at +* the character with index "istart" in the string "exprs". If it +* identifies the constant successfully, "*con" it will return its value +* and "*iend" will be set to the index of the final constant character +* in "exprs". +* +* If the characters encountered are clearly not part of a constant (it +* does not begin with a numeral or decimal point) the function returns +* with "*con" set to zero and "*iend" set to -1, but without reporting +* an error. However, if the first character appears to be a constant but +* its syntax proves to be invalid, then an error is reported. +* +* The expression must be in lower case with no embedded white space. +* The constant must not have a sign (+ or -) in front of it. + +* Parameters: +* method +* Pointer to a constant null-terminated character string +* containing the name of the method that invoked this function. +* This method name is used solely for constructing error messages. +* class +* Pointer to a constant null-terminated character string containing the +* class name of the Object being processed. This name is used solely +* for constructing error messages. +* exprs +* Pointer to a null-terminated string containing the expression +* to be parsed. +* istart +* Index of the first character in "exprs" to be considered by this +* function. +* iend +* Pointer to an int in which to return the index in "exprs" of the +* final character which forms part of the constant. If no constant is +* found, a value of -1 is returned. +* con +* Pointer to a double, in which the value of the constant, if found, +* will be returned. +* status +* Pointer to the inherited status variable. +*/ + +/* Local Variables: */ + char *str; /* Pointer to temporary string */ + char c; /* Single character from the expression */ + int dpoint; /* Decimal point encountered? */ + int expon; /* Exponent character encountered? */ + int i; /* Loop counter for characters */ + int iscon; /* Character is part of the constant? */ + int n; /* Number of values read by astSscanf */ + int nc; /* Number of characters read by astSscanf */ + int numer; /* Numeral encountered in current field? */ + int sign; /* Sign encountered? */ + int valid; /* Constant syntax valid? */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Initialise. */ + *con = 0.0; + *iend = -1; + +/* Check if the expression starts with a numeral or a decimal point. */ + c = exprs[ istart ]; + numer = isdigit( c ); + dpoint = ( c == '.' ); + +/* If it begins with any of these, the expression is clearly intended + to be a constant, so any failure beyond this point will result in an + error. Otherwise, failure to find a constant is not an error. */ + if ( numer || dpoint ) { + +/* Initialise remaining variables specifying the parser context. */ + expon = 0; + sign = 0; + valid = 1; + +/* Loop to increment the last constant character position until the + following character in the expression does not look like part of the + constant. */ + *iend = istart; + iscon = 1; + while ( ( c = exprs[ *iend + 1 ] ) && iscon ) { + iscon = 0; + +/* It may be part of a numerical constant if it is a numeral, wherever + it occurs. */ + if ( isdigit( c ) ) { + numer = 1; + iscon = 1; + +/* Or a decimal point, so long as it is the first one and is not in + the exponent field. Otherwise it is invalid. */ + } else if ( c == '.' ) { + if ( !( dpoint || expon ) ) { + dpoint = 1; + iscon = 1; + } else { + valid = 0; + } + +/* Or if it is a 'd' or 'e' exponent character, so long as it is the + first one and at least one numeral has been encountered first. + Otherwise it is invalid. */ + } else if ( ( c == 'd' ) || ( c == 'e' ) ) { + if ( !expon && numer ) { + expon = 1; + numer = 0; + iscon = 1; + } else { + valid = 0; + } + +/* Or if it is a sign, so long as it is in the exponent field and is + the first sign with no previous numerals in the same field. Otherwise + it is invalid (unless numerals have been encountered, in which case it + marks the end of the constant). */ + } else if ( ( c == '+' ) || ( c == '-' ) ) { + if ( expon && !sign && !numer ) { + sign = 1; + iscon = 1; + } else if ( !numer ) { + valid = 0; + } + } + +/* Increment the character count if the next character may be part of + the constant, or if it was invalid (it will then form part of the + error message). */ + if ( iscon || !valid ) ( *iend )++; + } + +/* Finally, check that the last field contained a numeral. */ + valid = ( valid && numer ); + +/* If the constant appears valid, allocate a temporary string to hold + it. */ + if ( valid ) { + str = astMalloc( (size_t) ( *iend - istart + 2 ) ); + if ( astOK ) { + +/* Copy the constant's characters, changing 'd' to 'e' so that + "astSscanf" will recognise it as an exponent character. */ + for ( i = istart; i <= *iend; i++ ) { + str[ i - istart ] = ( exprs[ i ] == 'd' ) ? 'e' : exprs[ i ]; + } + str[ *iend - istart + 1 ] = '\0'; + +/* Attempt to read the constant as a double, noting how many values + are read and how many characters consumed. */ + n = astSscanf( str, "%lf%n", con, &nc ); + +/* Check that one value was read and all the characters consumed. If + not, then the constant's syntax is invalid. */ + if ( ( n != 1 ) || ( nc < ( *iend - istart + 1 ) ) ) valid = 0; + } + +/* Free the temporary string. */ + str = astFree( str ); + } + +/* If the constant syntax is invalid, and no other error has occurred, + then report an error. */ + if ( astOK && !valid ) { + astError( AST__CONIN, + "%s(%s): Invalid constant syntax in the expression " + "\"%.*s\".", status, + method, class, *iend + 1, exprs ); + } + +/* If an error occurred, reset the output values. */ + if ( !astOK ) { + *iend = -1; + *con = 0.0; + } + } +} + +static void ParseName( const char *exprs, int istart, int *iend, int *status ) { +/* +* Name: +* ParseName + +* Purpose: +* Parse a name. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* void ParseName( const char *exprs, int istart, int *iend, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* This routine parses an expression, looking for a name starting at the +* character with index "istart" in the string "exprs". If it identifies +* a name successfully, "*iend" will return the index of the final name +* character in "exprs". A name must begin with an alphabetic character +* and subsequently contain only alphanumeric characters or underscores. +* +* If the expression does not contain a name at the specified location, +* "*iend" is set to -1. No error results. +* +* The expression should not contain embedded white space. + +* Parameters: +* exprs +* Pointer to a null-terminated string containing the expression +* to be parsed. +* istart +* Index of the first character in "exprs" to be considered by this +* function. +* iend +* Pointer to an int in which to return the index in "exprs" of the +* final character which forms part of the name. If no name is +* found, a value of -1 is returned. +* status +* Pointer to the inherited status variable. +*/ + +/* Local Variables: */ + char c; /* Single character from expression */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Initialise. */ + *iend = -1; + +/* Check the first character is valid for a name (alphabetic). */ + if ( isalpha( exprs[ istart ] ) ) { + +/* If so, loop to inspect each subsequent character until one is found + which is not part of a name (not alphanumeric or underscore). */ + for ( *iend = istart; ( c = exprs[ *iend + 1 ] ); ( *iend )++ ) { + if ( !( isalnum( c ) || ( c == '_' ) ) ) break; + } + } +} + +static void ParseVariable( const char *method, const char *class, + const char *exprs, int istart, int nvar, + const char *var[], int *ivar, int *iend, int *status ) { +/* +* Name: +* ParseVariable + +* Purpose: +* Parse a variable name. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* void ParseVariable( const char *method, const char *class, +* const char *exprs, int istart, int nvar, +* const char *var[], int *ivar, int *iend, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* This routine parses an expression, looking for a recognised variable +* name starting at the character with index "istart" in the string +* "exprs". If it identifies a variable name successfully, "*ivar" will +* return a value identifying it and "*iend" will return the index of the +* final variable name character in "exprs". To be recognised, a name +* must begin with an alphabetic character and subsequently contain only +* alphanumeric characters or underscores. It must also appear in the +* list of defined variable names supplied to this function. +* +* If the expression does not contain a name at the specified location, +* "*ivar" and "*iend" are set to -1 and no error results. However, if +* the expression contains a name but it is not in the list of defined +* variable names supplied, then an error is reported. +* +* This function is case sensitive. The expression should not contain +* embedded white space. + +* Parameters: +* method +* Pointer to a constant null-terminated character string +* containing the name of the method that invoked this function. +* This method name is used solely for constructing error messages. +* class +* Pointer to a constant null-terminated character string containing the +* class name of the Object being processed. This name is used solely +* for constructing error messages. +* exprs +* Pointer to a null-terminated string containing the expression +* to be parsed. +* istart +* Index of the first character in "exprs" to be considered by this +* function. +* nvar +* The number of defined variable names. +* var +* An array of pointers (with "nvar" elements) to null-terminated +* strings. Each of these should contain a variable name to be +* recognised. These strings are case sensitive and should contain +* no white space. +* ivar +* Pointer to an int in which to return the index in "vars" of the +* variable name found. If no variable name is found, a value of -1 +* is returned. +* iend +* Pointer to an int in which to return the index in "exprs" of the +* final character which forms part of the variable name. If no variable +* name is found, a value of -1 is returned. +* status +* Pointer to the inherited status variable. +*/ + +/* Local Variables: */ + int found; /* Variable name recognised? */ + int nc; /* Number of characters in variable name */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Initialise. */ + *ivar = -1; + *iend = -1; + +/* Determine if the characters in the expression starting at index + "istart" constitute a valid name. */ + ParseName( exprs, istart, iend, status ); + +/* If so, calculate the length of the name. */ + if ( *iend >= istart ) { + nc = *iend - istart + 1; + +/* Loop to compare the name with the list of variable names + supplied. */ + found = 0; + for ( *ivar = 0; *ivar < nvar; ( *ivar )++ ) { + found = ( nc == (int) strlen( var[ *ivar ] ) ) && + !strncmp( exprs + istart, var[ *ivar ], (size_t) nc ); + +/* Break if the name is recognised. */ + if ( found ) break; + } + +/* If it was not recognised, then report an error and reset the output + values. */ + if ( !found ) { + astError( AST__UDVOF, + "%s(%s): Undefined variable or function in the expression " + "\"%.*s\".", status, + method, class, *iend + 1, exprs ); + *ivar = -1; + *iend = -1; + } + } +} + +static double Poisson( Rcontext *context, double mean, int *status ) { +/* +* Name: +* Poisson + +* Purpose: +* Produce a pseudo-random sample from a Poisson distribution. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* double Poisson( Rcontext *context, double mean, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* On each invocation, this function returns a pseudo-random sample drawn +* from a Poisson distribution with a specified mean. A combination of +* methods is used, depending on the value of the mean. The algorithm is +* based on that given by Press et al. (Numerical Recipes), but +* re-implemented and extended. + +* Parameters: +* context +* Pointer to an Rcontext structure which holds the random number +* generator's context between invocations. +* mean +* The mean of the Poisson distribution, which should not be +* negative. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* A sample (which will only take integer values) from the Poisson +* distribution, or AST__BAD if the mean supplied is negative. + +* Notes: +* - The sequence of numbers returned is determined by the "seed" +* value in the Rcontext structure supplied. +* - If the seed value is changed, the "active" flag must also be cleared +* so that this function can re-initiallise the Rcontext structure before +* generating the next pseudo-random number. The "active" flag should +* also be clear to force initialisation the first time an Rcontext +* structure is used. +* - This function does not perform error checking and does not generate +* errors. It will execute even if the global error status is set. +*/ + +/* Local Constants: */ + const double small = 9.3; /* "Small" distribution mean value */ + +/* Local Variables: */ + double pfract; /* Probability of accepting sample */ + double product; /* Product of random samples */ + double ran; /* Sample from Lorentzian distribution */ + double result; /* Result value to return */ + static double beta; /* Constant for forming acceptance ratio */ + static double huge; /* Large mean where std. dev. is negligible */ + static double last_mean; /* Value of "mean" on last invocation */ + static double log_mean; /* Logarithm of "mean" */ + static double pi; /* Value of pi */ + static double ranmax; /* Maximum safe value of "ran" */ + static double root_2mean; /* sqrt( 2.0 * mean ) */ + static double sqrt_point9; /* Square root of 0.9 */ + static double thresh; /* Threshold for product of samples */ + static int init = 0; /* Local initialisation performed? */ + + LOCK_MUTEX6 + +/* If initialisation has not yet been performed, then perform it + now. */ + if ( !init ) { + +/* Initialise the mean value from the previous invocation. */ + last_mean = -1.0; + +/* Calculate simple constants. */ + pi = acos( -1.0 ); + sqrt_point9 = sqrt( 0.9 ); + +/* Calculate the value of the distribution mean for which the smallest + representable deviation from the mean permitted by the machine + precision is one thousand standard deviations. */ + huge = pow( 1.0e3 / DBL_EPSILON, 2.0 ); + +/* Calculate the largest value such that + (0.9+(sqrt_point9*ranmax)*(sqrt_point9*ranmax)) doesn't overflow, + allowing a small margin for rounding error. */ + ranmax = ( sqrt( DBL_MAX - 0.9 ) / sqrt( 0.9 ) ) * + ( 1.0 - 4.0 * DBL_EPSILON ); + +/* Note that initialisation has been performed. */ + init = 1; + } + +/* If the distribution mean is less than zero, then return a bad + result. */ + if ( mean < 0.0 ) { + result = AST__BAD; + +/* If the mean is zero, then the result can only be zero. */ + } else if ( mean == 0.0 ) { + result = 0.0; + +/* Otherwise, if the mean is sufficiently small, we can use the direct + method of summing a series of exponentially distributed random samples + and counting the number which occur before the mean is exceeded. This + is equivalent to multiplying a series of uniformly distributed + samples and counting the number which occur before the product + becomes less then an equivalent threshold. */ + } else if ( mean <= small ) { + +/* If the mean has changed since the last invocation, store the new + mean and calculate a new threshold. */ + if ( mean != last_mean ) { + last_mean = mean; + thresh = exp( -mean ); + } + +/* Initialise the product and the result. */ + product = 1.0; + result = -1.0; + +/* Multiply the random samples, counting the number needed to reach + the threshold. */ + do { + product *= Rand( context, status ); + result += 1.0; + } while ( product > thresh ); + +/* Otherwise, if the distribution mean is large (but not huge), we + must use an indirect rejection method. */ + } else if ( mean <= huge ) { + +/* If the mean has changed since the last invocation, then + re-calculate the constants required below. Note that because of the + restrictions we have placed on "mean", these calculations are safe + against overflow. */ + if ( mean != last_mean ) { + last_mean = mean; + log_mean = log( mean ); + root_2mean = sqrt( 2.0 * mean ); + beta = mean * log_mean - LogGamma( mean + 1.0, status ); + } + +/* Loop until a suitable random sample has been generated. */ + do { + do { + +/* First transform a sample from a uniform distribution to obtain a + sample from a Lorentzian distribution. Check that the result is not so + large as to cause overflow later. Also check for overflow in the maths + library. If necessary, obtain a new sample. */ + do { + errno = 0; + ran = tan( pi * Rand( context, status ) ); + } while ( ( ran > ranmax ) || + ( ( errno == ERANGE ) && + ( ( ( ran >= 0.0 ) ? ran : -ran ) == HUGE_VAL ) ) ); + +/* If OK, scale the sample and add a constant so that the sample's + distribution approximates the Poisson distribution we + require. Overflow is prevented by the check on "ran" above, together + with the restricted value of "mean". */ + result = ran * root_2mean + mean; + +/* If the result is less than zero (where the Poisson distribution has + value zero), then obtain a new sample. */ + } while ( result < 0.0 ); + +/* Round down to an integer, so that the sample is valid for a Poisson + distribution. */ + result = floor( result ); + +/* Calculate the ratio between the required Poisson distribution and + the Lorentzian from which we have sampled (the factor of 0.9 prevents + this exceeding 1.0, and overflow is again prevented by the checks + performed above). */ + ran *= sqrt_point9; + pfract = ( 0.9 + ran * ran ) * + exp( result * log_mean - LogGamma( result + 1.0, status ) - beta ); + +/* Accept the sample with this fractional probability, otherwise + obtain a new sample. */ + } while ( Rand( context, status ) > pfract ); + +/* If the mean is huge, the relative standard deviation will be + negligible compared to the machine precision. In such cases, the + probability of getting a result that differs from the mean is + effectively zero, so we can simply return the mean. */ + } else { + result = mean; + } + + UNLOCK_MUTEX6 + +/* Return the result. */ + return result; +} + +static double Rand( Rcontext *context, int *status ) { +/* +* Name: +* Rand + +* Purpose: +* Produce a uniformly distributed pseudo-random number. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* double Rand( Rcontext *context, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* On each invocation, this function returns a pseudo-random number +* uniformly distributed in the range 0.0 to 1.0 (inclusive). The +* underlying algorithm is that used by the "ran2" function of Press et +* al. (Numerical Recipes), which has a long period and good statistical +* properties. This independent implementation returns double precision +* values. + +* Parameters: +* context +* Pointer to an Rcontext structure which holds the random number +* generator's context between invocations. +* status +* Pointer to the inherited status variable. + +* Notes: +* - The sequence of numbers returned is determined by the "seed" +* value in the Rcontext structure supplied. +* - If the seed value is changed, the "active" flag must also be cleared +* so that this function can re-initiallise the Rcontext structure before +* generating the next pseudo-random number. The "active" flag should +* also be clear to force initialisation the first time an Rcontext +* structure is used. +* - This function does not perform error checking and does not generate +* errors. It will execute even if the global error status is set. +*/ + +/* Local Constants: */ + const long int a1 = 40014L; /* Random number generator constants... */ + const long int a2 = 40692L; + const long int m1 = 2147483563L; + const long int m2 = 2147483399L; + const long int q1 = 53668L; + const long int q2 = 52774L; + const long int r1 = 12211L; + const long int r2 = 3791L; + const int ntab = /* Size of shuffle table */ + AST_MATHMAP_RAND_CONTEXT_NTAB_; + const int nwarm = 8; /* Number of warm-up iterations */ + +/* Local Variables: */ + double result; /* Result value to return */ + double scale; /* Scale factor for random integers */ + double sum; /* Sum for forming normalisation constant */ + int dbits; /* Approximate bits in double mantissa */ + int irand; /* Loop counter for random integers */ + int itab; /* Loop counter for shuffle table */ + int lbits; /* Approximate bits used by generators */ + long int seed; /* Random number seed */ + long int tmp; /* Temporary variable */ + static double norm; /* Normalisation constant */ + static double scale0; /* Scale decrement for successive integers */ + static int init = 0; /* Local initialisation performed? */ + static int nrand; /* Number of random integers to use */ + +/* If the random number generator context is not active, then + initialise it. */ + if ( !context->active ) { + +/* First, perform local initialisation for this function, if not + already done. */ + LOCK_MUTEX4 + if ( !init ) { + +/* Obtain the approximate number of bits used by the random integer + generator from the value "m1". */ + (void) frexp( (double) m1, &lbits ); + +/* Obtain the approximate number of bits used by the mantissa of the + double value we want to produce, allowing for the (unlikely) + possibility that the mantissa's radix isn't 2. */ + dbits = (int) ceil( (double) DBL_MANT_DIG * + log( (double) FLT_RADIX ) / log( 2.0 ) ); + +/* Hence determine how many random integers we need to combine to + produce each double value, so that all the mantissa's bits will be + used. */ + nrand = ( dbits + lbits - 1 ) / lbits; + +/* Calculate the scale factor by which each successive random + integer's contribution to the result is reduced so as to generate + progressively less significant bits. */ + scale0 = 1.0 / (double) ( m1 - 1L ); + +/* Loop to sum the maximum contributions from each random integer + (assuming that each takes the largest possible value, of "m1-1", + from which we will later subtract 1). This produces the normalisation + factor by which the result must be scaled so as to lie between 0.0 and + 1.0 (inclusive). */ + sum = 0.0; + scale = 1.0; + for ( irand = 0; irand < nrand; irand++ ) { + scale *= scale0; + sum += scale; + } + norm = 1.0 / ( sum * (double) ( m1 - 2L ) ); + +/* Note that local initialisation has been done. */ + init = 1; + } + UNLOCK_MUTEX4 + +/* Obtain the seed value, enforcing positivity. */ + seed = (long int) context->seed; + if ( seed < 1 ) seed = seed + LONG_MAX; + if ( seed < 1 ) seed = LONG_MAX; + +/* Initialise the random number generators with this seed. */ + context->rand1 = context->rand2 = seed; + +/* Now loop to initialise the shuffle table with an initial set of + random values. We generate more values than required in order to "warm + up" the generator before recording values in the table. */ + for ( itab = ntab + nwarm - 1; itab >= 0; itab-- ) { + +/* Repeatedly update "rand1" from the expression "(rand1*a1)%m1" while + avoiding overflow. */ + tmp = context->rand1 / q1; + context->rand1 = a1 * ( context->rand1 - tmp * q1 ) - tmp * r1; + if ( context->rand1 < 0L ) context->rand1 += m1; + +/* After warming up, start recording values in the table. */ + if ( itab < ntab ) context->table[ itab ] = context->rand1; + } + +/* Record the last entry in the table as the "previous" random + integer. */ + context->random_int = context->table[ 0 ]; + +/* Note the random number generator context is active. */ + context->active = 1; + } + +/* Generate a random value. */ +/* ------------------------ */ +/* Initialise. */ + result = 0.0; + +/* Loop to generate sufficient random integers to combine into a + double value. */ + scale = norm; + for ( irand = 0; irand < nrand; irand++ ) { + +/* Update the first generator "rand1" from the expression + "(a1*rand1)%m1" while avoiding overflow. */ + tmp = context->rand1 / q1; + context->rand1 = a1 * ( context->rand1 - tmp * q1 ) - tmp * r1; + if ( context->rand1 < 0L ) context->rand1 += m1; + +/* Similarly, update the second generator "rand2" from the expression + "(a2*rand2)%m2". */ + tmp = context->rand2 / q2; + context->rand2 = a2 * ( context->rand2 - tmp * q2 ) - tmp * r2; + if ( context->rand2 < 0L ) context->rand2 += m2; + +/* Use the previous random integer to generate an index into the + shuffle table. */ + itab = (int) ( context->random_int / + ( 1L + ( m1 - 1L ) / (long int) ntab ) ); + +/* The algorithm left by RFWS seems to have a bug that "itab" can + sometimes be outside the range of [0.,ntab-1] causing the context->table + array to be addressed out of bounds. To avoid this, use the + following sticking plaster, since I'm not sure what the correct fix is. */ + if( itab < 0 ) itab = -itab; + itab = itab % ntab; + +/* Extract the table entry and replace it with a new random value from + the first generator "rand1". This is the Bays-Durham shuffle. */ + context->random_int = context->table[ itab ]; + context->table[ itab ] = context->rand1; + +/* Combine the extracted value with the latest value from the second + generator "rand2". */ + context->random_int -= context->rand2; + if ( context->random_int < 1L ) context->random_int += m1 - 1L; + +/* Update the scale factor to apply to the resulting random integer + and accumulate its contribution to the result. */ + scale *= scale0; + result += scale * (double) ( context->random_int - 1L ); + } + +/* Return the result. */ + return result; +} + +static void SetAttrib( AstObject *this_object, const char *setting, int *status ) { +/* +* Name: +* SetAttrib + +* Purpose: +* Set an attribute value for a MathMap. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* void SetAttrib( AstObject *this, const char *setting, int *status ) + +* Class Membership: +* MathMap member function (extends the astSetAttrib method inherited from +* the Mapping class). + +* Description: +* This function assigns an attribute value for a MathMap, the attribute +* and its value being specified by means of a string of the form: +* +* "attribute= value " +* +* Here, "attribute" specifies the attribute name and should be in lower +* case with no white space present. The value to the right of the "=" +* should be a suitable textual representation of the value to be assigned +* and this will be interpreted according to the attribute's data type. +* White space surrounding the value is only significant for string +* attributes. + +* Parameters: +* this +* Pointer to the MathMap. +* setting +* Pointer to a null terminated string specifying the new attribute +* value. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* void +*/ + +/* Local Vaiables: */ + AstMathMap *this; /* Pointer to the MathMap structure */ + int ival; /* Integer attribute value */ + int len; /* Length of setting string */ + int nc; /* Number of characters read by astSscanf */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Obtain a pointer to the MathMap structure. */ + this = (AstMathMap *) this_object; + +/* Obtain the length of the setting string. */ + len = strlen( setting ); + +/* Test for each recognised attribute in turn, using "astSscanf" to parse the + setting string and extract the attribute value (or an offset to it in the + case of string values). In each case, use the value set in "nc" to check + that the entire string was matched. Once a value has been obtained, use the + appropriate method to set it. */ + +/* Seed. */ +/* ----- */ + if ( nc = 0, + ( 1 == astSscanf( setting, "seed= %d %n", &ival, &nc ) ) + && ( nc >= len ) ) { + astSetSeed( this, ival ); + +/* SimpFI. */ +/* ------- */ + } else if ( nc = 0, + ( 1 == astSscanf( setting, "simpfi= %d %n", &ival, &nc ) ) + && ( nc >= len ) ) { + astSetSimpFI( this, ival ); + +/* SimpIF. */ +/* ------- */ + } else if ( nc = 0, + ( 1 == astSscanf( setting, "simpif= %d %n", &ival, &nc ) ) + && ( nc >= len ) ) { + astSetSimpIF( this, ival ); + +/* Pass any unrecognised setting to the parent method for further + interpretation. */ + } else { + (*parent_setattrib)( this_object, setting, status ); + } +} + +static int TestAttrib( AstObject *this_object, const char *attrib, int *status ) { +/* +* Name: +* TestAttrib + +* Purpose: +* Test if a specified attribute value is set for a MathMap. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* int TestAttrib( AstObject *this, const char *attrib, int *status ) + +* Class Membership: +* MathMap member function (over-rides the astTestAttrib protected +* method inherited from the Mapping class). + +* Description: +* This function returns a boolean result (0 or 1) to indicate whether +* a value has been set for one of a MathMap's attributes. + +* Parameters: +* this +* Pointer to the MathMap. +* attrib +* Pointer to a null terminated string specifying the attribute +* name. This should be in lower case with no surrounding white +* space. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* One if a value has been set, otherwise zero. + +* Notes: +* - A value of zero will be returned if this function is invoked +* with the global status set, or if it should fail for any reason. +*/ + +/* Local Variables: */ + AstMathMap *this; /* Pointer to the MathMap structure */ + int result; /* Result value to return */ + +/* Initialise. */ + result = 0; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Obtain a pointer to the MathMap structure. */ + this = (AstMathMap *) this_object; + +/* Check the attribute name and test the appropriate attribute. */ + +/* Seed. */ +/* ----- */ + if ( !strcmp( attrib, "seed" ) ) { + result = astTestSeed( this ); + +/* SimpFI. */ +/* ------- */ + } else if ( !strcmp( attrib, "simpfi" ) ) { + result = astTestSimpFI( this ); + +/* SimpIF. */ +/* ------- */ + } else if ( !strcmp( attrib, "simpif" ) ) { + result = astTestSimpIF( this ); + +/* If the attribute is not recognised, pass it on to the parent method + for further interpretation. */ + } else { + result = (*parent_testattrib)( this_object, attrib, status ); + } + +/* Return the result, */ + return result; +} + +static AstPointSet *Transform( AstMapping *map, AstPointSet *in, + int forward, AstPointSet *out, int *status ) { +/* +* Name: +* Transform + +* Purpose: +* Apply a MathMap to transform a set of points. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* AstPointSet *Transform( AstMapping *map, AstPointSet *in, +* int forward, AstPointSet *out, int *status ) + +* Class Membership: +* MathMap member function (over-rides the astTransform method inherited +* from the Mapping class). + +* Description: +* This function takes a MathMap and a set of points encapsulated in a +* PointSet and transforms the points so as to apply the required coordinate +* transformation. + +* Parameters: +* map +* Pointer to the MathMap. +* in +* Pointer to the PointSet holding the input coordinate data. +* forward +* A non-zero value indicates that the forward coordinate transformation +* should be applied, while a zero value requests the inverse +* transformation. +* out +* Pointer to a PointSet which will hold the transformed (output) +* coordinate values. A NULL value may also be given, in which case a +* new PointSet will be created by this function. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Pointer to the output (possibly new) PointSet. + +* Notes: +* - A null pointer will be returned if this function is invoked with the +* global error status set, or if it should fail for any reason. +* - The number of coordinate values per point in the input PointSet must +* match the number of coordinates for the MathMap being applied. +* - If an output PointSet is supplied, it must have space for sufficient +* number of points and coordinate values per point to accommodate the +* result. Any excess space will be ignored. +*/ + +/* Local Variables: */ + AstMathMap *this; /* Pointer to MathMap to be applied */ + AstPointSet *result; /* Pointer to output PointSet */ + double **data_ptr; /* Array of pointers to coordinate data */ + double **ptr_in; /* Pointer to input coordinate data */ + double **ptr_out; /* Pointer to output coordinate data */ + double *work; /* Workspace for intermediate results */ + int idata; /* Loop counter for data pointer elements */ + int ifun; /* Loop counter for functions */ + int ncoord_in; /* Number of coordinates per input point */ + int ncoord_out; /* Number of coordinates per output point */ + int ndata; /* Number of data pointer elements filled */ + int nfun; /* Number of functions to evaluate */ + int npoint; /* Number of points */ + +/* Check the global error status. */ + if ( !astOK ) return NULL; + +/* Initialise variables to avoid "used of uninitialised variable" + messages from dumb compilers. */ + work = NULL; + +/* Obtain a pointer to the MathMap. */ + this = (AstMathMap *) map; + +/* Apply the parent mapping using the stored pointer to the Transform member + function inherited from the parent Mapping class. This function validates + all arguments and generates an output PointSet if necessary, but does not + actually transform any coordinate values. */ + result = (*parent_transform)( map, in, forward, out, status ); + +/* We will now extend the parent astTransform method by performing the + transformation needed to generate the output coordinate values. */ + +/* Determine the numbers of points and coordinates per point from the input + and output PointSets and obtain pointers for accessing the input and output + coordinate values. */ + ncoord_in = astGetNcoord( in ); + ncoord_out = astGetNcoord( result ); + npoint = astGetNpoint( in ); + ptr_in = astGetPoints( in ); + ptr_out = astGetPoints( result ); + +/* Determine whether to apply the forward or inverse transformation, according + to the direction specified and whether the mapping has been inverted. */ + if ( astGetInvert( this ) ) forward = !forward; + +/* Obtain the number of transformation functions that must be + evaluated to perform the transformation. This will include any that + produce intermediate results from which the final results are + calculated. */ + nfun = forward ? this->nfwd : this->ninv; + +/* If intermediate results are to be calculated, then allocate + workspace to hold them (each intermediate result being a vector of + "npoint" double values). */ + if ( nfun > ncoord_out ) { + work = astMalloc( sizeof( double) * + (size_t) ( npoint * ( nfun - ncoord_out ) ) ); + } + +/* Also allocate space for an array to hold pointers to the input + data, intermediate results and output data. */ + data_ptr = astMalloc( sizeof( double * ) * (size_t) ( ncoord_in + nfun ) ); + +/* We now set up the "data_ptr" array to locate the data to be + processed. */ + if ( astOK ) { + +/* The first elements of this array point at the input data + vectors. */ + ndata = 0; + for ( idata = 0; idata < ncoord_in; idata++ ) { + data_ptr[ ndata++ ] = ptr_in[ idata ]; + } + +/* The following elements point at successive vectors within the + workspace array (if allocated). These vectors will act first as output + arrays for intermediate results, and then as input arrays for + subsequent calculations which use these results. */ + for ( idata = 0; idata < ( nfun - ncoord_out ); idata++ ) { + data_ptr[ ndata++ ] = work + ( idata * npoint ); + } + +/* The final elements point at the output coordinate data arrays into + which the final results will be written. */ + for ( idata = 0; idata < ncoord_out; idata++ ) { + data_ptr[ ndata++ ] = ptr_out[ idata ]; + } + +/* Perform coordinate transformation. */ +/* ---------------------------------- */ +/* Loop to evaluate each transformation function in turn. */ + for ( ifun = 0; ifun < nfun; ifun++ ) { + +/* Invoke the function that evaluates compiled expressions. Pass the + appropriate code and constants arrays, depending on the direction of + coordinate transformation, together with the required stack size. The + output array is the vector located by successive elements of the + "data_ptr" array (skipping the input data elements), while the + function has access to all previous elements of the "data_ptr" array + to locate the required input data. */ + EvaluateFunction( &this->rcontext, npoint, (const double **) data_ptr, + forward ? this->fwdcode[ ifun ] : + this->invcode[ ifun ], + forward ? this->fwdcon[ ifun ] : + this->invcon[ ifun ], + forward ? this->fwdstack : this->invstack, + data_ptr[ ifun + ncoord_in ], status ); + } + } + +/* Free the array of data pointers and any workspace allocated for + intermediate results. */ + data_ptr = astFree( data_ptr ); + if ( nfun > ncoord_out ) work = astFree( work ); + +/* If an error occurred, then return a NULL pointer. If no output + PointSet was supplied, also delete any new one that may have been + created. */ + if ( !astOK ) { + result = ( result == out ) ? NULL : astDelete( result ); + } + +/* Return a pointer to the output PointSet. */ + return result; +} + +static void ValidateSymbol( const char *method, const char *class, + const char *exprs, int iend, int sym, + int *lpar, int **argcount, int **opensym, + int *ncon, double **con, int *status ) { +/* +* Name: +* ValidateSymbol + +* Purpose: +* Validate a symbol in an expression. + +* Type: +* Private function. + +* Synopsis: +* #include "mathmap.h" +* void ValidateSymbol( const char *method, const char *class, +* const char *exprs, int iend, int sym, int *lpar, +* int **argcount, int **opensym, int *ncon, +* double **con, int *status ) + +* Class Membership: +* MathMap member function. + +* Description: +* This function validates an identified standard symbol during +* compilation of an expression. Its main task is to keep track of the +* level of parenthesis in the expression and to count the number of +* arguments supplied to functions at each level of parenthesis (for +* nested function calls). On this basis it is able to interpret and +* accept or reject symbols which represent function calls, parentheses +* and delimiters. Other symbols are accepted automatically. + +* Parameters: +* method +* Pointer to a constant null-terminated character string +* containing the name of the method that invoked this function. +* This method name is used solely for constructing error messages. +* class +* Pointer to a constant null-terminated character string containing the +* class name of the Object being processed. This name is used solely +* for constructing error messages. +* exprs +* Pointer to a null-terminated string containing the expression +* being parsed. This is only used for constructing error messages. +* iend +* Index in "exprs" of the last character belonging to the most +* recently identified symbol. This is only used for constructing error +* messages. +* sym +* Index in the static "symbol" array of the most recently identified +* symbol in the expression. This is the symbol to be verified. +* lpar +* Pointer to an int which holds the current level of parenthesis. On +* the first invocation, this should be zero. The returned value should +* be passed to subsequent invocations. +* argcount +* Address of a pointer to a dynamically allocated array of int in +* which argument count information is maintained for each level of +* parenthesis (e.g. for nested function calls). On the first invocation, +* "*argcount" should be NULL. This function will allocate the required +* space as needed and update this pointer. The returned pointer value +* should be passed to subsequent invocations. +* +* The allocated space must be freed by the caller (using astFree) when +* no longer required. +* opensym +* Address of a pointer to a dynamically allocated array of int, in which +* information is maintained about the functions associated with each +* level of parenthesis (e.g. for nested function calls). On the first +* invocation, "*opensym" should be NULL. This function will allocate the +* required space as needed and update this pointer. The returned pointer +* value should be passed to subsequent invocations. +* +* The allocated space must be freed by the caller (using astFree) when +* no longer required. +* ncon +* Pointer to an int which holds a count of the constants associated +* with the expression (and determines the size of the "*con" array). +* This function will update the count to reflect any new constants +* appended to the "*con" array and the returned value should be passed +* to subsequent invocations. +* con +* Address of a pointer to a dynamically allocated array of double, in +* which the constants associated with the expression being parsed are +* accumulated. On entry, "*con" should point at a dynamic array with +* at least "*ncon" elements containing existing constants (or may be +* NULL if no constants have yet been stored). This function will +* allocate the required space as needed and update this pointer (and +* "*ncon") appropriately. The returned pointer value should be passed +* to subsequent invocations. +* +* The allocated space must be freed by the caller (using astFree) when +* no longer required. +* status +* Pointer to the inherited status variable. + +* Notes: +* - The dynamically allocated arrays normally returned by this function +* will be freed and NULL pointers will be returned if this function is +* invoked with the global error status set, or if it should fail for any +* reason. +*/ + +/* Check the global error status, but do not return at this point + because dynamic arrays may require freeing. */ + if ( astOK ) { + +/* Check if the symbol is a comma. */ + if ( ( symbol[ sym ].text[ 0 ] == ',' ) && + ( symbol[ sym ].text[ 1 ] == '\0' ) ) { + +/* A comma is only used to delimit function arguments. If the current + level of parenthesis is zero, or the symbol which opened the current + level of parenthesis was not a function call (indicated by an argument + count of zero at the current level of parenthesis), then report an + error. */ + if ( ( *lpar <= 0 ) || ( ( *argcount )[ *lpar - 1 ] == 0 ) ) { + astError( AST__COMIN, + "%s(%s): Spurious comma encountered in the expression " + "\"%.*s\".", status, + method, class, iend + 1, exprs ); + +/* If a comma is valid, then increment the argument count at the + current level of parenthesis. */ + } else { + ( *argcount )[ *lpar - 1 ]++; + } + +/* If the symbol is not a comma, check if it increases the current + level of parenthesis. */ + } else if ( symbol[ sym ].parincrement > 0 ) { + +/* Increase the size of the arrays which hold parenthesis level + information and check for errors. */ + *argcount = astGrow( *argcount, *lpar + 1, sizeof( int ) ); + *opensym = astGrow( *opensym, *lpar + 1, sizeof( int ) ); + if ( astOK ) { + +/* Increment the level of parenthesis and initialise the argument + count at the new level. This count is set to zero if the symbol which + opens the parenthesis level is not a function call (indicated by a + zero "nargs" entry in the symbol data), and it subsequently remains at + zero. If the symbol is a function call, the argument count is + initially set to 1 and increments whenever a comma is encountered at + this parenthesis level. */ + ( *argcount )[ ++( *lpar ) - 1 ] = ( symbol[ sym ].nargs != 0 ); + +/* Remember the symbol which opened this parenthesis level. */ + ( *opensym )[ *lpar - 1 ] = sym; + } + +/* Check if the symbol decreases the current parenthesis level. */ + } else if ( symbol[ sym ].parincrement < 0 ) { + +/* Ensure that the parenthesis level is not already at zero. If it is, + then there is a missing left parenthesis in the expression being + compiled, so report an error. */ + if ( *lpar == 0 ) { + astError( AST__MLPAR, + "%s(%s): Missing left parenthesis in the expression " + "\"%.*s\".", status, + method, class, iend + 1, exprs ); + +/* If the parenthesis level is valid and the symbol which opened this + level of parenthesis was a function call with a fixed number of + arguments (indicated by a positive "nargs" entry in the symbol data), + then we must check the number of function arguments which have been + encountered. */ + } else if ( symbol[ ( *opensym )[ *lpar - 1 ] ].nargs > 0 ) { + +/* Report an error if the number of arguments is wrong. */ + if ( ( *argcount )[ *lpar - 1 ] != + symbol[ ( *opensym )[ *lpar - 1 ] ].nargs ) { + astError( AST__WRNFA, + "%s(%s): Wrong number of function arguments in the " + "expression \"%.*s\".", status, + method, class, iend + 1, exprs ); + +/* If the number of arguments is valid, decrement the parenthesis + level. */ + } else { + ( *lpar )--; + } + +/* If the symbol which opened this level of parenthesis was a function + call with a variable number of arguments (indicated by a negative + "nargs" entry in the symbol data), then we must check and process the + number of function arguments. */ + } else if ( symbol[ ( *opensym )[ *lpar - 1 ] ].nargs < 0 ) { + +/* Check that the minimum required number of arguments have been + supplied. Report an error if they have not. */ + if ( ( *argcount )[ *lpar - 1 ] < + ( -symbol[ ( *opensym )[ *lpar - 1 ] ].nargs ) ) { + astError( AST__WRNFA, + "%s(%s): Insufficient function arguments in the " + "expression \"%.*s\".", status, + method, class, iend + 1, exprs ); + +/* If the number of arguments is valid, increase the size of the + constants array and check for errors. */ + } else { + *con = astGrow( *con, *ncon + 1, sizeof( double ) ); + if ( astOK ) { + +/* Append the argument count to the end of the array of constants and + decrement the parenthesis level. */ + ( *con )[ ( *ncon )++ ] = + (double) ( *argcount )[ --( *lpar ) ]; + } + } + +/* Finally, if the symbol which opened this level of parenthesis was + not a function call ("nargs" entry in the symbol data is zero), then + decrement the parenthesis level. In this case there is no need to + check the argument count, because it will not have been + incremented. */ + } else { + ( *lpar )--; + } + } + } + +/* If an error occurred (or the global error status was set on entry), + then reset the parenthesis level and free any memory which may have + been allocated. */ + if ( !astOK ) { + *lpar = 0; + if ( *argcount ) *argcount = astFree( *argcount ); + if ( *opensym ) *opensym = astFree( *opensym ); + if ( *con ) *con = astFree( *con ); + } +} + +/* Functions which access class attributes. */ +/* ---------------------------------------- */ +/* Implement member functions to access the attributes associated with + this class using the macros defined for this purpose in the + "object.h" file. For a description of each attribute, see the class + interface (in the associated .h file). */ + +/* +*att++ +* Name: +* Seed + +* Purpose: +* Random number seed for a MathMap. + +* Type: +* Public attribute. + +* Synopsis: +* Integer. + +* Description: +* This attribute, which may take any integer value, determines the +* sequence of random numbers produced by the random number functions in +* MathMap expressions. It is set to an unpredictable default value when +* a MathMap is created, so that by default each MathMap uses a different +* set of random numbers. +* +* If required, you may set this Seed attribute to a value of your +* choosing in order to produce repeatable behaviour from the random +* number functions. You may also enquire the Seed value (e.g. if an +* initially unpredictable value has been used) and then use it to +* reproduce the resulting sequence of random numbers, either from the +* same MathMap or from another one. +* +* Clearing the Seed attribute gives it a new unpredictable default +* value. + +* Applicability: +* MathMap +* All MathMaps have this attribute. +*att-- +*/ +/* Clear the Seed value by setting it to a new unpredictable value + produced by DefaultSeed and clearing the "seed_set" flag in the + MathMap's random number generator context. Also clear the "active" + flag, so that the generator will be re-initialised to use this seed + when it is next invoked. */ +astMAKE_CLEAR(MathMap,Seed,rcontext.seed,( this->rcontext.seed_set = 0, + this->rcontext.active = 0, + DefaultSeed( &this->rcontext, status ) )) + +/* Return the "seed" value from the random number generator + context. */ +astMAKE_GET(MathMap,Seed,int,0,this->rcontext.seed) + +/* Store the new seed value in the MathMap's random number generator + context and set the context's "seed_set" flag. Also clear the "active" + flag, so that the generator will be re-initialised to use this seed + when it is next invoked. */ +astMAKE_SET(MathMap,Seed,int,rcontext.seed,( this->rcontext.seed_set = 1, + this->rcontext.active = 0, + value )) + +/* Test the "seed_set" flag in the random number generator context. */ +astMAKE_TEST(MathMap,Seed,( this->rcontext.seed_set )) + +/* +*att++ +* Name: +* SimpFI + +* Purpose: +* Forward-inverse MathMap pairs simplify? + +* Type: +* Public attribute. + +* Synopsis: +* Integer (boolean). + +* Description: +c This attribute should be set to a non-zero value if applying a +c MathMap's forward transformation, followed immediately by the matching +c inverse transformation will always restore the original set of +c coordinates. It indicates that AST may replace such a sequence of +c operations by an identity Mapping (a UnitMap) if it is encountered +c while simplifying a compound Mapping (e.g. using astSimplify). +f This attribute should be set to a non-zero value if applying a +f MathMap's forward transformation, followed immediately by the matching +f inverse transformation will always restore the original set of +f coordinates. It indicates that AST may replace such a sequence of +f operations by an identity Mapping (a UnitMap) if it is encountered +f while simplifying a compound Mapping (e.g. using AST_SIMPLIFY). +* +* By default, the SimpFI attribute is zero, so that AST will not perform +* this simplification unless you have set SimpFI to indicate that it is +* safe to do so. + +* Applicability: +* MathMap +* All MathMaps have this attribute. + +* Notes: +* - For simplification to occur, the two MathMaps must be in series and +* be identical (with textually identical transformation +* functions). Functional equivalence is not sufficient. +* - The consent of both MathMaps is required before simplification can +* take place. If either has a SimpFI value of zero, then simplification +* will not occur. +* - The SimpFI attribute controls simplification only in the case where +* a MathMap's forward transformation is followed by the matching inverse +* transformation. It does not apply if an inverse transformation is +* followed by a forward transformation. This latter case is controlled +* by the SimpIF attribute. +c - The "forward" and "inverse" transformations referred to are those +c defined when the MathMap is created (corresponding to the "fwd" and +c "inv" parameters of its constructor function). If the MathMap is +c inverted (i.e. its Invert attribute is non-zero), then the role of the +c SimpFI and SimpIF attributes will be interchanged. +f - The "forward" and "inverse" transformations referred to are those +f defined when the MathMap is created (corresponding to the FWD and +f INV arguments of its constructor function). If the MathMap is +f inverted (i.e. its Invert attribute is non-zero), then the role of the +f SimpFI and SimpIF attributes will be interchanged. +*att-- +*/ +/* Clear the SimpFI value by setting it to -INT_MAX. */ +astMAKE_CLEAR(MathMap,SimpFI,simp_fi,-INT_MAX) + +/* Supply a default of 0 if no SimpFI value has been set. */ +astMAKE_GET(MathMap,SimpFI,int,0,( ( this->simp_fi != -INT_MAX ) ? + this->simp_fi : 0 )) + +/* Set a SimpFI value of 1 if any non-zero value is supplied. */ +astMAKE_SET(MathMap,SimpFI,int,simp_fi,( value != 0 )) + +/* The SimpFI value is set if it is not -INT_MAX. */ +astMAKE_TEST(MathMap,SimpFI,( this->simp_fi != -INT_MAX )) + +/* +*att++ +* Name: +* SimpIF + +* Purpose: +* Inverse-forward MathMap pairs simplify? + +* Type: +* Public attribute. + +* Synopsis: +* Integer (boolean). + +* Description: +c This attribute should be set to a non-zero value if applying a +c MathMap's inverse transformation, followed immediately by the matching +c forward transformation will always restore the original set of +c coordinates. It indicates that AST may replace such a sequence of +c operations by an identity Mapping (a UnitMap) if it is encountered +c while simplifying a compound Mapping (e.g. using astSimplify). +f This attribute should be set to a non-zero value if applying a +f MathMap's inverse transformation, followed immediately by the matching +f forward transformation will always restore the original set of +f coordinates. It indicates that AST may replace such a sequence of +f operations by an identity Mapping (a UnitMap) if it is encountered +f while simplifying a compound Mapping (e.g. using AST_SIMPLIFY). +* +* By default, the SimpIF attribute is zero, so that AST will not perform +* this simplification unless you have set SimpIF to indicate that it is +* safe to do so. + +* Applicability: +* MathMap +* All MathMaps have this attribute. + +* Notes: +* - For simplification to occur, the two MathMaps must be in series and +* be identical (with textually identical transformation +* functions). Functional equivalence is not sufficient. +* - The consent of both MathMaps is required before simplification can +* take place. If either has a SimpIF value of zero, then simplification +* will not occur. +* - The SimpIF attribute controls simplification only in the case where +* a MathMap's inverse transformation is followed by the matching forward +* transformation. It does not apply if a forward transformation is +* followed by an inverse transformation. This latter case is controlled +* by the SimpFI attribute. +c - The "forward" and "inverse" transformations referred to are those +c defined when the MathMap is created (corresponding to the "fwd" and +c "inv" parameters of its constructor function). If the MathMap is +c inverted (i.e. its Invert attribute is non-zero), then the role of the +c SimpFI and SimpIF attributes will be interchanged. +f - The "forward" and "inverse" transformations referred to are those +f defined when the MathMap is created (corresponding to the FWD and +f INV arguments of its constructor function). If the MathMap is +f inverted (i.e. its Invert attribute is non-zero), then the role of the +f SimpFI and SimpIF attributes will be interchanged. +*att-- +*/ +/* Clear the SimpIF value by setting it to -INT_MAX. */ +astMAKE_CLEAR(MathMap,SimpIF,simp_if,-INT_MAX) + +/* Supply a default of 0 if no SimpIF value has been set. */ +astMAKE_GET(MathMap,SimpIF,int,0,( ( this->simp_if != -INT_MAX ) ? + this->simp_if : 0 )) + +/* Set a SimpIF value of 1 if any non-zero value is supplied. */ +astMAKE_SET(MathMap,SimpIF,int,simp_if,( value != 0 )) + +/* The SimpIF value is set if it is not -INT_MAX. */ +astMAKE_TEST(MathMap,SimpIF,( this->simp_if != -INT_MAX )) + +/* Copy constructor. */ +/* ----------------- */ +static void Copy( const AstObject *objin, AstObject *objout, int *status ) { +/* +* Name: +* Copy + +* Purpose: +* Copy constructor for MathMap objects. + +* Type: +* Private function. + +* Synopsis: +* void Copy( const AstObject *objin, AstObject *objout, int *status ) + +* Description: +* This function implements the copy constructor for MathMap objects. + +* Parameters: +* objin +* Pointer to the object to be copied. +* objout +* Pointer to the object being constructed. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* void + +* Notes: +* - This constructor makes a deep copy. +*/ + +/* Local Variables: */ + AstMathMap *in; /* Pointer to input MathMap */ + AstMathMap *out; /* Pointer to output MathMap */ + int ifun; /* Loop counter for functions */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Obtain pointers to the input and output MathMaps. */ + in = (AstMathMap *) objin; + out = (AstMathMap *) objout; + +/* For safety, first clear any references to the input memory from + the output MathMap. */ + out->fwdfun = NULL; + out->invfun = NULL; + out->fwdcode = NULL; + out->invcode = NULL; + out->fwdcon = NULL; + out->invcon = NULL; + +/* Now allocate and initialise each of the output pointer arrays + required. */ + if ( in->fwdfun ) { + MALLOC_POINTER_ARRAY( out->fwdfun, char *, out->nfwd ) + } + if ( in->invfun ) { + MALLOC_POINTER_ARRAY( out->invfun, char *, out->ninv ) + } + if ( in->fwdcode ) { + MALLOC_POINTER_ARRAY( out->fwdcode, int *, out->nfwd ) + } + if ( in->invcode ) { + MALLOC_POINTER_ARRAY( out->invcode, int *, out->ninv ) + } + if ( in->fwdcon ) { + MALLOC_POINTER_ARRAY( out->fwdcon, double *, out->nfwd ) + } + if ( in->invcon ) { + MALLOC_POINTER_ARRAY( out->invcon, double *, out->ninv ) + } + +/* If OK, loop to make copies of the data (where available) associated + with each forward transformation function, storing pointers to the + copy in the output pointer arrays allocated above. */ + if ( astOK ) { + for ( ifun = 0; ifun < out->nfwd; ifun++ ) { + if ( in->fwdfun && in->fwdfun[ ifun ] ) { + out->fwdfun[ ifun ] = astStore( NULL, in->fwdfun[ ifun ], + astSizeOf( in->fwdfun[ ifun ] ) ); + } + if ( in->fwdcode && in->fwdcode[ ifun ] ) { + out->fwdcode[ ifun ] = astStore( NULL, in->fwdcode[ ifun ], + astSizeOf( in->fwdcode[ ifun ] ) ); + } + if ( in->fwdcon && in->fwdcon[ ifun ] ) { + out->fwdcon[ ifun ] = astStore( NULL, in->fwdcon[ ifun ], + astSizeOf( in->fwdcon[ ifun ] ) ); + } + if ( !astOK ) break; + } + } + +/* Repeat this process for the inverse transformation functions. */ + if ( astOK ) { + for ( ifun = 0; ifun < out->ninv; ifun++ ) { + if ( in->invfun && in->invfun[ ifun ] ) { + out->invfun[ ifun ] = astStore( NULL, in->invfun[ ifun ], + astSizeOf( in->invfun[ ifun ] ) ); + } + if ( in->invcode && in->invcode[ ifun ] ) { + out->invcode[ ifun ] = astStore( NULL, in->invcode[ ifun ], + astSizeOf( in->invcode[ ifun ] ) ); + } + if ( in->invcon && in->invcon[ ifun ] ) { + out->invcon[ ifun ] = astStore( NULL, in->invcon[ ifun ], + astSizeOf( in->invcon[ ifun ] ) ); + } + if ( !astOK ) break; + } + } + +/* If an error occurred, clean up by freeing all output memory + allocated above. */ + if ( !astOK ) { + FREE_POINTER_ARRAY( out->fwdfun, out->nfwd ) + FREE_POINTER_ARRAY( out->invfun, out->ninv ) + FREE_POINTER_ARRAY( out->fwdcode, out->nfwd ) + FREE_POINTER_ARRAY( out->invcode, out->ninv ) + FREE_POINTER_ARRAY( out->fwdcon, out->nfwd ) + FREE_POINTER_ARRAY( out->invcon, out->ninv ) + } +} + +/* Destructor. */ +/* ----------- */ +static void Delete( AstObject *obj, int *status ) { +/* +* Name: +* Delete + +* Purpose: +* Destructor for MathMap objects. + +* Type: +* Private function. + +* Synopsis: +* void Delete( AstObject *obj, int *status ) + +* Description: +* This function implements the destructor for MathMap objects. + +* Parameters: +* obj +* Pointer to the object to be deleted. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* void + +* Notes: +* This function attempts to execute even if the global error status is +* set. +*/ + +/* Local Variables: */ + AstMathMap *this; /* Pointer to MathMap */ + +/* Obtain a pointer to the MathMap structure. */ + this = (AstMathMap *) obj; + +/* Free all memory allocated by the MathMap. */ + FREE_POINTER_ARRAY( this->fwdfun, this->nfwd ) + FREE_POINTER_ARRAY( this->invfun, this->ninv ) + FREE_POINTER_ARRAY( this->fwdcode, this->nfwd ) + FREE_POINTER_ARRAY( this->invcode, this->ninv ) + FREE_POINTER_ARRAY( this->fwdcon, this->nfwd ) + FREE_POINTER_ARRAY( this->invcon, this->ninv ) +} + +/* Dump function. */ +/* -------------- */ +static void Dump( AstObject *this_object, AstChannel *channel, int *status ) { +/* +* Name: +* Dump + +* Purpose: +* Dump function for MathMap objects. + +* Type: +* Private function. + +* Synopsis: +* void Dump( AstObject *this, AstChannel *channel, int *status ) + +* Description: +* This function implements the Dump function which writes out data +* for the MathMap class to an output Channel. + +* Parameters: +* this +* Pointer to the MathMap whose data are being written. +* channel +* Pointer to the Channel to which the data are being written. +* status +* Pointer to the inherited status variable. +*/ + +/* Local Constants: */ +#define COMMENT_LEN 150 /* Maximum length of a comment string */ +#define KEY_LEN 50 /* Maximum length of a keyword */ + +/* Local Variables: */ + AstMathMap *this; /* Pointer to the MathMap structure */ + char comment[ COMMENT_LEN + 1 ]; /* Buffer for comment strings */ + char key[ KEY_LEN + 1 ]; /* Buffer for keyword strings */ + int ifun; /* Loop counter for functions */ + int invert; /* MathMap inverted? */ + int ival; /* Integer attribute value */ + int nin; /* True number of input coordinates */ + int nout; /* True number of output coordinates */ + int set; /* Attribute value set? */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Obtain a pointer to the MathMap structure. */ + this = (AstMathMap *) this_object; + +/* Determine if the MathMap is inverted and obtain the "true" number + of input and output coordinates by un-doing the effects of any + inversion. */ + invert = astGetInvert( this ); + nin = !invert ? astGetNin( this ) : astGetNout( this ); + nout = !invert ? astGetNout( this ) : astGetNin( this ); + +/* Write out values representing the instance variables for the + MathMap class. Accompany these with appropriate comment strings, + possibly depending on the values being written.*/ + +/* In the case of attributes, we first use the appropriate (private) + Test... member function to see if they are set. If so, we then use + the (private) Get... function to obtain the value to be written + out. + + For attributes which are not set, we use the astGet... method to + obtain the value instead. This will supply a default value + (possibly provided by a derived class which over-rides this method) + which is more useful to a human reader as it corresponds to the + actual default attribute value. Since "set" will be zero, these + values are for information only and will not be read back. */ + +/* Number of forward transformation functions. */ +/* ------------------------------------------- */ +/* We regard this value as set if it differs from the number of output + coordinates for the MathMap. */ + set = ( this->nfwd != nout ); + astWriteInt( channel, "Nfwd", set, 0, this->nfwd, + "Number of forward transformation functions" ); + +/* Forward transformation functions. */ +/* --------------------------------- */ +/* Loop to write out each forward transformation function, generating + a suitable keyword and comment for each one. */ + for ( ifun = 0; ifun < this->nfwd; ifun++ ) { + (void) sprintf( key, "Fwd%d", ifun + 1 ); + (void) sprintf( comment, "Forward function %d", ifun + 1 ); + astWriteString( channel, key, 1, 1, this->fwdfun[ ifun ], comment ); + } + +/* Number of inverse transformation functions. */ +/* ------------------------------------------- */ +/* We regard this value as set if it differs from the number of input + coordinates for the MathMap. */ + set = ( this->ninv != nin ); + astWriteInt( channel, "Ninv", set, 0, this->ninv, + "Number of inverse transformation functions" ); + +/* Inverse transformation functions. */ +/* --------------------------------- */ +/* Similarly, loop to write out each inverse transformation + function. */ + for ( ifun = 0; ifun < this->ninv; ifun++ ) { + (void) sprintf( key, "Inv%d", ifun + 1 ); + (void) sprintf( comment, "Inverse function %d", ifun + 1 ); + astWriteString( channel, key, 1, 1, this->invfun[ ifun ], comment ); + } + +/* SimpFI. */ +/* ------- */ +/* Write out the forward-inverse simplification flag. */ + set = TestSimpFI( this, status ); + ival = set ? GetSimpFI( this, status ) : astGetSimpFI( this ); + astWriteInt( channel, "SimpFI", set, 0, ival, + ival ? "Forward-inverse pairs may simplify" : + "Forward-inverse pairs do not simplify" ); + +/* SimpIF. */ +/* ------- */ +/* Write out the inverse-forward simplification flag. */ + set = TestSimpIF( this, status ); + ival = set ? GetSimpIF( this, status ) : astGetSimpIF( this ); + astWriteInt( channel, "SimpIF", set, 0, ival, + ival ? "Inverse-forward pairs may simplify" : + "Inverse-forward pairs do not simplify" ); + +/* Seed. */ +/* ----- */ +/* Write out any random number seed value which is set. Prefix this with + a separate flag which indicates if the seed has been set. */ + set = TestSeed( this, status ); + ival = set ? GetSeed( this, status ) : astGetSeed( this ); + astWriteInt( channel, "Seeded", set, 0, set, + set? "Explicit random number seed set" : + "No random number seed set" ); + astWriteInt( channel, "Seed", set, 0, ival, + set ? "Random number seed value" : + "Default random number seed used" ); + +/* Undefine macros local to this function. */ +#undef COMMENT_LEN +#undef KEY_LEN +} + +/* Standard class functions. */ +/* ========================= */ +/* Implement the astIsAMathMap and astCheckMathMap functions using the macros + defined for this purpose in the "object.h" header file. */ +astMAKE_ISA(MathMap,Mapping) +astMAKE_CHECK(MathMap) + +AstMathMap *astMathMap_( int nin, int nout, + int nfwd, const char *fwd[], + int ninv, const char *inv[], + const char *options, int *status, ...) { +/* +*+ +* Name: +* astMathMap + +* Purpose: +* Create a MathMap. + +* Type: +* Protected function. + +* Synopsis: +* #include "mathmap.h" +* AstMathMap *astMathMap( int nin, int nout, +* int nfwd, const char *fwd[], +* int ninv, const char *inv[], +* const char *options, ..., int *status ) + +* Class Membership: +* MathMap constructor. + +* Description: +* This function creates a new MathMap and optionally initialises its +* attributes. + +* Parameters: +* nin +* Number of input variables for the MathMap. +* nout +* Number of output variables for the MathMap. +* nfwd +* The number of forward transformation functions being supplied. +* This must be at least equal to "nout". +* fwd +* Pointer to an array, with "nfwd" elements, of pointers to null +* terminated strings which contain each of the forward transformation +* functions. +* ninv +* The number of inverse transformation functions being supplied. +* This must be at least equal to "nin". +* inv +* Pointer to an array, with "ninv" elements, of pointers to null +* terminated strings which contain each of the inverse transformation +* functions. +* options +* Pointer to a null terminated string containing an optional +* comma-separated list of attribute assignments to be used for +* initialising the new MathMap. The syntax used is the same as +* for the astSet method and may include "printf" format +* specifiers identified by "%" symbols in the normal way. +* status +* Pointer to the inherited status variable. +* ... +* If the "options" string contains "%" format specifiers, then +* an optional list of arguments may follow it in order to +* supply values to be substituted for these specifiers. The +* rules for supplying these are identical to those for the +* astSet method (and for the C "printf" function). + +* Returned Value: +* A pointer to the new MathMap. + +* Notes: +* - A NULL pointer will be returned if this function is invoked +* with the global error status set, or if it should fail for any +* reason. +*- + +* Implementation Notes: +* - This function implements the basic MathMap constructor which is +* available via the protected interface to the MathMap class. A +* public interface is provided by the astMathMapId_ function. +*/ + +/* Local Variables: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ + AstMathMap *new; /* Pointer to new MathMap */ + va_list args; /* Variable argument list */ + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(NULL); + +/* Check the global status. */ + if ( !astOK ) return NULL; + +/* Initialise the MathMap, allocating memory and initialising the + virtual function table as well if necessary. */ + new = astInitMathMap( NULL, sizeof( AstMathMap ), !class_init, &class_vtab, + "MathMap", nin, nout, nfwd, fwd, ninv, inv ); + +/* If successful, note that the virtual function table has been + initialised. */ + if ( astOK ) { + class_init = 1; + +/* Obtain the variable argument list and pass it along with the options string + to the astVSet method to initialise the new MathMap's attributes. */ + va_start( args, status ); + astVSet( new, options, NULL, args ); + va_end( args ); + +/* If an error occurred, clean up by deleting the new object. */ + if ( !astOK ) new = astDelete( new ); + } + +/* Return a pointer to the new MathMap. */ + return new; +} + +AstMathMap *astMathMapId_( int nin, int nout, + int nfwd, const char *fwd[], + int ninv, const char *inv[], + const char *options, ... ) { +/* +*++ +* Name: +c astMathMap +f AST_MATHMAP + +* Purpose: +* Create a MathMap. + +* Type: +* Public function. + +* Synopsis: +c #include "mathmap.h" +c AstMathMap *astMathMap( int nin, int nout, +c int nfwd, const char *fwd[], +c int ninv, const char *inv[], +c const char *options, ... ) +f RESULT = AST_MATHMAP( NIN, NOUT, NFWD, FWD, NINV, INV, OPTIONS, STATUS ) + +* Class Membership: +* MathMap constructor. + +* Description: +* This function creates a new MathMap and optionally initialises its +* attributes. +* +c A MathMap is a Mapping which allows you to specify a set of forward +c and/or inverse transformation functions using arithmetic operations +c and mathematical functions similar to those available in C. The +c MathMap interprets these functions at run-time, whenever its forward +c or inverse transformation is required. Because the functions are not +c compiled in the normal sense (unlike an IntraMap), they may be used to +c describe coordinate transformations in a transportable manner. A +c MathMap therefore provides a flexible way of defining new types of +c Mapping whose descriptions may be stored as part of a dataset and +c interpreted by other programs. +f A MathMap is a Mapping which allows you to specify a set of forward +f and/or inverse transformation functions using arithmetic operations +f and mathematical functions similar to those available in Fortran. The +f MathMap interprets these functions at run-time, whenever its forward +f or inverse transformation is required. Because the functions are not +f compiled in the normal sense (unlike an IntraMap), they may be used to +f describe coordinate transformations in a transportable manner. A +f MathMap therefore provides a flexible way of defining new types of +f Mapping whose descriptions may be stored as part of a dataset and +f interpreted by other programs. + +* Parameters: +c nin +f NIN = INTEGER +* Number of input variables for the MathMap. This determines the +* value of its Nin attribute. +c nout +f NOUT = INTEGER +* Number of output variables for the MathMap. This determines the +* value of its Nout attribute. +c nfwd +f NFWD = INTEGER +* The number of forward transformation functions being supplied. +c This must be at least equal to "nout", but may be increased to +f This must be at least equal to NOUT, but may be increased to +* accommodate any additional expressions which define intermediate +* variables for the forward transformation (see the "Calculating +* Intermediate Values" section below). +c fwd +f FWD = CHARACTER * ( * )( NFWD ) +c An array (with "nfwd" elements) of pointers to null terminated strings +c which contain the expressions defining the forward transformation. +f An array which contains the expressions defining the forward +f transformation. +* The syntax of these expressions is described below. +c ninv +f NINV = INTEGER +* The number of inverse transformation functions being supplied. +c This must be at least equal to "nin", but may be increased to +f This must be at least equal to NIN, but may be increased to +* accommodate any additional expressions which define intermediate +* variables for the inverse transformation (see the "Calculating +* Intermediate Values" section below). +c inv +f INV = CHARACTER * ( * )( NINV ) +c An array (with "ninv" elements) of pointers to null terminated strings +c which contain the expressions defining the inverse transformation. +f An array which contains the expressions defining the inverse +f transformation. +* The syntax of these expressions is described below. +c options +f OPTIONS = CHARACTER * ( * ) (Given) +c Pointer to a null-terminated string containing an optional +c comma-separated list of attribute assignments to be used for +c initialising the new MathMap. The syntax used is identical to +c that for the astSet function and may include "printf" format +c specifiers identified by "%" symbols in the normal way. +c If no initialisation is required, a zero-length string may be +c supplied. +f A character string containing an optional comma-separated +f list of attribute assignments to be used for initialising the +f new MathMap. The syntax used is identical to that for the +f AST_SET routine. If no initialisation is required, a blank +f value may be supplied. +c ... +c If the "options" string contains "%" format specifiers, then +c an optional list of additional arguments may follow it in +c order to supply values to be substituted for these +c specifiers. The rules for supplying these are identical to +c those for the astSet function (and for the C "printf" +c function). +f STATUS = INTEGER (Given and Returned) +f The global status. + +* Returned Value: +c astMathMap() +f AST_MATHMAP = INTEGER +* A pointer to the new MathMap. + +* Defining Transformation Functions: +c A MathMap's transformation functions are supplied as a set of +c expressions in an array of character strings. Normally you would +c supply the same number of expressions for the forward transformation, +c via the "fwd" parameter, as there are output variables (given by the +c MathMap's Nout attribute). For instance, if Nout is 2 you might use: +c - "r = sqrt( x * x + y * y )" +c - "theta = atan2( y, x )" +c +c which defines a transformation from Cartesian to polar +c coordinates. Here, the variables that appear on the left of each +c expression ("r" and "theta") provide names for the output variables +c and those that appear on the right ("x" and "y") are references to +c input variables. +f A MathMap's transformation functions are supplied as a set of +f expressions in an array of character strings. Normally you would +f supply the same number of expressions for the forward transformation, +f via the FWD argument, as there are output variables (given by the +f MathMap's Nout attribute). For instance, if Nout is 2 you might use: +f - 'R = SQRT( X * X + Y * Y )' +f - 'THETA = ATAN2( Y, X )' +f +f which defines a transformation from Cartesian to polar +f coordinates. Here, the variables that appear on the left of each +f expression (R and THETA) provide names for the output variables and +f those that appear on the right (X and Y) are references to input +f variables. +* +c To complement this, you must also supply expressions for the inverse +c transformation via the "inv" parameter. In this case, the number of +c expressions given would normally match the number of MathMap input +c coordinates (given by the Nin attribute). If Nin is 2, you might use: +c - "x = r * cos( theta )" +c - "y = r * sin( theta )" +c +c which expresses the transformation from polar to Cartesian +c coordinates. Note that here the input variables ("x" and "y") are +c named on the left of each expression, and the output variables ("r" +c and "theta") are referenced on the right. +f To complement this, you must also supply expressions for the inverse +f transformation via the INV argument. In this case, the number of +f expressions given would normally match the number of MathMap input +f coordinates (given by the Nin attribute). If Nin is 2, you might use: +f - 'X = R * COS( THETA )' +f - 'Y = R * SIN( THETA )' +f +f which expresses the transformation from polar to Cartesian +f coordinates. Note that here the input variables (X and Y) are named on +f the left of each expression, and the output variables (R and THETA) +f are referenced on the right. +* +* Normally, you cannot refer to a variable on the right of an expression +* unless it is named on the left of an expression in the complementary +* set of functions. Therefore both sets of functions (forward and +* inverse) must be formulated using the same consistent set of variable +* names. This means that if you wish to leave one of the transformations +* undefined, you must supply dummy expressions which simply name each of +* the output (or input) variables. For example, you might use: +c - "x" +c - "y" +f - 'X' +f - 'Y' +* +* for the inverse transformation above, which serves to name the input +* variables but without defining an inverse transformation. + +* Calculating Intermediate Values: +c It is sometimes useful to calculate intermediate values and then to +c use these in the final expressions for the output (or input) +c variables. This may be done by supplying additional expressions for +c the forward (or inverse) transformation functions. For instance, the +c following array of five expressions describes 2-dimensional pin-cushion +c distortion: +c - "r = sqrt( xin * xin + yin * yin )" +c - "rout = r * ( 1 + 0.1 * r * r )" +c - "theta = atan2( yin, xin )" +c - "xout = rout * cos( theta )" +c - "yout = rout * sin( theta )" +f It is sometimes useful to calculate intermediate values and then to +f use these in the final expressions for the output (or input) +f variables. This may be done by supplying additional expressions for +f the forward (or inverse) transformation functions. For instance, the +f following array of five expressions describes 2-dimensional pin-cushion +f distortion: +f - 'R = SQRT( XIN * XIN + YIN * YIN )' +f - 'ROUT = R * ( 1 + 0.1 * R * R )' +f - 'THETA = ATAN2( YIN, XIN )', +f - 'XOUT = ROUT * COS( THETA )' +f - 'YOUT = ROUT * SIN( THETA )' +* +c Here, we first calculate three intermediate results ("r", "rout" +c and "theta") and then use these to calculate the final results ("xout" +c and "yout"). The MathMap knows that only the final two results +c constitute values for the output variables because its Nout attribute +c is set to 2. You may define as many intermediate variables in this +c way as you choose. Having defined a variable, you may then refer to it +c on the right of any subsequent expressions. +f Here, we first calculate three intermediate results (R, ROUT +f and THETA) and then use these to calculate the final results (XOUT +f and YOUT). The MathMap knows that only the final two results +f constitute values for the output variables because its Nout attribute +f is set to 2. You may define as many intermediate variables in this +f way as you choose. Having defined a variable, you may then refer to it +f on the right of any subsequent expressions. +* +c Note that when defining the inverse transformation you may only refer +c to the output variables "xout" and "yout". The intermediate variables +c "r", "rout" and "theta" (above) are private to the forward +c transformation and may not be referenced by the inverse +c transformation. The inverse transformation may, however, define its +c own private intermediate variables. +f Note that when defining the inverse transformation you may only refer +f to the output variables XOUT and YOUT. The intermediate variables R, +f ROUT and THETA (above) are private to the forward transformation and +f may not be referenced by the inverse transformation. The inverse +f transformation may, however, define its own private intermediate +f variables. + +* Expression Syntax: +c The expressions given for the forward and inverse transformations +c closely follow the syntax of the C programming language (with some +c extensions for compatibility with Fortran). They may contain +c references to variables and literal constants, together with +c arithmetic, boolean, relational and bitwise operators, and function +c invocations. A set of symbolic constants is also available. Each of +c these is described in detail below. Parentheses may be used to +c over-ride the normal order of evaluation. There is no built-in limit +c to the length of expressions and they are insensitive to case or the +c presence of additional white space. +f The expressions given for the forward and inverse transformations +f closely follow the syntax of Fortran (with some extensions for +f compatibility with the C language). They may contain references to +f variables and literal constants, together with arithmetic, logical, +f relational and bitwise operators, and function invocations. A set of +f symbolic constants is also available. Each of these is described in +f detail below. Parentheses may be used to over-ride the normal order of +f evaluation. There is no built-in limit to the length of expressions +f and they are insensitive to case or the presence of additional white +f space. + +* Variables: +* Variable names must begin with an alphabetic character and may contain +* only alphabetic characters, digits, and the underscore character +* "_". There is no built-in limit to the length of variable names. + +* Literal Constants: +c Literal constants, such as "0", "1", "0.007" or "2.505e-16" may appear +c in expressions, with the decimal point and exponent being optional (a +c "D" may also be used as an exponent character for compatibility with +c Fortran). A unary minus "-" may be used as a prefix. +f Literal constants, such as "0", "1", "0.007" or "2.505E-16" may appear +f in expressions, with the decimal point and exponent being optional (a +f "D" may also be used as an exponent character). A unary minus "-" may +f be used as a prefix. + +* Arithmetic Precision: +* All arithmetic is floating point, performed in double precision. + +* Propagation of Missing Data: +* Unless indicated otherwise, if any argument of a function or operator +* has the value AST__BAD (indicating missing data), then the result of +* that function or operation is also AST__BAD, so that such values are +* propagated automatically through all operations performed by MathMap +* transformations. The special value AST__BAD can be represented in +* expressions by the symbolic constant "<bad>". +* +* A <bad> result (i.e. equal to AST__BAD) is also produced in response +* to any numerical error (such as division by zero or numerical +* overflow), or if an invalid argument value is provided to a function +* or operator. + +* Arithmetic Operators: +* The following arithmetic operators are available: +c - x1 + x2: Sum of "x1" and "x2". +f - X1 + X2: Sum of X1 and X2. +c - x1 - x2: Difference of "x1" and "x2". +f - X1 - X2: Difference of X1 and X2. +c - x1 * x2: Product of "x1" and "x1". +f - X1 * X2: Product of X1 and X2. +c - x1 / x2: Ratio of "x1" and "x2". +f - X1 / X2: Ratio of X1 and X2. +c - x1 ** x2: "x1" raised to the power of "x2". +f - X1 ** X2: X1 raised to the power of X2. +c - + x: Unary plus, has no effect on its argument. +f - + X: Unary plus, has no effect on its argument. +c - - x: Unary minus, negates its argument. +f - - X: Unary minus, negates its argument. + +c Boolean Operators: +f Logical Operators: +c Boolean values are represented using zero to indicate false and +c non-zero to indicate true. In addition, the value AST__BAD is taken to +c mean "unknown". The values returned by boolean operators may therefore +c be 0, 1 or AST__BAD. Where appropriate, "tri-state" logic is +c implemented. For example, "a||b" may evaluate to 1 if "a" is non-zero, +c even if "b" has the value AST__BAD. This is because the result of the +c operation would not be affected by the value of "b", so long as "a" is +c non-zero. +f Logical values are represented using zero to indicate .FALSE. and +f non-zero to indicate .TRUE.. In addition, the value AST__BAD is taken to +f mean "unknown". The values returned by logical operators may therefore +f be 0, 1 or AST__BAD. Where appropriate, "tri-state" logic is +f implemented. For example, A.OR.B may evaluate to 1 if A is non-zero, +f even if B has the value AST__BAD. This is because the result of the +f operation would not be affected by the value of B, so long as A is +f non-zero. +* +c The following boolean operators are available: +f The following logical operators are available: +c - x1 && x2: Boolean AND between "x1" and "x2", returning 1 if both "x1" +c and "x2" are non-zero, and 0 otherwise. This operator implements +c tri-state logic. (The synonym ".and." is also provided for compatibility +c with Fortran.) +f - X1 .AND. X2: Logical AND between X1 and X2, returning 1 if both X1 +f and X2 are non-zero, and 0 otherwise. This operator implements +f tri-state logic. (The synonym "&&" is also provided for compatibility +f with C.) +c - x1 || x2: Boolean OR between "x1" and "x2", returning 1 if either "x1" +c or "x2" are non-zero, and 0 otherwise. This operator implements +c tri-state logic. (The synonym ".or." is also provided for compatibility +c with Fortran.) +f - X1 .OR. X2: Logical OR between X1 and X2, returning 1 if either X1 +f or X2 are non-zero, and 0 otherwise. This operator implements +f tri-state logic. (The synonym "||" is also provided for compatibility +f with C.) +c - x1 ^^ x2: Boolean exclusive OR (XOR) between "x1" and "x2", returning +c 1 if exactly one of "x1" and "x2" is non-zero, and 0 otherwise. Tri-state +c logic is not used with this operator. (The synonyms ".neqv." and ".xor." +c are also provided for compatibility with Fortran, although the second +c of these is not standard.) +f - X1 .NEQV. X2: Logical exclusive OR (XOR) between X1 and X2, +f returning 1 if exactly one of X1 and X2 is non-zero, and 0 +f otherwise. Tri-state logic is not used with this operator. (The +f synonym ".XOR." is also provided, although this is not standard +f Fortran. In addition, the C-like synonym "^^" may be used, although +f this is also not standard.) +c - x1 .eqv. x2: This is provided only for compatibility with Fortran +c and tests whether the boolean states of "x1" and "x2" (i.e. true/false) +c are equal. It is the negative of the exclusive OR (XOR) function. +c Tri-state logic is not used with this operator. +f - X1 .EQV. X2: Tests whether the logical states of X1 and X2 +f (i.e. .TRUE./.FALSE.) are equal. It is the negative of the exclusive OR +f (XOR) function. Tri-state logic is not used with this operator. +c - ! x: Boolean unary NOT operation, returning 1 if "x" is zero, and +c 0 otherwise. (The synonym ".not." is also provided for compatibility +c with Fortran.) +f - .NOT. X: Logical unary NOT operation, returning 1 if X is zero, and +f 0 otherwise. (The synonym "!" is also provided for compatibility with +f C.) + +* Relational Operators: +c Relational operators return the boolean result (0 or 1) of comparing +c the values of two floating point values for equality or inequality. The +c value AST__BAD may also be returned if either argument is <bad>. +f Relational operators return the logical result (0 or 1) of comparing +f the values of two floating point values for equality or inequality. The +f value AST__BAD may also be returned if either argument is <bad>. +* +* The following relational operators are available: +c - x1 == x2: Tests whether "x1" equals "x1". (The synonym ".eq." is +c also provided for compatibility with Fortran.) +f - X1 .EQ. X2: Tests whether X1 equals X2. (The synonym "==" is also +f provided for compatibility with C.) +c - x1 != x2: Tests whether "x1" is unequal to "x2". (The synonym ".ne." +c is also provided for compatibility with Fortran.) +f - X1 .NE. X2: Tests whether X1 is unequal to X2. (The synonym "!=" is +f also provided for compatibility with C.) +c - x1 > x2: Tests whether "x1" is greater than "x2". (The synonym +c ".gt." is also provided for compatibility with Fortran.) +f - X1 .GT. X2: Tests whether X1 is greater than X2. (The synonym ">" is +f also provided for compatibility with C.) +c - x1 >= x2: Tests whether "x1" is greater than or equal to "x2". (The +c synonym ".ge." is also provided for compatibility with Fortran.) +f - X1 .GE. X2: Tests whether X1 is greater than or equal to X2. (The +f synonym ">=" is also provided for compatibility with C.) +c - x1 < x2: Tests whether "x1" is less than "x2". (The synonym ".lt." +c is also provided for compatibility with Fortran.) +f - X1 .LT. X2: Tests whether X1 is less than X2. (The synonym "<" is also +f provided for compatibility with C.) +c - x1 <= x2: Tests whether "x1" is less than or equal to "x2". (The +c synonym ".le." is also provided for compatibility with Fortran.) +f - X1 .LE. X2: Tests whether X1 is less than or equal to X2. (The synonym +f "<=" is also provided for compatibility with C.) +* +c Note that relational operators cannot usefully be used to compare +c values with the <bad> value (representing missing data), because the +c result is always <bad>. The isbad() function should be used instead. +f Note that relational operators cannot usefully be used to compare +f values with the <bad> value (representing missing data), because the +f result is always <bad>. The ISBAD() function should be used instead. +f +f Note, also, that because logical operators can operate on floating +f point values, care must be taken to use parentheses in some cases +f where they would not normally be required in Fortran. For example, +f the expresssion: +f - .NOT. A .EQ. B +f +f must be written: +f - .NOT. ( A .EQ. B ) +f +f to prevent the .NOT. operator from associating with the variable A. + +* Bitwise Operators: +c The bitwise operators provided by C are often useful when operating on +c raw data (e.g. from instruments), so they are also provided for use in +c MathMap expressions. In this case, however, the values on which they +c operate are floating point values rather than pure integers. In order +c to produce results which match the pure integer case, the operands are +c regarded as fixed point binary numbers (i.e. with the binary +c equivalent of a decimal point) with negative numbers represented using +c twos-complement notation. For integer values, the resulting bit +c pattern corresponds to that of the equivalent signed integer (digits +c to the right of the point being zero). Operations on the bits +c representing the fractional part are also possible, however. +f Bitwise operators are often useful when operating on raw data +f (e.g. from instruments), so they are provided for use in MathMap +f expressions. In this case, however, the values on which they operate +f are floating point values rather than the more usual pure integers. In +f order to produce results which match the pure integer case, the +f operands are regarded as fixed point binary numbers (i.e. with the +f binary equivalent of a decimal point) with negative numbers +f represented using twos-complement notation. For integer values, the +f resulting bit pattern corresponds to that of the equivalent signed +f integer (digits to the right of the point being zero). Operations on +f the bits representing the fractional part are also possible, however. +* +* The following bitwise operators are available: +c - x1 >> x2: Rightward bit shift. The integer value of "x2" is taken +c (rounding towards zero) and the bits representing "x1" are then +c shifted this number of places to the right (or to the left if the +c number of places is negative). This is equivalent to dividing "x1" by +c the corresponding power of 2. +f - X1 >> X2: Rightward bit shift. The integer value of X2 is taken +f (rounding towards zero) and the bits representing X1 are then +f shifted this number of places to the right (or to the left if the +f number of places is negative). This is equivalent to dividing X1 by +f the corresponding power of 2. +c - x1 << x2: Leftward bit shift. The integer value of "x2" is taken +c (rounding towards zero), and the bits representing "x1" are then +c shifted this number of places to the left (or to the right if the +c number of places is negative). This is equivalent to multiplying "x1" +c by the corresponding power of 2. +f - X1 << X2: Leftward bit shift. The integer value of X2 is taken +f (rounding towards zero), and the bits representing X1 are then +f shifted this number of places to the left (or to the right if the +f number of places is negative). This is equivalent to multiplying X1 +f by the corresponding power of 2. +c - x1 & x2: Bitwise AND between the bits of "x1" and those of "x2" +c (equivalent to a boolean AND applied at each bit position in turn). +f - X1 & X2: Bitwise AND between the bits of X1 and those of X2 +f (equivalent to a logical AND applied at each bit position in turn). +c - x1 | x2: Bitwise OR between the bits of "x1" and those of "x2" +c (equivalent to a boolean OR applied at each bit position in turn). +f - X1 | X2: Bitwise OR between the bits of X1 and those of X2 +f (equivalent to a logical OR applied at each bit position in turn). +c - x1 ^ x2: Bitwise exclusive OR (XOR) between the bits of "x1" and +c those of "x2" (equivalent to a boolean XOR applied at each bit +c position in turn). +f - X1 ^ X2: Bitwise exclusive OR (XOR) between the bits of X1 and +f those of X2 (equivalent to a logical XOR applied at each bit +f position in turn). +* +c Note that no bit inversion operator ("~" in C) is provided. This is +c because inverting the bits of a twos-complement fixed point binary +c number is equivalent to simply negating it. This differs from the +c pure integer case because bits to the right of the binary point are +c also inverted. To invert only those bits to the left of the binary +c point, use a bitwise exclusive OR with the value -1 (i.e. "x^-1"). +f Note that no bit inversion operator is provided. This is +f because inverting the bits of a twos-complement fixed point binary +f number is equivalent to simply negating it. This differs from the +f pure integer case because bits to the right of the binary point are +f also inverted. To invert only those bits to the left of the binary +f point, use a bitwise exclusive OR with the value -1 (i.e. X^-1). + +* Functions: +* The following functions are available: +c - abs(x): Absolute value of "x" (sign removal), same as fabs(x). +f - ABS(X): Absolute value of X (sign removal), same as FABS(X). +c - acos(x): Inverse cosine of "x", in radians. +f - ACOS(X): Inverse cosine of X, in radians. +c - acosd(x): Inverse cosine of "x", in degrees. +f - ACOSD(X): Inverse cosine of X, in degrees. +c - acosh(x): Inverse hyperbolic cosine of "x". +f - ACOSH(X): Inverse hyperbolic cosine of X. +c - acoth(x): Inverse hyperbolic cotangent of "x". +f - ACOTH(X): Inverse hyperbolic cotangent of X. +c - acsch(x): Inverse hyperbolic cosecant of "x". +f - ACSCH(X): Inverse hyperbolic cosecant of X. +c - aint(x): Integer part of "x" (round towards zero), same as int(x). +f - AINT(X): Integer part of X (round towards zero), same as INT(X). +c - asech(x): Inverse hyperbolic secant of "x". +f - ASECH(X): Inverse hyperbolic secant of X. +c - asin(x): Inverse sine of "x", in radians. +f - ASIN(X): Inverse sine of X, in radians. +c - asind(x): Inverse sine of "x", in degrees. +f - ASIND(X): Inverse sine of X, in degrees. +c - asinh(x): Inverse hyperbolic sine of "x". +f - ASINH(X): Inverse hyperbolic sine of X. +c - atan(x): Inverse tangent of "x", in radians. +f - ATAN(X): Inverse tangent of X, in radians. +c - atand(x): Inverse tangent of "x", in degrees. +f - ATAND(X): Inverse tangent of X, in degrees. +c - atanh(x): Inverse hyperbolic tangent of "x". +f - ATANH(X): Inverse hyperbolic tangent of X. +c - atan2(x1, x2): Inverse tangent of "x1/x2", in radians. +f - ATAN2(X1, X2): Inverse tangent of X1/X2, in radians. +c - atan2d(x1, x2): Inverse tangent of "x1/x2", in degrees. +f - ATAN2D(X1, X2): Inverse tangent of X1/X2, in degrees. +c - ceil(x): Smallest integer value not less then "x" (round towards +c plus infinity). +f - CEIL(X): Smallest integer value not less then X (round towards +f plus infinity). +c - cos(x): Cosine of "x" in radians. +f - COS(X): Cosine of X in radians. +c - cosd(x): Cosine of "x" in degrees. +f - COSD(X): Cosine of X in degrees. +c - cosh(x): Hyperbolic cosine of "x". +f - COSH(X): Hyperbolic cosine of X. +c - coth(x): Hyperbolic cotangent of "x". +f - COTH(X): Hyperbolic cotangent of X. +c - csch(x): Hyperbolic cosecant of "x". +f - CSCH(X): Hyperbolic cosecant of X. +c - dim(x1, x2): Returns "x1-x2" if "x1" is greater than "x2", otherwise 0. +f - DIM(X1, X2): Returns X1-X2 if X1 is greater than X2, otherwise 0. +c - exp(x): Exponential function of "x". +f - EXP(X): Exponential function of X. +c - fabs(x): Absolute value of "x" (sign removal), same as abs(x). +f - FABS(X): Absolute value of X (sign removal), same as ABS(X). +c - floor(x): Largest integer not greater than "x" (round towards +c minus infinity). +f - FLOOR(X): Largest integer not greater than X (round towards +f minus infinity). +c - fmod(x1, x2): Remainder when "x1" is divided by "x2", same as +c mod(x1, x2). +f - FMOD(X1, X2): Remainder when X1 is divided by X2, same as +f MOD(X1, X2). +c - gauss(x1, x2): Random sample from a Gaussian distribution with mean +c "x1" and standard deviation "x2". +f - GAUSS(X1, X2): Random sample from a Gaussian distribution with mean +f X1 and standard deviation X2. +c - int(x): Integer part of "x" (round towards zero), same as aint(x). +f - INT(X): Integer part of X (round towards zero), same as AINT(X). +c - isbad(x): Returns 1 if "x" has the <bad> value (AST__BAD), otherwise 0. +f - ISBAD(X): Returns 1 if X has the <bad> value (AST__BAD), otherwise 0. +c - log(x): Natural logarithm of "x". +f - LOG(X): Natural logarithm of X. +c - log10(x): Logarithm of "x" to base 10. +f - LOG10(X): Logarithm of X to base 10. +c - max(x1, x2, ...): Maximum of two or more values. +f - MAX(X1, X2, ...): Maximum of two or more values. +c - min(x1, x2, ...): Minimum of two or more values. +f - MIN(X1, X2, ...): Minimum of two or more values. +c - mod(x1, x2): Remainder when "x1" is divided by "x2", same as +c fmod(x1, x2). +f - MOD(X1, X2): Remainder when X1 is divided by X2, same as +f FMOD(X1, X2). +c - nint(x): Nearest integer to "x" (round to nearest). +f - NINT(X): Nearest integer to X (round to nearest). +c - poisson(x): Random integer-valued sample from a Poisson +c distribution with mean "x". +f - POISSON(X): Random integer-valued sample from a Poisson +f distribution with mean X. +c - pow(x1, x2): "x1" raised to the power of "x2". +f - POW(X1, X2): X1 raised to the power of X2. +c - qif(x1, x2, x3): Returns "x2" if "x1" is true, and "x3" otherwise. +f - QIF(x1, x2, x3): Returns X2 if X1 is true, and X3 otherwise. +c - rand(x1, x2): Random sample from a uniform distribution in the +c range "x1" to "x2" inclusive. +f - RAND(X1, X2): Random sample from a uniform distribution in the +f range X1 to X2 inclusive. +c - sech(x): Hyperbolic secant of "x". +f - SECH(X): Hyperbolic secant of X. +c - sign(x1, x2): Absolute value of "x1" with the sign of "x2" +c (transfer of sign). +f - SIGN(X1, X2): Absolute value of X1 with the sign of X2 +f (transfer of sign). +c - sin(x): Sine of "x" in radians. +f - SIN(X): Sine of X in radians. +c - sinc(x): Sinc function of "x" [= "sin(x)/x"]. +f - SINC(X): Sinc function of X [= SIN(X)/X]. +c - sind(x): Sine of "x" in degrees. +f - SIND(X): Sine of X in degrees. +c - sinh(x): Hyperbolic sine of "x". +f - SINH(X): Hyperbolic sine of X. +c - sqr(x): Square of "x" (= "x*x"). +f - SQR(X): Square of X (= X*X). +c - sqrt(x): Square root of "x". +f - SQRT(X): Square root of X. +c - tan(x): Tangent of "x" in radians. +f - TAN(X): Tangent of X in radians. +c - tand(x): Tangent of "x" in degrees. +f - TAND(X): Tangent of X in degrees. +c - tanh(x): Hyperbolic tangent of "x". +f - TANH(X): Hyperbolic tangent of X. + +* Symbolic Constants: +* The following symbolic constants are available (the enclosing "<>" +* brackets must be included): +c - <bad>: The "bad" value (AST__BAD) used to flag missing data. Note +c that you cannot usefully compare values with this constant because the +c result is always <bad>. The isbad() function should be used instead. +f - <bad>: The "bad" value (AST__BAD) used to flag missing data. Note +f that you cannot usefully compare values with this constant because the +f result is always <bad>. The ISBAD() function should be used instead. +c - <dig>: Number of decimal digits of precision available in a +c floating point (double) value. +f - <dig>: Number of decimal digits of precision available in a +f floating point (double precision) value. +* - <e>: Base of natural logarithms. +* - <epsilon>: Smallest positive number such that 1.0+<epsilon> is +* distinguishable from unity. +c - <mant_dig>: The number of base <radix> digits stored in the +c mantissa of a floating point (double) value. +f - <mant_dig>: The number of base <radix> digits stored in the +f mantissa of a floating point (double precision) value. +c - <max>: Maximum representable floating point (double) value. +f - <max>: Maximum representable floating point (double precision) value. +c - <max_10_exp>: Maximum integer such that 10 raised to that power +c can be represented as a floating point (double) value. +f - <max_10_exp>: Maximum integer such that 10 raised to that power +f can be represented as a floating point (double precision) value. +c - <max_exp>: Maximum integer such that <radix> raised to that +c power minus 1 can be represented as a floating point (double) value. +f - <max_exp>: Maximum integer such that <radix> raised to that +f power minus 1 can be represented as a floating point (double precision) +f value. +c - <min>: Smallest positive number which can be represented as a +c normalised floating point (double) value. +f - <min>: Smallest positive number which can be represented as a +f normalised floating point (double precision) value. +c - <min_10_exp>: Minimum negative integer such that 10 raised to that +c power can be represented as a normalised floating point (double) value. +f - <min_10_exp>: Minimum negative integer such that 10 raised to that +f power can be represented as a normalised floating point (double +f precision) value. +c - <min_exp>: Minimum negative integer such that <radix> raised to +c that power minus 1 can be represented as a normalised floating point +c (double) value. +f - <min_exp>: Minimum negative integer such that <radix> raised to +f that power minus 1 can be represented as a normalised floating point +f (double precision) value. +* - <pi>: Ratio of the circumference of a circle to its diameter. +c - <radix>: The radix (number base) used to represent the mantissa of +c floating point (double) values. +f - <radix>: The radix (number base) used to represent the mantissa of +f floating point (double precision) values. +* - <rounds>: The mode used for rounding floating point results after +* addition. Possible values include: -1 (indeterminate), 0 (toward +* zero), 1 (to nearest), 2 (toward plus infinity) and 3 (toward minus +* infinity). Other values indicate machine-dependent behaviour. + +* Evaluation Precedence and Associativity: +* Items appearing in expressions are evaluated in the following order +* (highest precedence first): +* - Constants and variables +* - Function arguments and parenthesised expressions +* - Function invocations +* - Unary + - ! .not. +* - ** +* - * / +* - + - +* - << >> +* - < .lt. <= .le. > .gt. >= .ge. +* - == .eq. != .ne. +* - & +* - ^ +* - | +* - && .and. +* - ^^ +* - || .or +* - .eqv. .neqv. .xor. +* +* All operators associate from left-to-right, except for unary +, +* unary -, !, .not. and ** which associate from right-to-left. + +* Notes: +* - The sequence of numbers produced by the random number functions +* available within a MathMap is normally unpredictable and different for +* each MathMap. However, this behaviour may be controlled by means of +* the MathMap's Seed attribute. +c - Normally, compound Mappings (CmpMaps) which involve MathMaps will +c not be subject to simplification (e.g. using astSimplify) because AST +c cannot know how different MathMaps will interact. However, in the +c special case where a MathMap occurs in series with its own inverse, +c then simplification may be possible. Whether simplification does, in +c fact, occur under these circumstances is controlled by the MathMap's +c SimpFI and SimpIF attributes. +f - Normally, compound Mappings (CmpMaps) which involve MathMaps will +f not be subject to simplification (e.g. using AST_SIMPLIFY) because AST +f cannot know how different MathMaps will interact. However, in the +f special case where a MathMap occurs in series with its own inverse, +f then simplification may be possible. Whether simplification does, in +f fact, occur under these circumstances is controlled by the MathMap's +f SimpFI and SimpIF attributes. +* - A null Object pointer (AST__NULL) will be returned if this +c function is invoked with the AST error status set, or if it +f function is invoked with STATUS set to an error value, or if it +* should fail for any reason. +*-- + +* Implementation Notes: +* - This function implements the external (public) interface to +* the astMathMap constructor function. It returns an ID value +* (instead of a true C pointer) to external users, and must be +* provided because astMathMap_ has a variable argument list which +* cannot be encapsulated in a macro (where this conversion would +* otherwise occur). +* - The variable argument list also prevents this function from +* invoking astMathMap_ directly, so it must be a re-implementation +* of it in all respects, except for the final conversion of the +* result to an ID value. +*/ + +/* Local Variables: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ + AstMathMap *new; /* Pointer to new MathMap */ + va_list args; /* Variable argument list */ + + int *status; /* Pointer to inherited status value */ + +/* Get a pointer to the inherited status value. */ + status = astGetStatusPtr; + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(NULL); + +/* Check the global error status. */ + if ( !astOK ) return NULL; + +/* Initialise the MathMap, allocating memory and initialising the virtual + function table as well if necessary. */ + new = astInitMathMap( NULL, sizeof( AstMathMap ), !class_init, &class_vtab, + "MathMap", nin, nout, nfwd, fwd, ninv, inv ); + +/* If successful, note that the virtual function table has been initialised. */ + if ( astOK ) { + class_init = 1; + +/* Obtain the variable argument list and pass it along with the options string + to the astVSet method to initialise the new MathMap's attributes. */ + va_start( args, options ); + astVSet( new, options, NULL, args ); + va_end( args ); + +/* If an error occurred, clean up by deleting the new object. */ + if ( !astOK ) new = astDelete( new ); + } + +/* Return an ID value for the new MathMap. */ + return astMakeId( new ); +} + +AstMathMap *astInitMathMap_( void *mem, size_t size, int init, + AstMathMapVtab *vtab, const char *name, + int nin, int nout, + int nfwd, const char *fwd[], + int ninv, const char *inv[], int *status ) { +/* +*+ +* Name: +* astInitMathMap + +* Purpose: +* Initialise a MathMap. + +* Type: +* Protected function. + +* Synopsis: +* #include "mathmap.h" +* AstMathMap *astInitMathMap_( void *mem, size_t size, int init, +* AstMathMapVtab *vtab, const char *name, +* int nin, int nout, +* int nfwd, const char *fwd[], +* int ninv, const char *inv[] ) + +* Class Membership: +* MathMap initialiser. + +* Description: +* This function is provided for use by class implementations to initialise +* a new MathMap object. It allocates memory (if necessary) to accommodate +* the MathMap plus any additional data associated with the derived class. +* It then initialises a MathMap structure at the start of this memory. If +* the "init" flag is set, it also initialises the contents of a virtual +* function table for a MathMap at the start of the memory passed via the +* "vtab" parameter. + +* Parameters: +* mem +* A pointer to the memory in which the MathMap is to be initialised. +* This must be of sufficient size to accommodate the MathMap data +* (sizeof(MathMap)) plus any data used by the derived class. If a value +* of NULL is given, this function will allocate the memory itself using +* the "size" parameter to determine its size. +* size +* The amount of memory used by the MathMap (plus derived class data). +* This will be used to allocate memory if a value of NULL is given for +* the "mem" parameter. This value is also stored in the MathMap +* structure, so a valid value must be supplied even if not required for +* allocating memory. +* init +* A logical flag indicating if the MathMap's virtual function table is +* to be initialised. If this value is non-zero, the virtual function +* table will be initialised by this function. +* vtab +* Pointer to the start of the virtual function table to be associated +* with the new MathMap. +* name +* Pointer to a constant null-terminated character string which contains +* the name of the class to which the new object belongs (it is this +* pointer value that will subsequently be returned by the Object +* astClass function). +* nin +* Number of input variables for the MathMap. +* nout +* Number of output variables for the MathMap. +* nfwd +* The number of forward transformation functions being supplied. +* This must be at least equal to "nout". +* fwd +* Pointer to an array, with "nfwd" elements, of pointers to null +* terminated strings which contain each of the forward transformation +* functions. +* ninv +* The number of inverse transformation functions being supplied. +* This must be at least equal to "nin". +* inv +* Pointer to an array, with "ninv" elements, of pointers to null +* terminated strings which contain each of the inverse transformation +* functions. + +* Returned Value: +* A pointer to the new MathMap. + +* Notes: +* - This function does not attempt to ensure that the forward and inverse +* transformations performed by the resulting MathMap are consistent in any +* way. +* - This function makes a copy of the contents of the strings supplied. +* - A null pointer will be returned if this function is invoked with the +* global error status set, or if it should fail for any reason. +*- +*/ + +/* Local Variables: */ + AstMathMap *new; /* Pointer to new MathMap */ + char **fwdfun; /* Array of cleaned forward functions */ + char **invfun; /* Array of cleaned inverse functions */ + double **fwdcon; /* Constants for forward functions */ + double **invcon; /* Constants for inverse functions */ + int **fwdcode; /* Code for forward functions */ + int **invcode; /* Code for inverse functions */ + int fwdstack; /* Stack size for forward functions */ + int invstack; /* Stack size for inverse functions */ + +/* Initialise. */ + new = NULL; + +/* Check the global status. */ + if ( !astOK ) return new; + +/* If necessary, initialise the virtual function table. */ + if ( init ) astInitMathMapVtab( vtab, name ); + +/* Check the numbers of input and output variables for validity, + reporting an error if necessary. */ + if ( nin < 1 ) { + astError( AST__BADNI, + "astInitMathMap(%s): Bad number of input coordinates (%d).", status, + name, nin ); + astError( AST__BADNI, + "This number should be one or more." , status); + } else if ( nout < 1 ) { + astError( AST__BADNO, + "astInitMathMap(%s): Bad number of output coordinates (%d).", status, + name, nout ); + astError( AST__BADNI, + "This number should be one or more." , status); + +/* Check that sufficient number of forward and inverse transformation + functions have been supplied and report an error if necessary. */ + } else if ( nfwd < nout ) { + astError( AST__INNTF, + "astInitMathMap(%s): Too few forward transformation functions " + "given (%d).", status, + name, nfwd ); + astError( astStatus, + "At least %d forward transformation functions must be " + "supplied. ", status, + nout ); + } else if ( ninv < nin ) { + astError( AST__INNTF, + "astInitMathMap(%s): Too few inverse transformation functions " + "given (%d).", status, + name, ninv ); + astError( astStatus, + "At least %d inverse transformation functions must be " + "supplied. ", status, + nin ); + +/* Of OK, clean the forward and inverse functions provided. This makes + a lower-case copy with white space removed. */ + } else { + CleanFunctions( nfwd, fwd, &fwdfun, status ); + CleanFunctions( ninv, inv, &invfun, status ); + +/* Compile the cleaned functions. From the returned pointers (if + successful), we can now tell which transformations (forward and/or + inverse) are defined. */ + CompileMapping( "astInitMathMap", name, nin, nout, + nfwd, (const char **) fwdfun, + ninv, (const char **) invfun, + &fwdcode, &invcode, &fwdcon, &invcon, + &fwdstack, &invstack, status ); + +/* Initialise a Mapping structure (the parent class) as the first + component within the MathMap structure, allocating memory if + necessary. Specify that the Mapping should be defined in the required + directions. */ + new = (AstMathMap *) astInitMapping( mem, size, 0, + (AstMappingVtab *) vtab, name, + nin, nout, + ( fwdcode != NULL ), + ( invcode != NULL ) ); + + +/* If an error has occurred, free all the memory which may have been + allocated by the cleaning and compilation steps above. */ + if ( !astOK ) { + FREE_POINTER_ARRAY( fwdfun, nfwd ) + FREE_POINTER_ARRAY( invfun, ninv ) + FREE_POINTER_ARRAY( fwdcode, nfwd ) + FREE_POINTER_ARRAY( invcode, ninv ) + FREE_POINTER_ARRAY( fwdcon, nfwd ) + FREE_POINTER_ARRAY( invcon, ninv ) + } + +/* Initialise the MathMap data. */ +/* ---------------------------- */ +/* Store pointers to the compiled function information, together with + other MathMap data. */ + if ( new ) { + new->fwdfun = fwdfun; + new->invfun = invfun; + new->fwdcode = fwdcode; + new->invcode = invcode; + new->fwdcon = fwdcon; + new->invcon = invcon; + new->fwdstack = fwdstack; + new->invstack = invstack; + new->nfwd = nfwd; + new->ninv = ninv; + new->simp_fi = -INT_MAX; + new->simp_if = -INT_MAX; + +/* Initialise the random number generator context associated with the + MathMap, using an unpredictable default seed value. */ + new->rcontext.active = 0; + new->rcontext.random_int = 0; + new->rcontext.seed_set = 0; + new->rcontext.seed = DefaultSeed( &new->rcontext, status ); + +/* If an error occurred, clean up by deleting the new object. */ + if ( !astOK ) new = astDelete( new ); + } + } + +/* Return a pointer to the new object. */ + return new; +} + +AstMathMap *astLoadMathMap_( void *mem, size_t size, + AstMathMapVtab *vtab, const char *name, + AstChannel *channel, int *status ) { +/* +*+ +* Name: +* astLoadMathMap + +* Purpose: +* Load a MathMap. + +* Type: +* Protected function. + +* Synopsis: +* #include "mathmap.h" +* AstMathMap *astLoadMathMap( void *mem, size_t size, +* AstMathMapVtab *vtab, const char *name, +* AstChannel *channel ) + +* Class Membership: +* MathMap loader. + +* Description: +* This function is provided to load a new MathMap using data read +* from a Channel. It first loads the data used by the parent class +* (which allocates memory if necessary) and then initialises a +* MathMap structure in this memory, using data read from the input +* Channel. +* +* If the "init" flag is set, it also initialises the contents of a +* virtual function table for a MathMap at the start of the memory +* passed via the "vtab" parameter. + + +* Parameters: +* mem +* A pointer to the memory into which the MathMap is to be +* loaded. This must be of sufficient size to accommodate the +* MathMap data (sizeof(MathMap)) plus any data used by derived +* classes. If a value of NULL is given, this function will +* allocate the memory itself using the "size" parameter to +* determine its size. +* size +* The amount of memory used by the MathMap (plus derived class +* data). This will be used to allocate memory if a value of +* NULL is given for the "mem" parameter. This value is also +* stored in the MathMap structure, so a valid value must be +* supplied even if not required for allocating memory. +* +* If the "vtab" parameter is NULL, the "size" value is ignored +* and sizeof(AstMathMap) is used instead. +* vtab +* Pointer to the start of the virtual function table to be +* associated with the new MathMap. If this is NULL, a pointer +* to the (static) virtual function table for the MathMap class +* is used instead. +* name +* Pointer to a constant null-terminated character string which +* contains the name of the class to which the new object +* belongs (it is this pointer value that will subsequently be +* returned by the astGetClass method). +* +* If the "vtab" parameter is NULL, the "name" value is ignored +* and a pointer to the string "MathMap" is used instead. + +* Returned Value: +* A pointer to the new MathMap. + +* Notes: +* - A null pointer will be returned if this function is invoked +* with the global error status set, or if it should fail for any +* reason. +*- +*/ + +/* Local Constants: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ +#define KEY_LEN 50 /* Maximum length of a keyword */ + +/* Local Variables: */ + AstMathMap *new; /* Pointer to the new MathMap */ + char key[ KEY_LEN + 1 ]; /* Buffer for keyword strings */ + int ifun; /* Loop counter for functions */ + int invert; /* Invert attribute value */ + int nin; /* True number of input coordinates */ + int nout; /* True number of output coordinates */ + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(channel); + +/* Initialise. */ + new = NULL; + +/* Check the global error status. */ + if ( !astOK ) return new; + +/* If a NULL virtual function table has been supplied, then this is + the first loader to be invoked for this MathMap. In this case the + MathMap belongs to this class, so supply appropriate values to be + passed to the parent class loader (and its parent, etc.). */ + if ( !vtab ) { + size = sizeof( AstMathMap ); + vtab = &class_vtab; + name = "MathMap"; + +/* If required, initialise the virtual function table for this class. */ + if ( !class_init ) { + astInitMathMapVtab( vtab, name ); + class_init = 1; + } + } + +/* Invoke the parent class loader to load data for all the ancestral + classes of the current one, returning a pointer to the resulting + partly-built MathMap. */ + new = astLoadMapping( mem, size, (AstMappingVtab *) vtab, name, + channel ); + + if ( astOK ) { + +/* Read input data. */ +/* ================ */ +/* Request the input Channel to read all the input data appropriate to + this class into the internal "values list". */ + astReadClassData( channel, "MathMap" ); + +/* Determine if the MathMap is inverted and obtain the "true" number + of input and output coordinates by un-doing the effects of any + inversion. */ + invert = astGetInvert( new ); + nin = invert ? astGetNout( new ) : astGetNin( new ); + nout = invert ? astGetNin( new ) : astGetNout( new ); + +/* Now read each individual data item from this list and use it to + initialise the appropriate instance variable(s) for this class. */ + +/* In the case of attributes, we first read the "raw" input value, + supplying the "unset" value as the default. If a "set" value is + obtained, we then use the appropriate (private) Set... member + function to validate and set the value properly. */ + +/* Numbers of transformation functions. */ +/* ------------------------------------ */ +/* Read the numbers of forward and inverse transformation functions, + supplying appropriate defaults. */ + new->nfwd = astReadInt( channel, "nfwd", nout ); + new->ninv = astReadInt( channel, "ninv", nin ); + if ( astOK ) { + +/* Allocate memory for the MathMap's transformation function arrays. */ + MALLOC_POINTER_ARRAY( new->fwdfun, char *, new->nfwd ) + MALLOC_POINTER_ARRAY( new->invfun, char *, new->ninv ) + if ( astOK ) { + +/* Forward transformation functions. */ +/* --------------------------------- */ +/* Create a keyword for each forward transformation function and read + the function's value as a string. */ + for ( ifun = 0; ifun < new->nfwd; ifun++ ) { + (void) sprintf( key, "fwd%d", ifun + 1 ); + new->fwdfun[ ifun ] = astReadString( channel, key, "" ); + } + +/* Inverse transformation functions. */ +/* --------------------------------- */ +/* Repeat this process for the inverse transformation functions. */ + for ( ifun = 0; ifun < new->ninv; ifun++ ) { + (void) sprintf( key, "inv%d", ifun + 1 ); + new->invfun[ ifun ] = astReadString( channel, key, "" ); + } + +/* Forward-inverse simplification flag. */ +/* ------------------------------------ */ + new->simp_fi = astReadInt( channel, "simpfi", -INT_MAX ); + if ( TestSimpFI( new, status ) ) SetSimpFI( new, new->simp_fi, status ); + +/* Inverse-forward simplification flag. */ +/* ------------------------------------ */ + new->simp_if = astReadInt( channel, "simpif", -INT_MAX ); + if ( TestSimpIF( new, status ) ) SetSimpIF( new, new->simp_if, status ); + +/* Random number context. */ +/* ---------------------- */ +/* Initialise the random number generator context. */ + new->rcontext.active = 0; + new->rcontext.random_int = 0; + +/* Read the flag that determines if the Seed value is set, and the + Seed value itself. */ + new->rcontext.seed_set = astReadInt( channel, "seeded", 0 ); + if ( TestSeed( new, status ) ) { + new->rcontext.seed = astReadInt( channel, "seed", 0 ); + SetSeed( new, new->rcontext.seed, status ); + +/* Supply an unpredictable default Seed value if necessary. */ + } else { + new->rcontext.seed = DefaultSeed( &new->rcontext, status ); + } + +/* Compile the MathMap's transformation functions. */ + CompileMapping( "astLoadMathMap", name, nin, nout, + new->nfwd, (const char **) new->fwdfun, + new->ninv, (const char **) new->invfun, + &new->fwdcode, &new->invcode, + &new->fwdcon, &new->invcon, + &new->fwdstack, &new->invstack, status ); + } + +/* If an error occurred, clean up by deleting the new MathMap. */ + if ( !astOK ) new = astDelete( new ); + } + } + +/* Return the new MathMap pointer. */ + return new; + +/* Undefine macros local to this function. */ +#undef KEY_LEN +} + + + + + |