/cfitsio/eval.y
Happy | 5837 lines | 5275 code | 562 blank | 0 comment | 0 complexity | 05e094fffb8f04f7bfaeec5c1e661f58 MD5 | raw file
- %{
- /************************************************************************/
- /* */
- /* CFITSIO Lexical Parser */
- /* */
- /* This file is one of 3 files containing code which parses an */
- /* arithmetic expression and evaluates it in the context of an input */
- /* FITS file table extension. The CFITSIO lexical parser is divided */
- /* into the following 3 parts/files: the CFITSIO "front-end", */
- /* eval_f.c, contains the interface between the user/CFITSIO and the */
- /* real core of the parser; the FLEX interpreter, eval_l.c, takes the */
- /* input string and parses it into tokens and identifies the FITS */
- /* information required to evaluate the expression (ie, keywords and */
- /* columns); and, the BISON grammar and evaluation routines, eval_y.c, */
- /* receives the FLEX output and determines and performs the actual */
- /* operations. The files eval_l.c and eval_y.c are produced from */
- /* running flex and bison on the files eval.l and eval.y, respectively. */
- /* (flex and bison are available from any GNU archive: see www.gnu.org) */
- /* */
- /* The grammar rules, rather than evaluating the expression in situ, */
- /* builds a tree, or Nodal, structure mapping out the order of */
- /* operations and expression dependencies. This "compilation" process */
- /* allows for much faster processing of multiple rows. This technique */
- /* was developed by Uwe Lammers of the XMM Science Analysis System, */
- /* although the CFITSIO implementation is entirely code original. */
- /* */
- /* */
- /* Modification History: */
- /* */
- /* Kent Blackburn c1992 Original parser code developed for the */
- /* FTOOLS software package, in particular, */
- /* the fselect task. */
- /* Kent Blackburn c1995 BIT column support added */
- /* Peter D Wilson Feb 1998 Vector column support added */
- /* Peter D Wilson May 1998 Ported to CFITSIO library. User */
- /* interface routines written, in essence */
- /* making fselect, fcalc, and maketime */
- /* capabilities available to all tools */
- /* via single function calls. */
- /* Peter D Wilson Jun 1998 Major rewrite of parser core, so as to */
- /* create a run-time evaluation tree, */
- /* inspired by the work of Uwe Lammers, */
- /* resulting in a speed increase of */
- /* 10-100 times. */
- /* Peter D Wilson Jul 1998 gtifilter(a,b,c,d) function added */
- /* Peter D Wilson Aug 1998 regfilter(a,b,c,d) function added */
- /* Peter D Wilson Jul 1999 Make parser fitsfile-independent, */
- /* allowing a purely vector-based usage */
- /* Craig B Markwardt Jun 2004 Add MEDIAN() function */
- /* Craig B Markwardt Jun 2004 Add SUM(), and MIN/MAX() for bit arrays */
- /* Craig B Markwardt Jun 2004 Allow subscripting of nX bit arrays */
- /* Craig B Markwardt Jun 2004 Implement statistical functions */
- /* NVALID(), AVERAGE(), and STDDEV() */
- /* for integer and floating point vectors */
- /* Craig B Markwardt Jun 2004 Use NULL values for range errors instead*/
- /* of throwing a parse error */
- /* Craig B Markwardt Oct 2004 Add ACCUM() and SEQDIFF() functions */
- /* Craig B Markwardt Feb 2005 Add ANGSEP() function */
- /* Craig B Markwardt Aug 2005 CIRCLE, BOX, ELLIPSE, NEAR and REGFILTER*/
- /* functions now accept vector arguments */
- /* Craig B Markwardt Sum 2006 Add RANDOMN() and RANDOMP() functions */
- /* Craig B Markwardt Mar 2007 Allow arguments to RANDOM and RANDOMN to*/
- /* determine the output dimensions */
- /* Craig B Markwardt Aug 2009 Add substring STRMID() and string search*/
- /* STRSTR() functions; more overflow checks*/
- /* */
- /************************************************************************/
- #define APPROX 1.0e-7
- #include "eval_defs.h"
- #include "region.h"
- #include <time.h>
- #include <stdlib.h>
- #ifndef alloca
- #define alloca malloc
- #endif
- /* Shrink the initial stack depth to keep local data <32K (mac limit) */
- /* yacc will allocate more space if needed, though. */
- #define YYINITDEPTH 100
- /***************************************************************/
- /* Replace Bison's BACKUP macro with one that fixes a bug -- */
- /* must update state after popping the stack -- and allows */
- /* popping multiple terms at one time. */
- /***************************************************************/
- #define YYNEWBACKUP(token, value) \
- do \
- if (yychar == YYEMPTY ) \
- { yychar = (token); \
- memcpy( &yylval, &(value), sizeof(value) ); \
- yychar1 = YYTRANSLATE (yychar); \
- while (yylen--) YYPOPSTACK; \
- yystate = *yyssp; \
- goto yybackup; \
- } \
- else \
- { yyerror ("syntax error: cannot back up"); YYERROR; } \
- while (0)
- /***************************************************************/
- /* Useful macros for accessing/testing Nodes */
- /***************************************************************/
- #define TEST(a) if( (a)<0 ) YYERROR
- #define SIZE(a) gParse.Nodes[ a ].value.nelem
- #define TYPE(a) gParse.Nodes[ a ].type
- #define OPER(a) gParse.Nodes[ a ].operation
- #define PROMOTE(a,b) if( TYPE(a) > TYPE(b) ) \
- b = New_Unary( TYPE(a), 0, b ); \
- else if( TYPE(a) < TYPE(b) ) \
- a = New_Unary( TYPE(b), 0, a );
- /***** Internal functions *****/
- #ifdef __cplusplus
- extern "C" {
- #endif
- static int Alloc_Node ( void );
- static void Free_Last_Node( void );
- static void Evaluate_Node ( int thisNode );
- static int New_Const ( int returnType, void *value, long len );
- static int New_Column( int ColNum );
- static int New_Offset( int ColNum, int offset );
- static int New_Unary ( int returnType, int Op, int Node1 );
- static int New_BinOp ( int returnType, int Node1, int Op, int Node2 );
- static int New_Func ( int returnType, funcOp Op, int nNodes,
- int Node1, int Node2, int Node3, int Node4,
- int Node5, int Node6, int Node7 );
- static int New_FuncSize( int returnType, funcOp Op, int nNodes,
- int Node1, int Node2, int Node3, int Node4,
- int Node5, int Node6, int Node7, int Size);
- static int New_Deref ( int Var, int nDim,
- int Dim1, int Dim2, int Dim3, int Dim4, int Dim5 );
- static int New_GTI ( char *fname, int Node1, char *start, char *stop );
- static int New_REG ( char *fname, int NodeX, int NodeY, char *colNames );
- static int New_Vector( int subNode );
- static int Close_Vec ( int vecNode );
- static int Locate_Col( Node *this );
- static int Test_Dims ( int Node1, int Node2 );
- static void Copy_Dims ( int Node1, int Node2 );
- static void Allocate_Ptrs( Node *this );
- static void Do_Unary ( Node *this );
- static void Do_Offset ( Node *this );
- static void Do_BinOp_bit ( Node *this );
- static void Do_BinOp_str ( Node *this );
- static void Do_BinOp_log ( Node *this );
- static void Do_BinOp_lng ( Node *this );
- static void Do_BinOp_dbl ( Node *this );
- static void Do_Func ( Node *this );
- static void Do_Deref ( Node *this );
- static void Do_GTI ( Node *this );
- static void Do_REG ( Node *this );
- static void Do_Vector ( Node *this );
- static long Search_GTI ( double evtTime, long nGTI, double *start,
- double *stop, int ordered );
- static char saobox (double xcen, double ycen, double xwid, double ywid,
- double rot, double xcol, double ycol);
- static char ellipse(double xcen, double ycen, double xrad, double yrad,
- double rot, double xcol, double ycol);
- static char circle (double xcen, double ycen, double rad,
- double xcol, double ycol);
- static char bnear (double x, double y, double tolerance);
- static char bitcmp (char *bitstrm1, char *bitstrm2);
- static char bitlgte(char *bits1, int oper, char *bits2);
- static void bitand(char *result, char *bitstrm1, char *bitstrm2);
- static void bitor (char *result, char *bitstrm1, char *bitstrm2);
- static void bitnot(char *result, char *bits);
- static int cstrmid(char *dest_str, int dest_len,
- char *src_str, int src_len, int pos);
- static void yyerror(char *msg);
- #ifdef __cplusplus
- }
- #endif
- %}
- %union {
- int Node; /* Index of Node */
- double dbl; /* real value */
- long lng; /* integer value */
- char log; /* logical value */
- char str[MAX_STRLEN]; /* string value */
- }
- %token <log> BOOLEAN /* First 3 must be in order of */
- %token <lng> LONG /* increasing promotion for later use */
- %token <dbl> DOUBLE
- %token <str> STRING
- %token <str> BITSTR
- %token <str> FUNCTION
- %token <str> BFUNCTION /* Bit function */
- %token <str> IFUNCTION /* Integer function */
- %token <str> GTIFILTER
- %token <str> REGFILTER
- %token <lng> COLUMN
- %token <lng> BCOLUMN
- %token <lng> SCOLUMN
- %token <lng> BITCOL
- %token <lng> ROWREF
- %token <lng> NULLREF
- %token <lng> SNULLREF
- %type <Node> expr
- %type <Node> bexpr
- %type <Node> sexpr
- %type <Node> bits
- %type <Node> vector
- %type <Node> bvector
- %left ',' '=' ':' '{' '}'
- %right '?'
- %left OR
- %left AND
- %left EQ NE '~'
- %left GT LT LTE GTE
- %left '+' '-' '%'
- %left '*' '/'
- %left '|' '&'
- %right POWER
- %left NOT
- %left INTCAST FLTCAST
- %left UMINUS
- %left '['
- %right ACCUM DIFF
- %%
- lines: /* nothing ; was | lines line */
- | lines line
- ;
- line: '\n' {}
- | expr '\n'
- { if( $1<0 ) {
- yyerror("Couldn't build node structure: out of memory?");
- YYERROR; }
- gParse.resultNode = $1;
- }
- | bexpr '\n'
- { if( $1<0 ) {
- yyerror("Couldn't build node structure: out of memory?");
- YYERROR; }
- gParse.resultNode = $1;
- }
- | sexpr '\n'
- { if( $1<0 ) {
- yyerror("Couldn't build node structure: out of memory?");
- YYERROR; }
- gParse.resultNode = $1;
- }
- | bits '\n'
- { if( $1<0 ) {
- yyerror("Couldn't build node structure: out of memory?");
- YYERROR; }
- gParse.resultNode = $1;
- }
- | error '\n' { yyerrok; }
- ;
- bvector: '{' bexpr
- { $$ = New_Vector( $2 ); TEST($$); }
- | bvector ',' bexpr
- {
- if( gParse.Nodes[$1].nSubNodes >= MAXSUBS ) {
- $1 = Close_Vec( $1 ); TEST($1);
- $$ = New_Vector( $1 ); TEST($$);
- } else {
- $$ = $1;
- }
- gParse.Nodes[$$].SubNodes[ gParse.Nodes[$$].nSubNodes++ ]
- = $3;
- }
- ;
- vector: '{' expr
- { $$ = New_Vector( $2 ); TEST($$); }
- | vector ',' expr
- {
- if( TYPE($1) < TYPE($3) )
- TYPE($1) = TYPE($3);
- if( gParse.Nodes[$1].nSubNodes >= MAXSUBS ) {
- $1 = Close_Vec( $1 ); TEST($1);
- $$ = New_Vector( $1 ); TEST($$);
- } else {
- $$ = $1;
- }
- gParse.Nodes[$$].SubNodes[ gParse.Nodes[$$].nSubNodes++ ]
- = $3;
- }
- | vector ',' bexpr
- {
- if( gParse.Nodes[$1].nSubNodes >= MAXSUBS ) {
- $1 = Close_Vec( $1 ); TEST($1);
- $$ = New_Vector( $1 ); TEST($$);
- } else {
- $$ = $1;
- }
- gParse.Nodes[$$].SubNodes[ gParse.Nodes[$$].nSubNodes++ ]
- = $3;
- }
- | bvector ',' expr
- {
- TYPE($1) = TYPE($3);
- if( gParse.Nodes[$1].nSubNodes >= MAXSUBS ) {
- $1 = Close_Vec( $1 ); TEST($1);
- $$ = New_Vector( $1 ); TEST($$);
- } else {
- $$ = $1;
- }
- gParse.Nodes[$$].SubNodes[ gParse.Nodes[$$].nSubNodes++ ]
- = $3;
- }
- ;
- expr: vector '}'
- { $$ = Close_Vec( $1 ); TEST($$); }
- ;
- bexpr: bvector '}'
- { $$ = Close_Vec( $1 ); TEST($$); }
- ;
- bits: BITSTR
- {
- $$ = New_Const( BITSTR, $1, strlen($1)+1 ); TEST($$);
- SIZE($$) = strlen($1); }
- | BITCOL
- { $$ = New_Column( $1 ); TEST($$); }
- | BITCOL '{' expr '}'
- {
- if( TYPE($3) != LONG
- || OPER($3) != CONST_OP ) {
- yyerror("Offset argument must be a constant integer");
- YYERROR;
- }
- $$ = New_Offset( $1, $3 ); TEST($$);
- }
- | bits '&' bits
- { $$ = New_BinOp( BITSTR, $1, '&', $3 ); TEST($$);
- SIZE($$) = ( SIZE($1)>SIZE($3) ? SIZE($1) : SIZE($3) ); }
- | bits '|' bits
- { $$ = New_BinOp( BITSTR, $1, '|', $3 ); TEST($$);
- SIZE($$) = ( SIZE($1)>SIZE($3) ? SIZE($1) : SIZE($3) ); }
- | bits '+' bits
- {
- if (SIZE($1)+SIZE($3) >= MAX_STRLEN) {
- yyerror("Combined bit string size exceeds " MAX_STRLEN_S " bits");
- YYERROR;
- }
- $$ = New_BinOp( BITSTR, $1, '+', $3 ); TEST($$);
- SIZE($$) = SIZE($1) + SIZE($3);
- }
- | bits '[' expr ']'
- { $$ = New_Deref( $1, 1, $3, 0, 0, 0, 0 ); TEST($$); }
- | bits '[' expr ',' expr ']'
- { $$ = New_Deref( $1, 2, $3, $5, 0, 0, 0 ); TEST($$); }
- | bits '[' expr ',' expr ',' expr ']'
- { $$ = New_Deref( $1, 3, $3, $5, $7, 0, 0 ); TEST($$); }
- | bits '[' expr ',' expr ',' expr ',' expr ']'
- { $$ = New_Deref( $1, 4, $3, $5, $7, $9, 0 ); TEST($$); }
- | bits '[' expr ',' expr ',' expr ',' expr ',' expr ']'
- { $$ = New_Deref( $1, 5, $3, $5, $7, $9, $11 ); TEST($$); }
- | NOT bits
- { $$ = New_Unary( BITSTR, NOT, $2 ); TEST($$); }
- | '(' bits ')'
- { $$ = $2; }
- ;
- expr: LONG
- { $$ = New_Const( LONG, &($1), sizeof(long) ); TEST($$); }
- | DOUBLE
- { $$ = New_Const( DOUBLE, &($1), sizeof(double) ); TEST($$); }
- | COLUMN
- { $$ = New_Column( $1 ); TEST($$); }
- | COLUMN '{' expr '}'
- {
- if( TYPE($3) != LONG
- || OPER($3) != CONST_OP ) {
- yyerror("Offset argument must be a constant integer");
- YYERROR;
- }
- $$ = New_Offset( $1, $3 ); TEST($$);
- }
- | ROWREF
- { $$ = New_Func( LONG, row_fct, 0, 0, 0, 0, 0, 0, 0, 0 ); }
- | NULLREF
- { $$ = New_Func( LONG, null_fct, 0, 0, 0, 0, 0, 0, 0, 0 ); }
- | expr '%' expr
- { PROMOTE($1,$3); $$ = New_BinOp( TYPE($1), $1, '%', $3 );
- TEST($$); }
- | expr '+' expr
- { PROMOTE($1,$3); $$ = New_BinOp( TYPE($1), $1, '+', $3 );
- TEST($$); }
- | expr '-' expr
- { PROMOTE($1,$3); $$ = New_BinOp( TYPE($1), $1, '-', $3 );
- TEST($$); }
- | expr '*' expr
- { PROMOTE($1,$3); $$ = New_BinOp( TYPE($1), $1, '*', $3 );
- TEST($$); }
- | expr '/' expr
- { PROMOTE($1,$3); $$ = New_BinOp( TYPE($1), $1, '/', $3 );
- TEST($$); }
- | expr POWER expr
- { PROMOTE($1,$3); $$ = New_BinOp( TYPE($1), $1, POWER, $3 );
- TEST($$); }
- | '+' expr %prec UMINUS
- { $$ = $2; }
- | '-' expr %prec UMINUS
- { $$ = New_Unary( TYPE($2), UMINUS, $2 ); TEST($$); }
- | '(' expr ')'
- { $$ = $2; }
- | expr '*' bexpr
- { $3 = New_Unary( TYPE($1), 0, $3 );
- $$ = New_BinOp( TYPE($1), $1, '*', $3 );
- TEST($$); }
- | bexpr '*' expr
- { $1 = New_Unary( TYPE($3), 0, $1 );
- $$ = New_BinOp( TYPE($3), $1, '*', $3 );
- TEST($$); }
- | bexpr '?' expr ':' expr
- {
- PROMOTE($3,$5);
- if( ! Test_Dims($3,$5) ) {
- yyerror("Incompatible dimensions in '?:' arguments");
- YYERROR;
- }
- $$ = New_Func( 0, ifthenelse_fct, 3, $3, $5, $1,
- 0, 0, 0, 0 );
- TEST($$);
- if( SIZE($3)<SIZE($5) ) Copy_Dims($$, $5);
- TYPE($1) = TYPE($3);
- if( ! Test_Dims($1,$$) ) {
- yyerror("Incompatible dimensions in '?:' condition");
- YYERROR;
- }
- TYPE($1) = BOOLEAN;
- if( SIZE($$)<SIZE($1) ) Copy_Dims($$, $1);
- }
- | bexpr '?' bexpr ':' expr
- {
- PROMOTE($3,$5);
- if( ! Test_Dims($3,$5) ) {
- yyerror("Incompatible dimensions in '?:' arguments");
- YYERROR;
- }
- $$ = New_Func( 0, ifthenelse_fct, 3, $3, $5, $1,
- 0, 0, 0, 0 );
- TEST($$);
- if( SIZE($3)<SIZE($5) ) Copy_Dims($$, $5);
- TYPE($1) = TYPE($3);
- if( ! Test_Dims($1,$$) ) {
- yyerror("Incompatible dimensions in '?:' condition");
- YYERROR;
- }
- TYPE($1) = BOOLEAN;
- if( SIZE($$)<SIZE($1) ) Copy_Dims($$, $1);
- }
- | bexpr '?' expr ':' bexpr
- {
- PROMOTE($3,$5);
- if( ! Test_Dims($3,$5) ) {
- yyerror("Incompatible dimensions in '?:' arguments");
- YYERROR;
- }
- $$ = New_Func( 0, ifthenelse_fct, 3, $3, $5, $1,
- 0, 0, 0, 0 );
- TEST($$);
- if( SIZE($3)<SIZE($5) ) Copy_Dims($$, $5);
- TYPE($1) = TYPE($3);
- if( ! Test_Dims($1,$$) ) {
- yyerror("Incompatible dimensions in '?:' condition");
- YYERROR;
- }
- TYPE($1) = BOOLEAN;
- if( SIZE($$)<SIZE($1) ) Copy_Dims($$, $1);
- }
- | FUNCTION ')'
- { if (FSTRCMP($1,"RANDOM(") == 0) { /* Scalar RANDOM() */
- srand( (unsigned int) time(NULL) );
- $$ = New_Func( DOUBLE, rnd_fct, 0, 0, 0, 0, 0, 0, 0, 0 );
- } else if (FSTRCMP($1,"RANDOMN(") == 0) {/*Scalar RANDOMN()*/
- srand( (unsigned int) time(NULL) );
- $$ = New_Func( DOUBLE, gasrnd_fct, 0, 0, 0, 0, 0, 0, 0, 0 );
- } else {
- yyerror("Function() not supported");
- YYERROR;
- }
- TEST($$);
- }
- | FUNCTION bexpr ')'
- { if (FSTRCMP($1,"SUM(") == 0) {
- $$ = New_Func( LONG, sum_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- } else if (FSTRCMP($1,"NELEM(") == 0) {
- $$ = New_Const( LONG, &( SIZE($2) ), sizeof(long) );
- } else if (FSTRCMP($1,"ACCUM(") == 0) {
- long zero = 0;
- $$ = New_BinOp( LONG , $2, ACCUM, New_Const( LONG, &zero, sizeof(zero) ));
- } else {
- yyerror("Function(bool) not supported");
- YYERROR;
- }
- TEST($$);
- }
- | FUNCTION sexpr ')'
- { if (FSTRCMP($1,"NELEM(") == 0) {
- $$ = New_Const( LONG, &( SIZE($2) ), sizeof(long) );
- } else if (FSTRCMP($1,"NVALID(") == 0) {
- $$ = New_Func( LONG, nonnull_fct, 1, $2,
- 0, 0, 0, 0, 0, 0 );
- } else {
- yyerror("Function(str) not supported");
- YYERROR;
- }
- TEST($$);
- }
- | FUNCTION bits ')'
- { if (FSTRCMP($1,"NELEM(") == 0) {
- $$ = New_Const( LONG, &( SIZE($2) ), sizeof(long) );
- } else if (FSTRCMP($1,"NVALID(") == 0) { /* Bit arrays do not have NULL */
- $$ = New_Const( LONG, &( SIZE($2) ), sizeof(long) );
- } else if (FSTRCMP($1,"SUM(") == 0) {
- $$ = New_Func( LONG, sum_fct, 1, $2,
- 0, 0, 0, 0, 0, 0 );
- } else if (FSTRCMP($1,"MIN(") == 0) {
- $$ = New_Func( TYPE($2), /* Force 1D result */
- min1_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- /* Note: $2 is a vector so the result can never
- be a constant. Therefore it will never be set
- inside New_Func(), and it is safe to set SIZE() */
- SIZE($$) = 1;
- } else if (FSTRCMP($1,"ACCUM(") == 0) {
- long zero = 0;
- $$ = New_BinOp( LONG , $2, ACCUM, New_Const( LONG, &zero, sizeof(zero) ));
- } else if (FSTRCMP($1,"MAX(") == 0) {
- $$ = New_Func( TYPE($2), /* Force 1D result */
- max1_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- /* Note: $2 is a vector so the result can never
- be a constant. Therefore it will never be set
- inside New_Func(), and it is safe to set SIZE() */
- SIZE($$) = 1;
- } else {
- yyerror("Function(bits) not supported");
- YYERROR;
- }
- TEST($$);
- }
- | FUNCTION expr ')'
- { if (FSTRCMP($1,"SUM(") == 0)
- $$ = New_Func( TYPE($2), sum_fct, 1, $2,
- 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"AVERAGE(") == 0)
- $$ = New_Func( DOUBLE, average_fct, 1, $2,
- 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"STDDEV(") == 0)
- $$ = New_Func( DOUBLE, stddev_fct, 1, $2,
- 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"MEDIAN(") == 0)
- $$ = New_Func( TYPE($2), median_fct, 1, $2,
- 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"NELEM(") == 0)
- $$ = New_Const( LONG, &( SIZE($2) ), sizeof(long) );
- else if (FSTRCMP($1,"NVALID(") == 0)
- $$ = New_Func( LONG, nonnull_fct, 1, $2,
- 0, 0, 0, 0, 0, 0 );
- else if ((FSTRCMP($1,"ACCUM(") == 0) && (TYPE($2) == LONG)) {
- long zero = 0;
- $$ = New_BinOp( LONG , $2, ACCUM, New_Const( LONG, &zero, sizeof(zero) ));
- } else if ((FSTRCMP($1,"ACCUM(") == 0) && (TYPE($2) == DOUBLE)) {
- double zero = 0;
- $$ = New_BinOp( DOUBLE , $2, ACCUM, New_Const( DOUBLE, &zero, sizeof(zero) ));
- } else if ((FSTRCMP($1,"SEQDIFF(") == 0) && (TYPE($2) == LONG)) {
- long zero = 0;
- $$ = New_BinOp( LONG , $2, DIFF, New_Const( LONG, &zero, sizeof(zero) ));
- } else if ((FSTRCMP($1,"SEQDIFF(") == 0) && (TYPE($2) == DOUBLE)) {
- double zero = 0;
- $$ = New_BinOp( DOUBLE , $2, DIFF, New_Const( DOUBLE, &zero, sizeof(zero) ));
- } else if (FSTRCMP($1,"ABS(") == 0)
- $$ = New_Func( 0, abs_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"MIN(") == 0)
- $$ = New_Func( TYPE($2), /* Force 1D result */
- min1_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"MAX(") == 0)
- $$ = New_Func( TYPE($2), /* Force 1D result */
- max1_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"RANDOM(") == 0) { /* Vector RANDOM() */
- srand( (unsigned int) time(NULL) );
- $$ = New_Func( 0, rnd_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- TEST($$);
- TYPE($$) = DOUBLE;
- } else if (FSTRCMP($1,"RANDOMN(") == 0) {
- srand( (unsigned int) time(NULL) ); /* Vector RANDOMN() */
- $$ = New_Func( 0, gasrnd_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- TEST($$);
- TYPE($$) = DOUBLE;
- }
- else { /* These all take DOUBLE arguments */
- if( TYPE($2) != DOUBLE ) $2 = New_Unary( DOUBLE, 0, $2 );
- if (FSTRCMP($1,"SIN(") == 0)
- $$ = New_Func( 0, sin_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"COS(") == 0)
- $$ = New_Func( 0, cos_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"TAN(") == 0)
- $$ = New_Func( 0, tan_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"ARCSIN(") == 0
- || FSTRCMP($1,"ASIN(") == 0)
- $$ = New_Func( 0, asin_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"ARCCOS(") == 0
- || FSTRCMP($1,"ACOS(") == 0)
- $$ = New_Func( 0, acos_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"ARCTAN(") == 0
- || FSTRCMP($1,"ATAN(") == 0)
- $$ = New_Func( 0, atan_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"SINH(") == 0)
- $$ = New_Func( 0, sinh_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"COSH(") == 0)
- $$ = New_Func( 0, cosh_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"TANH(") == 0)
- $$ = New_Func( 0, tanh_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"EXP(") == 0)
- $$ = New_Func( 0, exp_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"LOG(") == 0)
- $$ = New_Func( 0, log_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"LOG10(") == 0)
- $$ = New_Func( 0, log10_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"SQRT(") == 0)
- $$ = New_Func( 0, sqrt_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"ROUND(") == 0)
- $$ = New_Func( 0, round_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"FLOOR(") == 0)
- $$ = New_Func( 0, floor_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"CEIL(") == 0)
- $$ = New_Func( 0, ceil_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
- else if (FSTRCMP($1,"RANDOMP(") == 0) {
- srand( (unsigned int) time(NULL) );
- $$ = New_Func( 0, poirnd_fct, 1, $2,
- 0, 0, 0, 0, 0, 0 );
- TYPE($$) = LONG;
- } else {
- yyerror("Function(expr) not supported");
- YYERROR;
- }
- }
- TEST($$);
- }
- | IFUNCTION sexpr ',' sexpr ')'
- {
- if (FSTRCMP($1,"STRSTR(") == 0) {
- $$ = New_Func( LONG, strpos_fct, 2, $2, $4, 0,
- 0, 0, 0, 0 );
- TEST($$);
- }
- }
- | FUNCTION expr ',' expr ')'
- {
- if (FSTRCMP($1,"DEFNULL(") == 0) {
- if( SIZE($2)>=SIZE($4) && Test_Dims( $2, $4 ) ) {
- PROMOTE($2,$4);
- $$ = New_Func( 0, defnull_fct, 2, $2, $4, 0,
- 0, 0, 0, 0 );
- TEST($$);
- } else {
- yyerror("Dimensions of DEFNULL arguments "
- "are not compatible");
- YYERROR;
- }
- } else if (FSTRCMP($1,"ARCTAN2(") == 0) {
- if( TYPE($2) != DOUBLE ) $2 = New_Unary( DOUBLE, 0, $2 );
- if( TYPE($4) != DOUBLE ) $4 = New_Unary( DOUBLE, 0, $4 );
- if( Test_Dims( $2, $4 ) ) {
- $$ = New_Func( 0, atan2_fct, 2, $2, $4, 0, 0, 0, 0, 0 );
- TEST($$);
- if( SIZE($2)<SIZE($4) ) Copy_Dims($$, $4);
- } else {
- yyerror("Dimensions of arctan2 arguments "
- "are not compatible");
- YYERROR;
- }
- } else if (FSTRCMP($1,"MIN(") == 0) {
- PROMOTE( $2, $4 );
- if( Test_Dims( $2, $4 ) ) {
- $$ = New_Func( 0, min2_fct, 2, $2, $4, 0, 0, 0, 0, 0 );
- TEST($$);
- if( SIZE($2)<SIZE($4) ) Copy_Dims($$, $4);
- } else {
- yyerror("Dimensions of min(a,b) arguments "
- "are not compatible");
- YYERROR;
- }
- } else if (FSTRCMP($1,"MAX(") == 0) {
- PROMOTE( $2, $4 );
- if( Test_Dims( $2, $4 ) ) {
- $$ = New_Func( 0, max2_fct, 2, $2, $4, 0, 0, 0, 0, 0 );
- TEST($$);
- if( SIZE($2)<SIZE($4) ) Copy_Dims($$, $4);
- } else {
- yyerror("Dimensions of max(a,b) arguments "
- "are not compatible");
- YYERROR;
- }
- #if 0
- } else if (FSTRCMP($1,"STRSTR(") == 0) {
- if( TYPE($2) != STRING || TYPE($4) != STRING) {
- yyerror("Arguments to strstr(s,r) must be strings");
- YYERROR;
- }
- $$ = New_Func( LONG, strpos_fct, 2, $2, $4, 0,
- 0, 0, 0, 0 );
- TEST($$);
- #endif
- } else {
- yyerror("Function(expr,expr) not supported");
- YYERROR;
- }
- }
- | FUNCTION expr ',' expr ',' expr ',' expr ')'
- {
- if (FSTRCMP($1,"ANGSEP(") == 0) {
- if( TYPE($2) != DOUBLE ) $2 = New_Unary( DOUBLE, 0, $2 );
- if( TYPE($4) != DOUBLE ) $4 = New_Unary( DOUBLE, 0, $4 );
- if( TYPE($6) != DOUBLE ) $6 = New_Unary( DOUBLE, 0, $6 );
- if( TYPE($8) != DOUBLE ) $8 = New_Unary( DOUBLE, 0, $8 );
- if( Test_Dims( $2, $4 ) && Test_Dims( $4, $6 ) &&
- Test_Dims( $6, $8 ) ) {
- $$ = New_Func( 0, angsep_fct, 4, $2, $4, $6, $8,0,0,0 );
- TEST($$);
- if( SIZE($2)<SIZE($4) ) Copy_Dims($$, $4);
- if( SIZE($4)<SIZE($6) ) Copy_Dims($$, $6);
- if( SIZE($6)<SIZE($8) ) Copy_Dims($$, $8);
- } else {
- yyerror("Dimensions of ANGSEP arguments "
- "are not compatible");
- YYERROR;
- }
- } else {
- yyerror("Function(expr,expr,expr,expr) not supported");
- YYERROR;
- }
- }
- | expr '[' expr ']'
- { $$ = New_Deref( $1, 1, $3, 0, 0, 0, 0 ); TEST($$); }
- | expr '[' expr ',' expr ']'
- { $$ = New_Deref( $1, 2, $3, $5, 0, 0, 0 ); TEST($$); }
- | expr '[' expr ',' expr ',' expr ']'
- { $$ = New_Deref( $1, 3, $3, $5, $7, 0, 0 ); TEST($$); }
- | expr '[' expr ',' expr ',' expr ',' expr ']'
- { $$ = New_Deref( $1, 4, $3, $5, $7, $9, 0 ); TEST($$); }
- | expr '[' expr ',' expr ',' expr ',' expr ',' expr ']'
- { $$ = New_Deref( $1, 5, $3, $5, $7, $9, $11 ); TEST($$); }
- | INTCAST expr
- { $$ = New_Unary( LONG, INTCAST, $2 ); TEST($$); }
- | INTCAST bexpr
- { $$ = New_Unary( LONG, INTCAST, $2 ); TEST($$); }
- | FLTCAST expr
- { $$ = New_Unary( DOUBLE, FLTCAST, $2 ); TEST($$); }
- | FLTCAST bexpr
- { $$ = New_Unary( DOUBLE, FLTCAST, $2 ); TEST($$); }
- ;
- bexpr: BOOLEAN
- { $$ = New_Const( BOOLEAN, &($1), sizeof(char) ); TEST($$); }
- | BCOLUMN
- { $$ = New_Column( $1 ); TEST($$); }
- | BCOLUMN '{' expr '}'
- {
- if( TYPE($3) != LONG
- || OPER($3) != CONST_OP ) {
- yyerror("Offset argument must be a constant integer");
- YYERROR;
- }
- $$ = New_Offset( $1, $3 ); TEST($$);
- }
- | bits EQ bits
- { $$ = New_BinOp( BOOLEAN, $1, EQ, $3 ); TEST($$);
- SIZE($$) = 1; }
- | bits NE bits
- { $$ = New_BinOp( BOOLEAN, $1, NE, $3 ); TEST($$);
- SIZE($$) = 1; }
- | bits LT bits
- { $$ = New_BinOp( BOOLEAN, $1, LT, $3 ); TEST($$);
- SIZE($$) = 1; }
- | bits LTE bits
- { $$ = New_BinOp( BOOLEAN, $1, LTE, $3 ); TEST($$);
- SIZE($$) = 1; }
- | bits GT bits
- { $$ = New_BinOp( BOOLEAN, $1, GT, $3 ); TEST($$);
- SIZE($$) = 1; }
- | bits GTE bits
- { $$ = New_BinOp( BOOLEAN, $1, GTE, $3 ); TEST($$);
- SIZE($$) = 1; }
- | expr GT expr
- { PROMOTE($1,$3); $$ = New_BinOp( BOOLEAN, $1, GT, $3 );
- TEST($$); }
- | expr LT expr
- { PROMOTE($1,$3); $$ = New_BinOp( BOOLEAN, $1, LT, $3 );
- TEST($$); }
- | expr GTE expr
- { PROMOTE($1,$3); $$ = New_BinOp( BOOLEAN, $1, GTE, $3 );
- TEST($$); }
- | expr LTE expr
- { PROMOTE($1,$3); $$ = New_BinOp( BOOLEAN, $1, LTE, $3 );
- TEST($$); }
- | expr '~' expr
- { PROMOTE($1,$3); $$ = New_BinOp( BOOLEAN, $1, '~', $3 );
- TEST($$); }
- | expr EQ expr
- { PROMOTE($1,$3); $$ = New_BinOp( BOOLEAN, $1, EQ, $3 );
- TEST($$); }
- | expr NE expr
- { PROMOTE($1,$3); $$ = New_BinOp( BOOLEAN, $1, NE, $3 );
- TEST($$); }
- | sexpr EQ sexpr
- { $$ = New_BinOp( BOOLEAN, $1, EQ, $3 ); TEST($$);
- SIZE($$) = 1; }
- | sexpr NE sexpr
- { $$ = New_BinOp( BOOLEAN, $1, NE, $3 ); TEST($$);
- SIZE($$) = 1; }
- | sexpr GT sexpr
- { $$ = New_BinOp( BOOLEAN, $1, GT, $3 ); TEST($$);
- SIZE($$) = 1; }
- | sexpr GTE sexpr
- { $$ = New_BinOp( BOOLEAN, $1, GTE, $3 ); TEST($$);
- SIZE($$) = 1; }
- | sexpr LT sexpr
- { $$ = New_BinOp( BOOLEAN, $1, LT, $3 ); TEST($$);
- SIZE($$) = 1; }
- | sexpr LTE sexpr
- { $$ = New_BinOp( BOOLEAN, $1, LTE, $3 ); TEST($$);
- SIZE($$) = 1; }
- | bexpr AND bexpr
- { $$ = New_BinOp( BOOLEAN, $1, AND, $3 ); TEST($$); }
- | bexpr OR bexpr
- { $$ = New_BinOp( BOOLEAN, $1, OR, $3 ); TEST($$); }
- | bexpr EQ bexpr
- { $$ = New_BinOp( BOOLEAN, $1, EQ, $3 ); TEST($$); }
- | bexpr NE bexpr
- { $$ = New_BinOp( BOOLEAN, $1, NE, $3 ); TEST($$); }
- | expr '=' expr ':' expr
- { PROMOTE($1,$3); PROMOTE($1,$5); PROMOTE($3,$5);
- $3 = New_BinOp( BOOLEAN, $3, LTE, $1 );
- $5 = New_BinOp( BOOLEAN, $1, LTE, $5 );
- $$ = New_BinOp( BOOLEAN, $3, AND, $5 );
- TEST($$); }
- | bexpr '?' bexpr ':' bexpr
- {
- if( ! Test_Dims($3,$5) ) {
- yyerror("Incompatible dimensions in '?:' arguments");
- YYERROR;
- }
- $$ = New_Func( 0, ifthenelse_fct, 3, $3, $5, $1,
- 0, 0, 0, 0 );
- TEST($$);
- if( SIZE($3)<SIZE($5) ) Copy_Dims($$, $5);
- if( ! Test_Dims($1,$$) ) {
- yyerror("Incompatible dimensions in '?:' condition");
- YYERROR;
- }
- if( SIZE($$)<SIZE($1) ) Copy_Dims($$, $1);
- }
- | BFUNCTION expr ')'
- {
- if (FSTRCMP($1,"ISNULL(") == 0) {
- $$ = New_Func( 0, isnull_fct, 1, $2, 0, 0,
- 0, 0, 0, 0 );
- TEST($$);
- /* Use expression's size, but return BOOLEAN */
- TYPE($$) = BOOLEAN;
- } else {
- yyerror("Boolean Function(expr) not supported");
- YYERROR;
- }
- }
- | BFUNCTION bexpr ')'
- {
- if (FSTRCMP($1,"ISNULL(") == 0) {
- $$ = New_Func( 0, isnull_fct, 1, $2, 0, 0,
- 0, 0, 0, 0 );
- TEST($$);
- /* Use expression's size, but return BOOLEAN */
- TYPE($$) = BOOLEAN;
- } else {
- yyerror("Boolean Function(expr) not supported");
- YYERROR;
- }
- }
- | BFUNCTION sexpr ')'
- {
- if (FSTRCMP($1,"ISNULL(") == 0) {
- $$ = New_Func( BOOLEAN, isnull_fct, 1, $2, 0, 0,
- 0, 0, 0, 0 );
- TEST($$);
- } else {
- yyerror("Boolean Function(expr) not supported");
- YYERROR;
- }
- }
- | FUNCTION bexpr ',' bexpr ')'
- {
- if (FSTRCMP($1,"DEFNULL(") == 0) {
- if( SIZE($2)>=SIZE($4) && Test_Dims( $2, $4 ) ) {
- $$ = New_Func( 0, defnull_fct, 2, $2, $4, 0,
- 0, 0, 0, 0 );
- TEST($$);
- } else {
- yyerror("Dimensions of DEFNULL arguments are not compatible");
- YYERROR;
- }
- } else {
- yyerror("Boolean Function(expr,expr) not supported");
- YYERROR;
- }
- }
- | BFUNCTION expr ',' expr ',' expr ')'
- {
- if( TYPE($2) != DOUBLE ) $2 = New_Unary( DOUBLE, 0, $2 );
- if( TYPE($4) != DOUBLE ) $4 = New_Unary( DOUBLE, 0, $4 );
- if( TYPE($6) != DOUBLE ) $6 = New_Unary( DOUBLE, 0, $6 );
- if( ! (Test_Dims( $2, $4 ) && Test_Dims( $4, $6 ) ) ) {
- yyerror("Dimensions of NEAR arguments "
- "are not compatible");
- YYERROR;
- } else {
- if (FSTRCMP($1,"NEAR(") == 0) {
- $$ = New_Func( BOOLEAN, near_fct, 3, $2, $4, $6,
- 0, 0, 0, 0 );
- } else {
- yyerror("Boolean Function not supported");
- YYERROR;
- }
- TEST($$);
- if( SIZE($$)<SIZE($2) ) Copy_Dims($$, $2);
- if( SIZE($2)<SIZE($4) ) Copy_Dims($$, $4);
- if( SIZE($4)<SIZE($6) ) Copy_Dims($$, $6);
- }
- }
- | BFUNCTION expr ',' expr ',' expr ',' expr ',' expr ')'
- {
- if( TYPE($2) != DOUBLE ) $2 = New_Unary( DOUBLE, 0, $2 );
- if( TYPE($4) != DOUBLE ) $4 = New_Unary( DOUBLE, 0, $4 );
- if( TYPE($6) != DOUBLE ) $6 = New_Unary( DOUBLE, 0, $6 );
- if( TYPE($8) != DOUBLE ) $8 = New_Unary( DOUBLE, 0, $8 );
- if( TYPE($10)!= DOUBLE ) $10= New_Unary( DOUBLE, 0, $10);
- if( ! (Test_Dims( $2, $4 ) && Test_Dims( $4, $6 ) &&
- Test_Dims( $6, $8 ) && Test_Dims( $8, $10 )) ) {
- yyerror("Dimensions of CIRCLE arguments "
- "are not compatible");
- YYERROR;
- } else {
- if (FSTRCMP($1,"CIRCLE(") == 0) {
- $$ = New_Func( BOOLEAN, circle_fct, 5, $2, $4, $6, $8,
- $10, 0, 0 );
- } else {
- yyerror("Boolean Function not supported");
- YYERROR;
- }
- TEST($$);
- if( SIZE($$)<SIZE($2) ) Copy_Dims($$, $2);
- if( SIZE($2)<SIZE($4) ) Copy_Dims($$, $4);
- if( SIZE($4)<SIZE($6) ) Copy_Dims($$, $6);
- if( SIZE($6)<SIZE($8) ) Copy_Dims($$, $8);
- if( SIZE($8)<SIZE($10) ) Copy_Dims($$, $10);
- }
- }
- | BFUNCTION expr ',' expr ',' expr ',' expr ',' expr ',' expr ',' expr ')'
- {
- if( TYPE($2) != DOUBLE ) $2 = New_Unary( DOUBLE, 0, $2 );
- if( TYPE($4) != DOUBLE ) $4 = New_Unary( DOUBLE, 0, $4 );
- if( TYPE($6) != DOUBLE ) $6 = New_Unary( DOUBLE, 0, $6 );
- if( TYPE($8) != DOUBLE ) $8 = New_Unary( DOUBLE, 0, $8 );
- if( TYPE($10)!= DOUBLE ) $10= New_Unary( DOUBLE, 0, $10);
- if( TYPE($12)!= DOUBLE ) $12= New_Unary( DOUBLE, 0, $12);
- if( TYPE($14)!= DOUBLE ) $14= New_Unary( DOUBLE, 0, $14);
- if( ! (Test_Dims( $2, $4 ) && Test_Dims( $4, $6 ) &&
- Test_Dims( $6, $8 ) && Test_Dims( $8, $10 ) &&
- Test_Dims($10,$12 ) && Test_Dims($12, $14 ) ) ) {
- yyerror("Dimensions of BOX or ELLIPSE arguments "
- "are not compatible");
- YYERROR;
- } else {
- if (FSTRCMP($1,"BOX(") == 0) {
- $$ = New_Func( BOOLEAN, box_fct, 7, $2, $4, $6, $8,
- $10, $12, $14 );
- } else if (FSTRCMP($1,"ELLIPSE(") == 0) {
- $$ = New_Func( BOOLEAN, elps_fct, 7, $2, $4, $6, $8,
- $10, $12, $14 );
- } else {
- yyerror("SAO Image Function not supported");
- YYERROR;
- }
- TEST($$);
- if( SIZE($$)<SIZE($2) ) Copy_Dims($$, $2);
- if( SIZE($2)<SIZE($4) ) Copy_Dims($$, $4);
- if( SIZE($4)<SIZE($6) ) Copy_Dims($$, $6);
- if( SIZE($6)<SIZE($8) ) Copy_Dims($$, $8);
- if( SIZE($8)<SIZE($10) ) Copy_Dims($$, $10);
- if( SIZE($10)<SIZE($12) ) Copy_Dims($$, $12);
- if( SIZE($12)<SIZE($14) ) Copy_Dims($$, $14);
- }
- }
- | GTIFILTER ')'
- { /* Use defaults for all elements */
- $$ = New_GTI( "", -99, "*START*", "*STOP*" );
- TEST($$); }
- | GTIFILTER STRING ')'
- { /* Use defaults for all except filename */
- $$ = New_GTI( $2, -99, "*START*", "*STOP*" );
- TEST($$); }
- | GTIFILTER STRING ',' expr ')'
- { $$ = New_GTI( $2, $4, "*START*", "*STOP*" );
- TEST($$); }
- | GTIFILTER STRING ',' expr ',' STRING ',' STRING ')'
- { $$ = New_GTI( $2, $4, $6, $8 );
- TEST($$); }
- | REGFILTER STRING ')'
- { /* Use defaults for all except filename */
- $$ = New_REG( $2, -99, -99, "" );
- TEST($$); }
- | REGFILTER STRING ',' expr ',' expr ')'
- { $$ = New_REG( $2, $4, $6, "" );
- TEST($$); }
- | REGFILTER STRING ',' expr ',' expr ',' STRING ')'
- { $$ = New_REG( $2, $4, $6, $8 );
- TEST($$); }
- | bexpr '[' expr ']'
- { $$ = New_Deref( $1, 1, $3, 0, 0, 0, 0 ); TEST($$); }
- | bexpr '[' expr ',' expr ']'
- { $$ = New_Deref( $1, 2, $3, $5, 0, 0, 0 ); TEST($$); }
- | bexpr '[' expr ',' expr ',' expr ']'
- { $$ = New_Deref( $1, 3, $3, $5, $7, 0, 0 ); TEST($$); }
- | bexpr '[' expr ',' expr ',' expr ',' expr ']'
- { $$ = New_Deref( $1, 4, $3, $5, $7, $9, 0 ); TEST($$); }
- | bexpr '[' expr ',' expr ',' expr ',' expr ',' expr ']'
- { $$ = New_Deref( $1, 5, $3, $5, $7, $9, $11 ); TEST($$); }
- | NOT bexpr
- { $$ = New_Unary( BOOLEAN, NOT, $2 ); TEST($$); }
- | '(' bexpr ')'
- { $$ = $2; }
- ;
- sexpr: STRING
- { $$ = New_Const( STRING, $1, strlen($1)+1 ); TEST($$);
- SIZE($$) = strlen($1); }
- | SCOLUMN
- { $$ = New_Column( $1 ); TEST($$); }
- | SCOLUMN '{' expr '}'
- {
- if( TYPE($3) != LONG
- || OPER($3) != CONST_OP ) {
- yyerror("Offset argument must be a constant integer");
- YYERROR;
- }
- $$ = New_Offset( $1, $3 ); TEST($$);
- }
- | SNULLREF
- { $$ = New_Func( STRING, null_fct, 0, 0, 0, 0, 0, 0, 0, 0 ); }
- | '(' sexpr ')'
- { $$ = $2; }
- | sexpr '+' sexpr
- {
- if (SIZE($1)+SIZE($3) >= MAX_STRLEN) {
- yyerror("Combined string size exceeds " MAX_STRLEN_S " characters");
- YYERROR;
- }
- $$ = New_BinOp( STRING, $1, '+', $3 ); TEST($$);
- SIZE($$) = SIZE($1) + SIZE($3);
- }
- | bexpr '?' sexpr ':' sexpr
- {
- int outSize;
- if( SIZE($1)!=1 ) {
- yyerror("Cannot have a vector string column");
- YYERROR;
- }
- /* Since the output can be calculated now, as a constant
- scalar, we must precalculate the output size, in
- order to avoid an overflow. */
- outSize = SIZE($3);
- if (SIZE($5) > outSize) outSize = SIZE($5);
- $$ = New_FuncSize( 0, ifthenelse_fct, 3, $3, $5, $1,
- 0, 0, 0, 0, outSize);
-
- TEST($$);
- if( SIZE($3)<SIZE($5) ) Copy_Dims($$, $5);
- }
- | FUNCTION sexpr ',' sexpr ')'
- {
- if (FSTRCMP($1,"DEFNULL(") == 0) {
- int outSize;
- /* Since the output can be calculated now, as a constant
- scalar, we must precalculate the output size, in
- order to avoid an overflow. */
- outSize = SIZE($2);
- if (SIZE($4) > outSize) outSize = SIZE($4);
-
- $$ = New_FuncSize( 0, defnull_fct, 2, $2, $4, 0,
- 0, 0, 0, 0, outSize );
- TEST($$);
- if( SIZE($4)>SIZE($2) ) SIZE($$) = SIZE($4);
- } else {
- yyerror("Function(string,string) not supported");
- YYERROR;
- }
- }
- | FUNCTION sexpr ',' expr ',' expr ')'
- {
- if (FSTRCMP($1,"STRMID(") == 0) {
- int len;
- if( TYPE($4) != LONG || SIZE($4) != 1 ||
- TYPE($6) != LONG || SIZE($6) != 1) {
- yyerror("When using STRMID(S,P,N), P and N must be integers (and not vector columns)");
- YYERROR;
- }
- if (OPER($6) == CONST_OP) {
- /* Constant value: use that directly */
- len = (gParse.Nodes[$6].value.data.lng);
- } else {
- /* Variable value: use the maximum possible (from $2) */
- len = SIZE($2);
- }
- if (len <= 0 || len >= MAX_STRLEN) {
- yyerror("STRMID(S,P,N), N must be 1-" MAX_STRLEN_S);
- YYERROR;
- }
- $$ = New_FuncSize( 0, strmid_fct, 3, $2, $4,$6,0,0,0,0,len);
- TEST($$);
- } else {
- yyerror("Function(string,expr,expr) not supported");
- YYERROR;
- }
- }
- ;
- %%
- /*************************************************************************/
- /* Start of "New" routines which build the expression Nodal structure */
- /*************************************************************************/
- static int Alloc_Node( void )
- {
- /* Use this for allocation to guarantee *Nodes */
- Node *newNodePtr; /* survives on failure, making it still valid */
- /* while working our way out of this error */
- if( gParse.nNodes == gParse.nNodesAlloc ) {
- if( gParse.Nodes ) {
- gParse.nNodesAlloc += gParse.nNodesAlloc;
- newNodePtr = (Node *)realloc( gParse.Nodes,
- sizeof(Node)*gParse.nNodesAlloc );
- } else {
- gParse.nNodesAlloc = 100;
- newNodePtr = (Node *)malloc ( sizeof(Node)*gParse.nNodesAlloc );
- }
- if( newNodePtr ) {
- gParse.Nodes = newNodePtr;
- } else {
- gParse.status = MEMORY_ALLOCATION;
- return( -1 );
- }
- }
- return ( gParse.nNodes++ );
- }
- static void Free_Last_Node( void )
- {
- if( gParse.nNodes ) gParse.nNodes--;
- }
- static int New_Const( int returnType, void *value, long len )
- {
- Node *this;
- int n;
- n = Alloc_Node();
- if( n>=0 ) {
- this = gParse.Nodes + n;
- this->operation = CONST_OP; /* Flag a constant */
- this->DoOp = NULL;
- this->nSubNodes = 0;
- this->type = returnType;
- memcpy( &(this->value.data), value, len );
- this->value.undef = NULL;
- this->value.nelem = 1;
- this->value.naxis = 1;
- this->value.naxes[0] = 1;
- }
- return(n);
- }
- static int New_Column( int ColNum )
- {
- Node *this;
- int n, i;
- n = Alloc_Node();
- if( n>=0 ) {
- this = gParse.Nodes + n;
- this->operation = -ColNum;
- this->DoOp = NULL;
- this->nSubNodes = 0;
- this->type = gParse.varData[ColNum].type;
- this->value.nelem = gParse.varData[ColNum].nelem;
- this->value.naxis = gParse.varData[ColNum].naxis;
- for( i=0; i<gParse.varData[ColNum].naxis; i++ )
- this->value.naxes[i] = gParse.varData[ColNum].naxes[i];
- }
- return(n);
- }
- static int New_Offset( int ColNum, int offsetNode )
- {
- Node *this;
- int n, i, colNode;
- colNode = New_Column( ColNum );
- if( colNode<0 ) return(-1);
- n = Alloc_Node();
- if( n>=0 ) {
- this = gParse.Nodes + n;
- this->operation = '{';
- this->DoOp = Do_Offset;
- this->nSubNodes = 2;
- this->SubNodes[0] = colNode;
- this->SubNodes[1] = offsetNode;
- this->type = gParse.varData[ColNum].type;
- this->value.nelem = gParse.varData[ColNum].nelem;
- this->value.naxis = gParse.varData[ColNum].naxis;
- for( i=0; i<gParse.varData[ColNum].naxis; i++ )
- this->value.naxes[i] = gParse.varData[ColNum].naxes[i];
- }
- return(n);
- }
- static int New_Unary( int returnType, int Op, int Node1 )
- {
- Node *this, *that;
- int i,n;
- if( Node1<0 ) return(-1);
- that = gParse.Nodes + Node1;
- if( !Op ) Op = returnType;
- if( (Op==DOUBLE || Op==FLTCAST) && that->type==DOUBLE ) return( Node1 );
- if( (Op==LONG || Op==INTCAST) && that->type==LONG ) return( Node1 );
- if( (Op==BOOLEAN ) && that->type==BOOLEAN ) return( Node1 );
-
- n = Alloc_Node();
- if( n>=0 ) {
- this = gParse.Nodes + n;
- this->operation = Op;
- this->DoOp = Do_Unary;
- this->nSubNodes = 1;
- this->SubNodes[0] = Node1;
- this->type = returnType;
- that = gParse.Nodes + Node1; /* Reset in case .Nodes mv'd */
- this->value.nelem = that->value.nelem;
- this->value.naxis = that->value.naxis;
- for( i=0; i<that->value.naxis; i++ )
- this->value.naxes[i] = that->value.naxes[i];
- if( that->operation==CONST_OP ) this->DoOp( this );
- }
- return( n );
- }
- static int New_BinOp( int returnType, int Node1, int Op, int Node2 )
- {
- Node *this,*that1,*that2;
- int n,i,constant;
- if( Node1<0 || Node2<0 ) return(-1);
- n = Alloc_Node();
- if( n>=0 ) {
- this = gParse.Nodes + n;
- this->operation = Op;
- this->nSubNodes = 2;
- this->SubNodes[0]= Node1;
- this->SubNodes[1]= Node2;
- this->type = returnType;
- that1 = gParse.Nodes + Node1;
- that2 = gParse.Nodes + Node2;
- constant = (that1->operation==CONST_OP
- && that2->operation==CONST_OP);
- if( that1->type!=STRING && that1->type!=BITSTR )
- if( !Test_Dims( Node1, Node2 ) ) {
- Free_Last_Node();
- yyerror("Array sizes/dims do not match for binary operator");
- return(-1);
- }
- if( that1->value.nelem == 1 ) that1 = that2;
- this->value.nelem = that1->value.nelem;
- this->value.naxis = that1->value.naxis;
- for( i=0; i<that1->value.naxis; i++ )
- this->value.naxes[i] = that1->value.naxes[i];
- if ( Op == ACCUM && that1->type == BITSTR ) {
- /* ACCUM is rank-reducing on bit strings */
- this->value.nelem = 1;
- this->value.naxis = 1;
- this->value.naxes[0] = 1;
- }
- /* Both subnodes should be of same time */
- switch( that1->type ) {
- case BITSTR: this->DoOp = Do_BinOp_bit; break;
- case STRING: this->DoOp = Do_BinOp_str; break;
- case BOOLEAN: this->DoOp = Do_BinOp_log; break;
- case LONG: this->DoOp = Do_BinOp_lng; break;
- case DOUBLE: this->DoOp = Do_BinOp_dbl; break;
- }
- if( constant ) this->DoOp( this );
- }
- return( n );
- }
- static int New_Func( int returnType, funcOp Op, int nNodes,
- int Node1, int Node2, int Node3, int Node4,
- int Node5, int Node6, int Node7 )
- {
- return New_FuncSize(returnType, Op, nNodes,
- Node1, Node2, Node3, Node4,
- Node5, Node6, Node7, 0);
- }
- static int New_FuncSize( int returnType, funcOp Op, int nNodes,
- int Node1, int Node2, int Node3, int Node4,
- int Node5, int Node6, int Node7, int Size )
- /* If returnType==0 , use Node1's type and vector sizes as returnType, */
- /* else return a single value of type returnType */
- {
- Node *this, *that;
- int i,n,constant;
- if( Node1<0 || Node2<0 || Node3<0 || Node4<0 ||
- Node5<0 || Node6<0 || Node7<0 ) return(-1);
- n = Alloc_Node();
- if( n>=0 ) {
- this = gParse.Nodes + n;
- this->operation = (int)Op;
- this->DoOp = Do_Func;
- this->nSubNodes = nNodes;
- this->SubNodes[0] = Node1;
- this->SubNodes[1] = Node2;
- this->SubNodes[2] = Node3;
- this->SubNodes[3] = Node4;
- this->SubNodes[4] = Node5;
- this->SubNodes[5] = Node6;
- this->SubNodes[6] = Node7;
- i = constant = nNodes; /* Functions with zero params are not const */
- if (Op == poirnd_fct) constant = 0; /* Nor is Poisson deviate */
- while( i-- )
- constant = ( constant && OPER(this->SubNodes[i]) == CONST_OP );
-
- if( returnType ) {
- this->type = returnType;
- this->value.nelem = 1;
- this->value.naxis = 1;
- this->value.naxes[0] = 1;
- } else {
- that = gParse.Nodes + Node1;
- this->type = that->type;
- this->value.nelem = that->value.nelem;
- this->value.naxis = that->value.naxis;
- for( i=0; i<that->value.naxis; i++ )
- this->value.naxes[i] = that->value.naxes[i];
- }
- /* Force explicit size before evaluating */
- if (Size > 0) this->value.nelem = Size;
- if( constant ) this->DoOp( this );
- }
- return( n );
- }
- static int New_Deref( int Var, int nDim,
- int Dim1, int Dim2, int Dim3, int Dim4, int Dim5 )
- {
- int n, idx, constant;
- long elem=0;
- Node *this, *theVar, *theDim[MAXDIMS];
- if( Var<0 || Dim1<0 || Dim2<0 || Dim3<0 || Dim4<0 || Dim5<0 ) return(-1);
- theVar = gParse.Nodes + Var;
- if( theVar->operation==CONST_OP || theVar->value.nelem==1 ) {
- yyerror("Cannot index a scalar value");
- return(-1);
- }
- n = Alloc_Node();
- if( n>=0 ) {
- this = gParse.Nodes + n;
- this->nSubNodes = nDim+1;
- theVar = gParse.Nodes + (this->SubNodes[0]=Var);
- theDim[0] = gParse.Nodes + (this->SubNodes[1]=Dim1);
- theDim[1] = gParse.Nodes + (this->SubNodes[2]=Dim2);
- theDim[2] = gParse.Nodes + (this->SubNodes[3]=Dim3);
- theDim[3] = gParse.Nodes + (this->SubNodes[4]=Dim4);
- theDim[4] = gParse.Nodes + (this->SubNodes[5]=Dim5);
- constant = theVar->operation==CONST_OP;
- for( idx=0; idx<nDim; idx++ )
- constant = (constant && theDim[idx]->operation==CONST_OP);
- for( idx=0; idx<nDim; idx++ )
- if( theDim[idx]->value.nelem>1 ) {
- Free_Last_Node();
- yyerror("Cannot use an array as an index value");
- return(-1);
- } else if( theDim[idx]->type!=LONG ) {
- Free_Last_Node();
- yyerror("Index value must be an integer type");
- return(-1);
- }
- this->operation = '[';
- this->DoOp = Do_Deref;
- this->type = theVar->type;
- if( theVar->value.naxis == nDim ) { /* All dimensions specified */
- this->value.nelem = 1;
- this->value.naxis = 1;
- this->value.naxes[0] = 1;
- } else if( nDim==1 ) { /* Dereference only one dimension */
- elem=1;
- this->value.naxis = theVar->value.naxis-1;
- for( idx=0; idx<this->value.naxis; idx++ ) {
- elem *= ( this->value.naxes[idx] = theVar->value.naxes[idx] );
- }
- this->value.nelem = elem;
- } else {
- Free_Last_Node();
- yyerror("Must specify just one or all indices for vector");
- return(-1);
- }
- if( constant ) this->DoOp( this );
- }
- return(n);
- }
- extern int yyGetVariable( char *varName, YYSTYPE *varVal );
- static int New_GTI( char *fname, int Node1, char *start, char *stop )
- {
- fitsfile *fptr;
- Node *this, *that0, *that1;
- int type,i,n, startCol, stopCol, Node0;
- int hdutype, hdunum, evthdu, samefile, extvers, movetotype, tstat;
- char extname[100];
- long nrows;
- double timeZeroI[2], timeZeroF[2], dt, timeSpan;
- char xcol[20], xexpr[20];
- YYSTYPE colVal;
- if( Node1==-99 ) {
- type = yyGetVariable( "TIME", &colVal );
- if( type==COLUMN ) {
- Node1 = New_Column( (int)colVal.lng );
- } else {
- yyerror("Could not build TIME column for GTIFILTER");
- return(-1);
- }
- }
- Node1 = New_Unary( DOUBLE, 0, Node1 );
- Node0 = Alloc_Node(); /* This will hold the START/STOP times */
- if( Node1<0 || Node0<0 ) return(-1);
- /* Record current HDU number in case we need to move within this file */
- fptr = gParse.def_fptr;
- ffghdn( fptr, &evthdu );
- /* Look for TIMEZERO keywords in current extension */
- tstat = 0;
- if( ffgkyd( fptr, "TIMEZERO", timeZeroI, NULL, &tstat ) ) {
- tstat = 0;
- if( ffgkyd( fptr, "TIMEZERI", timeZeroI, NULL, &tstat ) ) {
- timeZeroI[0] = timeZeroF[0] = 0.0;
- } else if( ffgkyd( fptr, "TIMEZERF", timeZeroF, NULL, &tstat ) ) {
- timeZeroF[0] = 0.0;
- }
- } else {
- timeZeroF[0] = 0.0;
- }
- /* Resolve filename parameter */
- switch( fname[0] ) {
- case '\0':
- samefile = 1;
- hdunum = 1;
- break;
- case '[':
- samefile = 1;
- i = 1;
- while( fname[i] != '\0' && fname[i] != ']' ) i++;
- if( fname[i] ) {
- fname[i] = '\0';
- fname++;
- ffexts( fname, &hdunum, extname, &extvers, &movetotype,
- xcol, xexpr, &gParse.status );
- if( *extname ) {
- ffmnhd( fptr, movetotype, extname, extvers, &gParse.status );
- ffghdn( fptr, &hdunum );
- } else if( hdunum ) {
- ffmahd( fptr, ++hdunum, &hdutype, &gParse.status );
- } else if( !gParse.status ) {
- yyerror("Cannot use primary array for GTI filter");
- return( -1 );
- }
- } else {
- yyerror("File extension specifier lacks closing ']'");
- return( -1 );
- }
- break;
- case '+':
- samefile = 1;
- hdunum = atoi( fname ) + 1;
- if( hdunum>1 )
- ffmahd( fptr, hdunum, &hdutype, &gParse.status );
- else {
- yyerror("Cannot use primary array for GTI filter");
- return( -1 );
- }
- break;
- default:
- samefile = 0;
- if( ! ffopen( &fptr, fname, READONLY, &gParse.status ) )
- ffghdn( fptr, &hdunum );
- break;
- }
- if( gParse.status ) return(-1);
- /* If at primary, search for GTI extension */
- if( hdunum==1 ) {
- while( 1 ) {
- hdunum++;
- if( ffmahd( fptr, hdunum, &hdutype, &gParse.status ) ) break;
- if( hdutype==IMAGE_HDU ) continue;
- tstat = 0;
- if( ffgkys( fptr, "EXTNAME", extname, NULL, &tstat ) ) continue;
- ffupch( extname );
- if( strstr( extname, "GTI" ) ) break;
- }
- if( gParse.status ) {
- if( gParse.status==END_OF_FILE )
- yyerror("GTI extension not found in this file");
- return(-1);
- }
- }
- /* Locate START/STOP Columns */
- ffgcno( fptr, CASEINSEN, start, &startCol, &gParse.status );
- ffgcno( fptr, CASEINSEN, stop, &stopCol, &gParse.status );
- if( gParse.status ) return(-1);
- /* Look for TIMEZERO keywords in GTI extension */
- tstat = 0;
- if( ffgkyd( fptr, "TIMEZERO", timeZeroI+1, NULL, &tstat ) ) {
- tstat = 0;
- if( ffgkyd( fptr, "TIMEZERI", timeZeroI+1, NULL, &tstat ) ) {
- timeZeroI[1] = timeZeroF[1] = 0.0;
- } else if( ffgkyd( fptr, "TIMEZERF", timeZeroF+1, NULL, &tstat ) ) {
- timeZeroF[1] = 0.0;
- }
- } else {
- timeZeroF[1] = 0.0;
- }
- n = Alloc_Node();
- if( n >= 0 ) {
- this = gParse.Nodes + n;
- this->nSubNodes = 2;
- this->SubNodes[1] = Node1;
- this->operation = (int)gtifilt_fct;
- this->DoOp = Do_GTI;
- this->type = BOOLEAN;
- that1 = gParse.Nodes + Node1;
- this->value.nelem = that1->value.nelem;
- this->value.naxis = that1->value.naxis;
- for( i=0; i < that1->value.naxis; i++ )
- this->value.naxes[i] = that1->value.naxes[i];
- /* Init START/STOP node to be treated as a "constant" */
- this->SubNodes[0] = Node0;
- that0 = gParse.Nodes + Node0;
- that0->operation = CONST_OP;
- that0->DoOp = NULL;
- that0->value.data.ptr= NULL;
- /* Read in START/STOP times */
- if( ffgkyj( fptr, "NAXIS2", &nrows, NULL, &gParse.status ) )
- return(-1);
- that0->value.nelem = nrows;
- if( nrows ) {
- that0->value.data.dblptr = (double*)malloc( 2*nrows*sizeof(double) );
- if( !that0->value.data.dblptr ) {
- gParse.status = MEMORY_ALLOCATION;
- return(-1);
- }
-
- ffgcvd( fptr, startCol, 1L, 1L, nrows, 0.0,
- that0->value.data.dblptr, &i, &gParse.status );
- ffgcvd( fptr, stopCol, 1L, 1L, nrows, 0.0,
- that0->value.data.dblptr+nrows, &i, &gParse.status );
- if( gParse.status ) {
- free( that0->value.data.dblptr );
- return(-1);
- }
- /* Test for fully time-ordered GTI... both START && STOP */
- that0->type = 1; /* Assume yes */
- i = nrows;
- while( --i )
- if( that0->value.data.dblptr[i-1]
- >= that0->value.data.dblptr[i]
- || that0->value.data.dblptr[i-1+nrows]
- >= that0->value.data.dblptr[i+nrows] ) {
- that0->type = 0;
- break;
- }
-
- /* Handle TIMEZERO offset, if any */
-
- dt = (timeZeroI[1] - timeZeroI[0]) + (timeZeroF[1] - timeZeroF[0]);
- timeSpan = that0->value.data.dblptr[nrows+nrows-1]
- - that0->value.data.dblptr[0];
-
- if( fabs( dt / timeSpan ) > 1e-12 ) {
- for( i=0; i<(nrows+nrows); i++ )
- that0->value.data.dblptr[i] += dt;
- }
- }
- if( OPER(Node1)==CONST_OP )
- this->DoOp( this );
- }
- if( samefile )
- ffmahd( fptr, evthdu, &hdutype, &gParse.status );
- else
- ffclos( fptr, &gParse.status );
- return( n );
- }
- static int New_REG( char *fname, int NodeX, int NodeY, char *colNames )
- {
- Node *this, *that0;
- int type, n, Node0;
- int Xcol, Ycol, tstat;
- WCSdata wcs;
- SAORegion *Rgn;
- char *cX, *cY;
- YYSTYPE colVal;
- if( NodeX==-99 ) {
- type = yyGetVariable( "X", &colVal );
- if( type==COLUMN ) {
- NodeX = New_Column( (int)colVal.lng );
- } else {
- yyerror("Could not build X column for REGFILTER");
- return(-1);
- }
- }
- if( NodeY==-99 ) {
- type = yyGetVariable( "Y", &colVal );
- if( type==COLUMN ) {
- NodeY = New_Column( (int)colVal.lng );
- } else {
- yyerror("Could not build Y column for REGFILTER");
- return(-1);
- }
- }
- NodeX = New_Unary( DOUBLE, 0, NodeX );
- NodeY = New_Unary( DOUBLE, 0, NodeY );
- Node0 = Alloc_Node(); /* This will hold the Region Data */
- if( NodeX<0 || NodeY<0 || Node0<0 ) return(-1);
- if( ! (Test_Dims( NodeX, NodeY ) ) ) {
- yyerror("Dimensions of REGFILTER arguments are not compatible");
- return (-1);
- }
- n = Alloc_Node();
- if( n >= 0 ) {
- this = gParse.Nodes + n;
- this->nSubNodes = 3;
- this->SubNodes[0] = Node0;
- this->SubNodes[1] = NodeX;
- this->SubNodes[2] = NodeY;
- this->operation = (int)regfilt_fct;
- this->DoOp = Do_REG;
- this->type = BOOLEAN;
- this->value.nelem = 1;
- this->value.naxis = 1;
- this->value.naxes[0] = 1;
-
- Copy_Dims(n, NodeX);
- if( SIZE(NodeX)<SIZE(NodeY) ) Copy_Dims(n, NodeY);
- /* Init Region node to be treated as a "constant" */
- that0 = gParse.Nodes + Node0;
- that0->operation = CONST_OP;
- that0->DoOp = NULL;
- /* Identify what columns to use for WCS information */
- Xcol = Ycol = 0;
- if( *colNames ) {
- /* Use the column names in this string for WCS info */
- while( *colNames==' ' ) colNames++;
- cX = cY = colNames;
- while( *cY && *cY!=' ' && *cY!=',' ) cY++;
- if( *cY )
- *(cY++) = '\0';
- while( *cY==' ' ) cY++;
- if( !*cY ) {
- yyerror("Could not extract valid pair of column names from REGFILTER");
- Free_Last_Node();
- return( -1 );
- }
- fits_get_colnum( gParse.def_fptr, CASEINSEN, cX, &Xcol,
- &gParse.status );
- fits_get_colnum( gParse.def_fptr, CASEINSEN, cY, &Ycol,
- &gParse.status );
- if( gParse.status ) {
- yyerror("Could not locate columns indicated for WCS info");
- Free_Last_Node();
- return( -1 );
- }
- } else {
- /* Try to find columns used in X/Y expressions */
- Xcol = Locate_Col( gParse.Nodes + NodeX );
- Ycol = Locate_Col( gParse.Nodes + NodeY );
- if( Xcol<0 || Ycol<0 ) {
- yyerror("Found multiple X/Y column references in REGFILTER");
- Free_Last_Node();
- return( -1 );
- }
- }
- /* Now, get the WCS info, if it exists, from the indicated columns */
- wcs.exists = 0;
- if( Xcol>0 && Ycol>0 ) {
- tstat = 0;
- ffgtcs( gParse.def_fptr, Xcol, Ycol,
- &wcs.xrefval, &wcs.yrefval,
- &wcs.xrefpix, &wcs.yrefpix,
- &wcs.xinc, &wcs.yinc,
- &wcs.rot, wcs.type,
- &tstat );
- if( tstat==NO_WCS_KEY ) {
- wcs.exists = 0;
- } else if( tstat ) {
- gParse.status = tstat;
- Free_Last_Node();
- return( -1 );
- } else {
- wcs.exists = 1;
- }
- }
- /* Read in Region file */
- fits_read_rgnfile( fname, &wcs, &Rgn, &gParse.status );
- if( gParse.status ) {
- Free_Last_Node();
- return( -1 );
- }
- that0->value.data.ptr = Rgn;
- if( OPER(NodeX)==CONST_OP && OPER(NodeY)==CONST_OP )
- this->DoOp( this );
- }
- return( n );
- }
- static int New_Vector( int subNode )
- {
- Node *this, *that;
- int n;
- n = Alloc_Node();
- if( n >= 0 ) {
- this = gParse.Nodes + n;
- that = gParse.Nodes + subNode;
- this->type = that->type;
- this->nSubNodes = 1;
- this->SubNodes[0] = subNode;
- this->operation = '{';
- this->DoOp = Do_Vector;
- }
- return( n );
- }
- static int Close_Vec( int vecNode )
- {
- Node *this;
- int n, nelem=0;
- this = gParse.Nodes + vecNode;
- for( n=0; n < this->nSubNodes; n++ ) {
- if( TYPE( this->SubNodes[n] ) != this->type ) {
- this->SubNodes[n] = New_Unary( this->type, 0, this->SubNodes[n] );
- if( this->SubNodes[n]<0 ) return(-1);
- }
- nelem += SIZE(this->SubNodes[n]);
- }
- this->value.naxis = 1;
- this->value.nelem = nelem;
- this->value.naxes[0] = nelem;
- return( vecNode );
- }
- static int Locate_Col( Node *this )
- /* Locate the TABLE column number of any columns in "this" calculation. */
- /* Return ZERO if none found, or negative if more than 1 found. */
- {
- Node *that;
- int i, col=0, newCol, nfound=0;
-
- if( this->nSubNodes==0
- && this->operation<=0 && this->operation!=CONST_OP )
- return gParse.colData[ - this->operation].colnum;
- for( i=0; i<this->nSubNodes; i++ ) {
- that = gParse.Nodes + this->SubNodes[i];
- if( that->operation>0 ) {
- newCol = Locate_Col( that );
- if( newCol<=0 ) {
- nfound += -newCol;
- } else {
- if( !nfound ) {
- col = newCol;
- nfound++;
- } else if( col != newCol ) {
- nfound++;
- }
- }
- } else if( that->operation!=CONST_OP ) {
- /* Found a Column */
- newCol = gParse.colData[- that->operation].colnum;
- if( !nfound ) {
- col = newCol;
- nfound++;
- } else if( col != newCol ) {
- nfound++;
- }
- }
- }
- if( nfound!=1 )
- return( - nfound );
- else
- return( col );
- }
- static int Test_Dims( int Node1, int Node2 )
- {
- Node *that1, *that2;
- int valid, i;
- if( Node1<0 || Node2<0 ) return(0);
- that1 = gParse.Nodes + Node1;
- that2 = gParse.Nodes + Node2;
- if( that1->value.nelem==1 || that2->value.nelem==1 )
- valid = 1;
- else if( that1->type==that2->type
- && that1->value.nelem==that2->value.nelem
- && that1->value.naxis==that2->value.naxis ) {
- valid = 1;
- for( i=0; i<that1->value.naxis; i++ ) {
- if( that1->value.naxes[i]!=that2->value.naxes[i] )
- valid = 0;
- }
- } else
- valid = 0;
- return( valid );
- }
- static void Copy_Dims( int Node1, int Node2 )
- {
- Node *that1, *that2;
- int i;
- if( Node1<0 || Node2<0 ) return;
- that1 = gParse.Nodes + Node1;
- that2 = gParse.Nodes + Node2;
- that1->value.nelem = that2->value.nelem;
- that1->value.naxis = that2->value.naxis;
- for( i=0; i<that2->value.naxis; i++ )
- that1->value.naxes[i] = that2->value.naxes[i];
- }
- /********************************************************************/
- /* Routines for actually evaluating the expression start here */
- /********************************************************************/
- void Evaluate_Parser( long firstRow, long nRows )
- /***********************************************************************/
- /* Reset the parser for processing another batch of data... */
- /* firstRow: Row number of the first element to evaluate */
- /* nRows: Number of rows to be processed */
- /* Initialize each COLUMN node so that its UNDEF and DATA pointers */
- /* point to the appropriate column arrays. */
- /* Finally, call Evaluate_Node for final node. */
- /***********************************************************************/
- {
- int i, column;
- long offset, rowOffset;
- gParse.firstRow = firstRow;
- gParse.nRows = nRows;
- /* Reset Column Nodes' pointers to point to right data and UNDEF arrays */
- rowOffset = firstRow - gParse.firstDataRow;
- for( i=0; i<gParse.nNodes; i++ ) {
- if( OPER(i) > 0 || OPER(i) == CONST_OP ) continue;
- column = -OPER(i);
- offset = gParse.varData[column].nelem * rowOffset;
- gParse.Nodes[i].value.undef = gParse.varData[column].undef + offset;
- switch( gParse.Nodes[i].type ) {
- case BITSTR:
- gParse.Nodes[i].value.data.strptr =
- (char**)gParse.varData[column].data + rowOffset;
- gParse.Nodes[i].value.undef = NULL;
- break;
- case STRING:
- gParse.Nodes[i].value.data.strptr =
- (char**)gParse.varData[column].data + rowOffset;
- gParse.Nodes[i].value.undef = gParse.varData[column].undef + rowOffset;
- break;
- case BOOLEAN:
- gParse.Nodes[i].value.data.logptr =
- (char*)gParse.varData[column].data + offset;
- break;
- case LONG:
- gParse.Nodes[i].value.data.lngptr =
- (long*)gParse.varData[column].data + offset;
- break;
- case DOUBLE:
- gParse.Nodes[i].value.data.dblptr =
- (double*)gParse.varData[column].data + offset;
- break;
- }
- }
- Evaluate_Node( gParse.resultNode );
- }
- static void Evaluate_Node( int thisNode )
- /**********************************************************************/
- /* Recursively evaluate thisNode's subNodes, then call one of the */
- /* Do_<Action> functions pointed to by thisNode's DoOp element. */
- /**********************************************************************/
- {
- Node *this;
- int i;
-
- if( gParse.status ) return;
- this = gParse.Nodes + thisNode;
- if( this->operation>0 ) { /* <=0 indicate constants and columns */
- i = this->nSubNodes;
- while( i-- ) {
- Evaluate_Node( this->SubNodes[i] );
- if( gParse.status ) return;
- }
- this->DoOp( this );
- }
- }
- static void Allocate_Ptrs( Node *this )
- {
- long elem, row, size;
- if( this->type==BITSTR || this->type==STRING ) {
- this->value.data.strptr = (char**)malloc( gParse.nRows
- * sizeof(char*) );
- if( this->value.data.strptr ) {
- this->value.data.strptr[0] = (char*)malloc( gParse.nRows
- * (this->value.nelem+2)
- * sizeof(char) );
- if( this->value.data.strptr[0] ) {
- row = 0;
- while( (++row)<gParse.nRows ) {
- this->value.data.strptr[row] =
- this->value.data.strptr[row-1] + this->value.nelem+1;
- }
- if( this->type==STRING ) {
- this->value.undef = this->value.data.strptr[row-1]
- + this->value.nelem+1;
- } else {
- this->value.undef = NULL; /* BITSTRs don't use undef array */
- }
- } else {
- gParse.status = MEMORY_ALLOCATION;
- free( this->value.data.strptr );
- }
- } else {
- gParse.status = MEMORY_ALLOCATION;
- }
- } else {
- elem = this->value.nelem * gParse.nRows;
- switch( this->type ) {
- case DOUBLE: size = sizeof( double ); break;
- case LONG: size = sizeof( long ); break;
- case BOOLEAN: size = sizeof( char ); break;
- default: size = 1; break;
- }
- this->value.data.ptr = calloc(size+1, elem);
- if( this->value.data.ptr==NULL ) {
- gParse.status = MEMORY_ALLOCATION;
- } else {
- this->value.undef = (char *)this->value.data.ptr + elem*size;
- }
- }
- }
- static void Do_Unary( Node *this )
- {
- Node *that;
- long elem;
- that = gParse.Nodes + this->SubNodes[0];
- if( that->operation==CONST_OP ) { /* Operating on a constant! */
- switch( this->operation ) {
- case DOUBLE:
- case FLTCAST:
- if( that->type==LONG )
- this->value.data.dbl = (double)that->value.data.lng;
- else if( that->type==BOOLEAN )
- this->value.data.dbl = ( that->value.data.log ? 1.0 : 0.0 );
- break;
- case LONG:
- case INTCAST:
- if( that->type==DOUBLE )
- this->value.data.lng = (long)that->value.data.dbl;
- else if( that->type==BOOLEAN )
- this->value.data.lng = ( that->value.data.log ? 1L : 0L );
- break;
- case BOOLEAN:
- if( that->type==DOUBLE )
- this->value.data.log = ( that->value.data.dbl != 0.0 );
- else if( that->type==LONG )
- this->value.data.log = ( that->value.data.lng != 0L );
- break;
- case UMINUS:
- if( that->type==DOUBLE )
- this->value.data.dbl = - that->value.data.dbl;
- else if( that->type==LONG )
- this->value.data.lng = - that->value.data.lng;
- break;
- case NOT:
- if( that->type==BOOLEAN )
- this->value.data.log = ( ! that->value.data.log );
- else if( that->type==BITSTR )
- bitnot( this->value.data.str, that->value.data.str );
- break;
- }
- this->operation = CONST_OP;
- } else {
- Allocate_Ptrs( this );
- if( !gParse.status ) {
- if( this->type!=BITSTR ) {
- elem = gParse.nRows;
- if( this->type!=STRING )
- elem *= this->value.nelem;
- while( elem-- )
- this->value.undef[elem] = that->value.undef[elem];
- }
- elem = gParse.nRows * this->value.nelem;
- switch( this->operation ) {
- case BOOLEAN:
- if( that->type==DOUBLE )
- while( elem-- )
- this->value.data.logptr[elem] =
- ( that->value.data.dblptr[elem] != 0.0 );
- else if( that->type==LONG )
- while( elem-- )
- this->value.data.logptr[elem] =
- ( that->value.data.lngptr[elem] != 0L );
- break;
- case DOUBLE:
- case FLTCAST:
- if( that->type==LONG )
- while( elem-- )
- this->value.data.dblptr[elem] =
- (double)that->value.data.lngptr[elem];
- else if( that->type==BOOLEAN )
- while( elem-- )
- this->value.data.dblptr[elem] =
- ( that->value.data.logptr[elem] ? 1.0 : 0.0 );
- break;
- case LONG:
- case INTCAST:
- if( that->type==DOUBLE )
- while( elem-- )
- this->value.data.lngptr[elem] =
- (long)that->value.data.dblptr[elem];
- else if( that->type==BOOLEAN )
- while( elem-- )
- this->value.data.lngptr[elem] =
- ( that->value.data.logptr[elem] ? 1L : 0L );
- break;
- case UMINUS:
- if( that->type==DOUBLE ) {
- while( elem-- )
- this->value.data.dblptr[elem] =
- - that->value.data.dblptr[elem];
- } else if( that->type==LONG ) {
- while( elem-- )
- this->value.data.lngptr[elem] =
- - that->value.data.lngptr[elem];
- }
- break;
- case NOT:
- if( that->type==BOOLEAN ) {
- while( elem-- )
- this->value.data.logptr[elem] =
- ( ! that->value.data.logptr[elem] );
- } else if( that->type==BITSTR ) {
- elem = gParse.nRows;
- while( elem-- )
- bitnot( this->value.data.strptr[elem],
- that->value.data.strptr[elem] );
- }
- break;
- }
- }
- }
- if( that->operation>0 ) {
- free( that->value.data.ptr );
- }
- }
- static void Do_Offset( Node *this )
- {
- Node *col;
- long fRow, nRowOverlap, nRowReload, rowOffset;
- long nelem, elem, offset, nRealElem;
- int status;
- col = gParse.Nodes + this->SubNodes[0];
- rowOffset = gParse.Nodes[ this->SubNodes[1] ].value.data.lng;
- Allocate_Ptrs( this );
- fRow = gParse.firstRow + rowOffset;
- if( this->type==STRING || this->type==BITSTR )
- nRealElem = 1;
- else
- nRealElem = this->value.nelem;
- nelem = nRealElem;
- if( fRow < gParse.firstDataRow ) {
- /* Must fill in data at start of array */
- nRowReload = gParse.firstDataRow - fRow;
- if( nRowReload > gParse.nRows ) nRowReload = gParse.nRows;
- nRowOverlap = gParse.nRows - nRowReload;
- offset = 0;
- /* NULLify any values falling out of bounds */
- while( fRow<1 && nRowReload>0 ) {
- if( this->type == BITSTR ) {
- nelem = this->value.nelem;
- this->value.data.strptr[offset][ nelem ] = '\0';
- while( nelem-- ) this->value.data.strptr[offset][nelem] = '0';
- offset++;
- } else {
- while( nelem-- )
- this->value.undef[offset++] = 1;
- }
- nelem = nRealElem;
- fRow++;
- nRowReload--;
- }
- } else if( fRow + gParse.nRows > gParse.firstDataRow + gParse.nDataRows ) {
- /* Must fill in data at end of array */
- nRowReload = (fRow+gParse.nRows) - (gParse.firstDataRow+gParse.nDataRows);
- if( nRowReload>gParse.nRows ) {
- nRowReload = gParse.nRows;
- } else {
- fRow = gParse.firstDataRow + gParse.nDataRows;
- }
- nRowOverlap = gParse.nRows - nRowReload;
- offset = nRowOverlap * nelem;
- /* NULLify any values falling out of bounds */
- elem = gParse.nRows * nelem;
- while( fRow+nRowReload>gParse.totalRows && nRowReload>0 ) {
- if( this->type == BITSTR ) {
- nelem = this->value.nelem;
- elem--;
- this->value.data.strptr[elem][ nelem ] = '\0';
- while( nelem-- ) this->value.data.strptr[elem][nelem] = '0';
- } else {
- while( nelem-- )
- this->value.undef[--elem] = 1;
- }
- nelem = nRealElem;
- nRowReload--;
- }
- } else {
- nRowReload = 0;
- nRowOverlap = gParse.nRows;
- offset = 0;
- }
- if( nRowReload>0 ) {
- switch( this->type ) {
- case BITSTR:
- case STRING:
- status = (*gParse.loadData)( -col->operation, fRow, nRowReload,
- this->value.data.strptr+offset,
- this->value.undef+offset );
- break;
- case BOOLEAN:
- status = (*gParse.loadData)( -col->operation, fRow, nRowReload,
- this->value.data.logptr+offset,
- this->value.undef+offset );
- break;
- case LONG:
- status = (*gParse.loadData)( -col->operation, fRow, nRowReload,
- this->value.data.lngptr+offset,
- this->value.undef+offset );
- break;
- case DOUBLE:
- status = (*gParse.loadData)( -col->operation, fRow, nRowReload,
- this->value.data.dblptr+offset,
- this->value.undef+offset );
- break;
- }
- }
- /* Now copy over the overlapping region, if any */
- if( nRowOverlap <= 0 ) return;
- if( rowOffset>0 )
- elem = nRowOverlap * nelem;
- else
- elem = gParse.nRows * nelem;
- offset = nelem * rowOffset;
- while( nRowOverlap-- && !gParse.status ) {
- while( nelem-- && !gParse.status ) {
- elem--;
- if( this->type != BITSTR )
- this->value.undef[elem] = col->value.undef[elem+offset];
- switch( this->type ) {
- case BITSTR:
- strcpy( this->value.data.strptr[elem ],
- col->value.data.strptr[elem+offset] );
- break;
- case STRING:
- strcpy( this->value.data.strptr[elem ],
- col->value.data.strptr[elem+offset] );
- break;
- case BOOLEAN:
- this->value.data.logptr[elem] = col->value.data.logptr[elem+offset];
- break;
- case LONG:
- this->value.data.lngptr[elem] = col->value.data.lngptr[elem+offset];
- break;
- case DOUBLE:
- this->value.data.dblptr[elem] = col->value.data.dblptr[elem+offset];
- break;
- }
- }
- nelem = nRealElem;
- }
- }
- static void Do_BinOp_bit( Node *this )
- {
- Node *that1, *that2;
- char *sptr1=NULL, *sptr2=NULL;
- int const1, const2;
- long rows;
- that1 = gParse.Nodes + this->SubNodes[0];
- that2 = gParse.Nodes + this->SubNodes[1];
- const1 = ( that1->operation==CONST_OP );
- const2 = ( that2->operation==CONST_OP );
- sptr1 = ( const1 ? that1->value.data.str : NULL );
- sptr2 = ( const2 ? that2->value.data.str : NULL );
- if( const1 && const2 ) {
- switch( this->operation ) {
- case NE:
- this->value.data.log = !bitcmp( sptr1, sptr2 );
- break;
- case EQ:
- this->value.data.log = bitcmp( sptr1, sptr2 );
- break;
- case GT:
- case LT:
- case LTE:
- case GTE:
- this->value.data.log = bitlgte( sptr1, this->operation, sptr2 );
- break;
- case '|':
- bitor( this->value.data.str, sptr1, sptr2 );
- break;
- case '&':
- bitand( this->value.data.str, sptr1, sptr2 );
- break;
- case '+':
- strcpy( this->value.data.str, sptr1 );
- strcat( this->value.data.str, sptr2 );
- break;
- case ACCUM:
- this->value.data.lng = 0;
- while( *sptr1 ) {
- if ( *sptr1 == '1' ) this->value.data.lng ++;
- sptr1 ++;
- }
- break;
-
- }
- this->operation = CONST_OP;
- } else {
- Allocate_Ptrs( this );
- if( !gParse.status ) {
- rows = gParse.nRows;
- switch( this->operation ) {
- /* BITSTR comparisons */
- case NE:
- case EQ:
- case GT:
- case LT:
- case LTE:
- case GTE:
- while( rows-- ) {
- if( !const1 )
- sptr1 = that1->value.data.strptr[rows];
- if( !const2 )
- sptr2 = that2->value.data.strptr[rows];
- switch( this->operation ) {
- case NE: this->value.data.logptr[rows] =
- !bitcmp( sptr1, sptr2 );
- break;
- case EQ: this->value.data.logptr[rows] =
- bitcmp( sptr1, sptr2 );
- break;
- case GT:
- case LT:
- case LTE:
- case GTE: this->value.data.logptr[rows] =
- bitlgte( sptr1, this->operation, sptr2 );
- break;
- }
- this->value.undef[rows] = 0;
- }
- break;
-
- /* BITSTR AND/ORs ... no UNDEFS in or out */
-
- case '|':
- case '&':
- case '+':
- while( rows-- ) {
- if( !const1 )
- sptr1 = that1->value.data.strptr[rows];
- if( !const2 )
- sptr2 = that2->value.data.strptr[rows];
- if( this->operation=='|' )
- bitor( this->value.data.strptr[rows], sptr1, sptr2 );
- else if( this->operation=='&' )
- bitand( this->value.data.strptr[rows], sptr1, sptr2 );
- else {
- strcpy( this->value.data.strptr[rows], sptr1 );
- strcat( this->value.data.strptr[rows], sptr2 );
- }
- }
- break;
- /* Accumulate 1 bits */
- case ACCUM:
- {
- long i, previous, curr;
- previous = that2->value.data.lng;
-
- /* Cumulative sum of this chunk */
- for (i=0; i<rows; i++) {
- sptr1 = that1->value.data.strptr[i];
- for (curr = 0; *sptr1; sptr1 ++) {
- if ( *sptr1 == '1' ) curr ++;
- }
- previous += curr;
- this->value.data.lngptr[i] = previous;
- this->value.undef[i] = 0;
- }
-
- /* Store final cumulant for next pass */
- that2->value.data.lng = previous;
- }
- }
- }
- }
- if( that1->operation>0 ) {
- free( that1->value.data.strptr[0] );
- free( that1->value.data.strptr );
- }
- if( that2->operation>0 ) {
- free( that2->value.data.strptr[0] );
- free( that2->value.data.strptr );
- }
- }
- static void Do_BinOp_str( Node *this )
- {
- Node *that1, *that2;
- char *sptr1, *sptr2, null1=0, null2=0;
- int const1, const2, val;
- long rows;
- that1 = gParse.Nodes + this->SubNodes[0];
- that2 = gParse.Nodes + this->SubNodes[1];
- const1 = ( that1->operation==CONST_OP );
- const2 = ( that2->operation==CONST_OP );
- sptr1 = ( const1 ? that1->value.data.str : NULL );
- sptr2 = ( const2 ? that2->value.data.str : NULL );
- if( const1 && const2 ) { /* Result is a constant */
- switch( this->operation ) {
- /* Compare Strings */
- case NE:
- case EQ:
- val = ( FSTRCMP( sptr1, sptr2 ) == 0 );
- this->value.data.log = ( this->operation==EQ ? val : !val );
- break;
- case GT:
- this->value.data.log = ( FSTRCMP( sptr1, sptr2 ) > 0 );
- break;
- case LT:
- this->value.data.log = ( FSTRCMP( sptr1, sptr2 ) < 0 );
- break;
- case GTE:
- this->value.data.log = ( FSTRCMP( sptr1, sptr2 ) >= 0 );
- break;
- case LTE:
- this->value.data.log = ( FSTRCMP( sptr1, sptr2 ) <= 0 );
- break;
- /* Concat Strings */
- case '+':
- strcpy( this->value.data.str, sptr1 );
- strcat( this->value.data.str, sptr2 );
- break;
- }
- this->operation = CONST_OP;
- } else { /* Not a constant */
- Allocate_Ptrs( this );
- if( !gParse.status ) {
- rows = gParse.nRows;
- switch( this->operation ) {
- /* Compare Strings */
- case NE:
- case EQ:
- while( rows-- ) {
- if( !const1 ) null1 = that1->value.undef[rows];
- if( !const2 ) null2 = that2->value.undef[rows];
- this->value.undef[rows] = (null1 || null2);
- if( ! this->value.undef[rows] ) {
- if( !const1 ) sptr1 = that1->value.data.strptr[rows];
- if( !const2 ) sptr2 = that2->value.data.strptr[rows];
- val = ( FSTRCMP( sptr1, sptr2 ) == 0 );
- this->value.data.logptr[rows] =
- ( this->operation==EQ ? val : !val );
- }
- }
- break;
-
- case GT:
- case LT:
- while( rows-- ) {
- if( !const1 ) null1 = that1->value.undef[rows];
- if( !const2 ) null2 = that2->value.undef[rows];
- this->value.undef[rows] = (null1 || null2);
- if( ! this->value.undef[rows] ) {
- if( !const1 ) sptr1 = that1->value.data.strptr[rows];
- if( !const2 ) sptr2 = that2->value.data.strptr[rows];
- val = ( FSTRCMP( sptr1, sptr2 ) );
- this->value.data.logptr[rows] =
- ( this->operation==GT ? val>0 : val<0 );
- }
- }
- break;
- case GTE:
- case LTE:
- while( rows-- ) {
- if( !const1 ) null1 = that1->value.undef[rows];
- if( !const2 ) null2 = that2->value.undef[rows];
- this->value.undef[rows] = (null1 || null2);
- if( ! this->value.undef[rows] ) {
- if( !const1 ) sptr1 = that1->value.data.strptr[rows];
- if( !const2 ) sptr2 = that2->value.data.strptr[rows];
- val = ( FSTRCMP( sptr1, sptr2 ) );
- this->value.data.logptr[rows] =
- ( this->operation==GTE ? val>=0 : val<=0 );
- }
- }
- break;
- /* Concat Strings */
-
- case '+':
- while( rows-- ) {
- if( !const1 ) null1 = that1->value.undef[rows];
- if( !const2 ) null2 = that2->value.undef[rows];
- this->value.undef[rows] = (null1 || null2);
- if( ! this->value.undef[rows] ) {
- if( !const1 ) sptr1 = that1->value.data.strptr[rows];
- if( !const2 ) sptr2 = that2->value.data.strptr[rows];
- strcpy( this->value.data.strptr[rows], sptr1 );
- strcat( this->value.data.strptr[rows], sptr2 );
- }
- }
- break;
- }
- }
- }
- if( that1->operation>0 ) {
- free( that1->value.data.strptr[0] );
- free( that1->value.data.strptr );
- }
- if( that2->operation>0 ) {
- free( that2->value.data.strptr[0] );
- free( that2->value.data.strptr );
- }
- }
- static void Do_BinOp_log( Node *this )
- {
- Node *that1, *that2;
- int vector1, vector2;
- char val1=0, val2=0, null1=0, null2=0;
- long rows, nelem, elem;
- that1 = gParse.Nodes + this->SubNodes[0];
- that2 = gParse.Nodes + this->SubNodes[1];
- vector1 = ( that1->operation!=CONST_OP );
- if( vector1 )
- vector1 = that1->value.nelem;
- else {
- val1 = that1->value.data.log;
- }
- vector2 = ( that2->operation!=CONST_OP );
- if( vector2 )
- vector2 = that2->value.nelem;
- else {
- val2 = that2->value.data.log;
- }
- if( !vector1 && !vector2 ) { /* Result is a constant */
- switch( this->operation ) {
- case OR:
- this->value.data.log = (val1 || val2);
- break;
- case AND:
- this->value.data.log = (val1 && val2);
- break;
- case EQ:
- this->value.data.log = ( (val1 && val2) || (!val1 && !val2) );
- break;
- case NE:
- this->value.data.log = ( (val1 && !val2) || (!val1 && val2) );
- break;
- case ACCUM:
- this->value.data.lng = val1;
- break;
- }
- this->operation=CONST_OP;
- } else if (this->operation == ACCUM) {
- long i, previous, curr;
- rows = gParse.nRows;
- nelem = this->value.nelem;
- elem = this->value.nelem * rows;
-
- Allocate_Ptrs( this );
-
- if( !gParse.status ) {
- previous = that2->value.data.lng;
-
- /* Cumulative sum of this chunk */
- for (i=0; i<elem; i++) {
- if (!that1->value.undef[i]) {
- curr = that1->value.data.logptr[i];
- previous += curr;
- }
- this->value.data.lngptr[i] = previous;
- this->value.undef[i] = 0;
- }
-
- /* Store final cumulant for next pass */
- that2->value.data.lng = previous;
- }
-
- } else {
- rows = gParse.nRows;
- nelem = this->value.nelem;
- elem = this->value.nelem * rows;
- Allocate_Ptrs( this );
- if( !gParse.status ) {
-
- if (this->operation == ACCUM) {
- long i, previous, curr;
-
- previous = that2->value.data.lng;
-
- /* Cumulative sum of this chunk */
- for (i=0; i<elem; i++) {
- if (!that1->value.undef[i]) {
- curr = that1->value.data.logptr[i];
- previous += curr;
- }
- this->value.data.lngptr[i] = previous;
- this->value.undef[i] = 0;
- }
-
- /* Store final cumulant for next pass */
- that2->value.data.lng = previous;
- }
-
- while( rows-- ) {
- while( nelem-- ) {
- elem--;
- if( vector1>1 ) {
- val1 = that1->value.data.logptr[elem];
- null1 = that1->value.undef[elem];
- } else if( vector1 ) {
- val1 = that1->value.data.logptr[rows];
- null1 = that1->value.undef[rows];
- }
- if( vector2>1 ) {
- val2 = that2->value.data.logptr[elem];
- null2 = that2->value.undef[elem];
- } else if( vector2 ) {
- val2 = that2->value.data.logptr[rows];
- null2 = that2->value.undef[rows];
- }
- this->value.undef[elem] = (null1 || null2);
- switch( this->operation ) {
- case OR:
- /* This is more complicated than others to suppress UNDEFs */
- /* in those cases where the other argument is DEF && TRUE */
- if( !null1 && !null2 ) {
- this->value.data.logptr[elem] = (val1 || val2);
- } else if( (null1 && !null2 && val2)
- || ( !null1 && null2 && val1 ) ) {
- this->value.data.logptr[elem] = 1;
- this->value.undef[elem] = 0;
- }
- break;
- case AND:
- /* This is more complicated than others to suppress UNDEFs */
- /* in those cases where the other argument is DEF && FALSE */
- if( !null1 && !null2 ) {
- this->value.data.logptr[elem] = (val1 && val2);
- } else if( (null1 && !null2 && !val2)
- || ( !null1 && null2 && !val1 ) ) {
- this->value.data.logptr[elem] = 0;
- this->value.undef[elem] = 0;
- }
- break;
- case EQ:
- this->value.data.logptr[elem] =
- ( (val1 && val2) || (!val1 && !val2) );
- break;
- case NE:
- this->value.data.logptr[elem] =
- ( (val1 && !val2) || (!val1 && val2) );
- break;
- }
- }
- nelem = this->value.nelem;
- }
- }
- }
- if( that1->operation>0 ) {
- free( that1->value.data.ptr );
- }
- if( that2->operation>0 ) {
- free( that2->value.data.ptr );
- }
- }
- static void Do_BinOp_lng( Node *this )
- {
- Node *that1, *that2;
- int vector1, vector2;
- long val1=0, val2=0;
- char null1=0, null2=0;
- long rows, nelem, elem;
- that1 = gParse.Nodes + this->SubNodes[0];
- that2 = gParse.Nodes + this->SubNodes[1];
- vector1 = ( that1->operation!=CONST_OP );
- if( vector1 )
- vector1 = that1->value.nelem;
- else {
- val1 = that1->value.data.lng;
- }
- vector2 = ( that2->operation!=CONST_OP );
- if( vector2 )
- vector2 = that2->value.nelem;
- else {
- val2 = that2->value.data.lng;
- }
- if( !vector1 && !vector2 ) { /* Result is a constant */
- switch( this->operation ) {
- case '~': /* Treat as == for LONGS */
- case EQ: this->value.data.log = (val1 == val2); break;
- case NE: this->value.data.log = (val1 != val2); break;
- case GT: this->value.data.log = (val1 > val2); break;
- case LT: this->value.data.log = (val1 < val2); break;
- case LTE: this->value.data.log = (val1 <= val2); break;
- case GTE: this->value.data.log = (val1 >= val2); break;
- case '+': this->value.data.lng = (val1 + val2); break;
- case '-': this->value.data.lng = (val1 - val2); break;
- case '*': this->value.data.lng = (val1 * val2); break;
- case '%':
- if( val2 ) this->value.data.lng = (val1 % val2);
- else yyerror("Divide by Zero");
- break;
- case '/':
- if( val2 ) this->value.data.lng = (val1 / val2);
- else yyerror("Divide by Zero");
- break;
- case POWER:
- this->value.data.lng = (long)pow((double)val1,(double)val2);
- break;
- case ACCUM:
- this->value.data.lng = val1;
- break;
- case DIFF:
- this->value.data.lng = 0;
- break;
- }
- this->operation=CONST_OP;
- } else if ((this->operation == ACCUM) || (this->operation == DIFF)) {
- long i, previous, curr;
- long undef;
- rows = gParse.nRows;
- nelem = this->value.nelem;
- elem = this->value.nelem * rows;
-
- Allocate_Ptrs( this );
-
- if( !gParse.status ) {
- previous = that2->value.data.lng;
- undef = (long) that2->value.undef;
-
- if (this->operation == ACCUM) {
- /* Cumulative sum of this chunk */
- for (i=0; i<elem; i++) {
- if (!that1->value.undef[i]) {
- curr = that1->value.data.lngptr[i];
- previous += curr;
- }
- this->value.data.lngptr[i] = previous;
- this->value.undef[i] = 0;
- }
- } else {
- /* Sequential difference for this chunk */
- for (i=0; i<elem; i++) {
- curr = that1->value.data.lngptr[i];
- if (that1->value.undef[i] || undef) {
- /* Either this, or previous, value was undefined */
- this->value.data.lngptr[i] = 0;
- this->value.undef[i] = 1;
- } else {
- /* Both defined, we are okay! */
- this->value.data.lngptr[i] = curr - previous;
- this->value.undef[i] = 0;
- }
- previous = curr;
- undef = that1->value.undef[i];
- }
- }
-
- /* Store final cumulant for next pass */
- that2->value.data.lng = previous;
- that2->value.undef = (char *) undef; /* XXX evil, but no harm here */
- }
-
- } else {
- rows = gParse.nRows;
- nelem = this->value.nelem;
- elem = this->value.nelem * rows;
- Allocate_Ptrs( this );
- while( rows-- && !gParse.status ) {
- while( nelem-- && !gParse.status ) {
- elem--;
- if( vector1>1 ) {
- val1 = that1->value.data.lngptr[elem];
- null1 = that1->value.undef[elem];
- } else if( vector1 ) {
- val1 = that1->value.data.lngptr[rows];
- null1 = that1->value.undef[rows];
- }
- if( vector2>1 ) {
- val2 = that2->value.data.lngptr[elem];
- null2 = that2->value.undef[elem];
- } else if( vector2 ) {
- val2 = that2->value.data.lngptr[rows];
- null2 = that2->value.undef[rows];
- }
- this->value.undef[elem] = (null1 || null2);
- switch( this->operation ) {
- case '~': /* Treat as == for LONGS */
- case EQ: this->value.data.logptr[elem] = (val1 == val2); break;
- case NE: this->value.data.logptr[elem] = (val1 != val2); break;
- case GT: this->value.data.logptr[elem] = (val1 > val2); break;
- case LT: this->value.data.logptr[elem] = (val1 < val2); break;
- case LTE: this->value.data.logptr[elem] = (val1 <= val2); break;
- case GTE: this->value.data.logptr[elem] = (val1 >= val2); break;
-
- case '+': this->value.data.lngptr[elem] = (val1 + val2); break;
- case '-': this->value.data.lngptr[elem] = (val1 - val2); break;
- case '*': this->value.data.lngptr[elem] = (val1 * val2); break;
- case '%':
- if( val2 ) this->value.data.lngptr[elem] = (val1 % val2);
- else {
- this->value.data.lngptr[elem] = 0;
- this->value.undef[elem] = 1;
- }
- break;
- case '/':
- if( val2 ) this->value.data.lngptr[elem] = (val1 / val2);
- else {
- this->value.data.lngptr[elem] = 0;
- this->value.undef[elem] = 1;
- }
- break;
- case POWER:
- this->value.data.lngptr[elem] = (long)pow((double)val1,(double)val2);
- break;
- }
- }
- nelem = this->value.nelem;
- }
- }
- if( that1->operation>0 ) {
- free( that1->value.data.ptr );
- }
- if( that2->operation>0 ) {
- free( that2->value.data.ptr );
- }
- }
- static void Do_BinOp_dbl( Node *this )
- {
- Node *that1, *that2;
- int vector1, vector2;
- double val1=0.0, val2=0.0;
- char null1=0, null2=0;
- long rows, nelem, elem;
- that1 = gParse.Nodes + this->SubNodes[0];
- that2 = gParse.Nodes + this->SubNodes[1];
- vector1 = ( that1->operation!=CONST_OP );
- if( vector1 )
- vector1 = that1->value.nelem;
- else {
- val1 = that1->value.data.dbl;
- }
- vector2 = ( that2->operation!=CONST_OP );
- if( vector2 )
- vector2 = that2->value.nelem;
- else {
- val2 = that2->value.data.dbl;
- }
- if( !vector1 && !vector2 ) { /* Result is a constant */
- switch( this->operation ) {
- case '~': this->value.data.log = ( fabs(val1-val2) < APPROX ); break;
- case EQ: this->value.data.log = (val1 == val2); break;
- case NE: this->value.data.log = (val1 != val2); break;
- case GT: this->value.data.log = (val1 > val2); break;
- case LT: this->value.data.log = (val1 < val2); break;
- case LTE: this->value.data.log = (val1 <= val2); break;
- case GTE: this->value.data.log = (val1 >= val2); break;
- case '+': this->value.data.dbl = (val1 + val2); break;
- case '-': this->value.data.dbl = (val1 - val2); break;
- case '*': this->value.data.dbl = (val1 * val2); break;
- case '%':
- if( val2 ) this->value.data.dbl = val1 - val2*((int)(val1/val2));
- else yyerror("Divide by Zero");
- break;
- case '/':
- if( val2 ) this->value.data.dbl = (val1 / val2);
- else yyerror("Divide by Zero");
- break;
- case POWER:
- this->value.data.dbl = (double)pow(val1,val2);
- break;
- case ACCUM:
- this->value.data.dbl = val1;
- break;
- case DIFF:
- this->value.data.dbl = 0;
- break;
- }
- this->operation=CONST_OP;
- } else if ((this->operation == ACCUM) || (this->operation == DIFF)) {
- long i;
- long undef;
- double previous, curr;
- rows = gParse.nRows;
- nelem = this->value.nelem;
- elem = this->value.nelem * rows;
-
- Allocate_Ptrs( this );
-
- if( !gParse.status ) {
- previous = that2->value.data.dbl;
- undef = (long) that2->value.undef;
-
- if (this->operation == ACCUM) {
- /* Cumulative sum of this chunk */
- for (i=0; i<elem; i++) {
- if (!that1->value.undef[i]) {
- curr = that1->value.data.dblptr[i];
- previous += curr;
- }
- this->value.data.dblptr[i] = previous;
- this->value.undef[i] = 0;
- }
- } else {
- /* Sequential difference for this chunk */
- for (i=0; i<elem; i++) {
- curr = that1->value.data.dblptr[i];
- if (that1->value.undef[i] || undef) {
- /* Either this, or previous, value was undefined */
- this->value.data.dblptr[i] = 0;
- this->value.undef[i] = 1;
- } else {
- /* Both defined, we are okay! */
- this->value.data.dblptr[i] = curr - previous;
- this->value.undef[i] = 0;
- }
- previous = curr;
- undef = that1->value.undef[i];
- }
- }
-
- /* Store final cumulant for next pass */
- that2->value.data.dbl = previous;
- that2->value.undef = (char *) undef; /* XXX evil, but no harm here */
- }
-
- } else {
- rows = gParse.nRows;
- nelem = this->value.nelem;
- elem = this->value.nelem * rows;
- Allocate_Ptrs( this );
- while( rows-- && !gParse.status ) {
- while( nelem-- && !gParse.status ) {
- elem--;
- if( vector1>1 ) {
- val1 = that1->value.data.dblptr[elem];
- null1 = that1->value.undef[elem];
- } else if( vector1 ) {
- val1 = that1->value.data.dblptr[rows];
- null1 = that1->value.undef[rows];
- }
- if( vector2>1 ) {
- val2 = that2->value.data.dblptr[elem];
- null2 = that2->value.undef[elem];
- } else if( vector2 ) {
- val2 = that2->value.data.dblptr[rows];
- null2 = that2->value.undef[rows];
- }
- this->value.undef[elem] = (null1 || null2);
- switch( this->operation ) {
- case '~': this->value.data.logptr[elem] =
- ( fabs(val1-val2) < APPROX ); break;
- case EQ: this->value.data.logptr[elem] = (val1 == val2); break;
- case NE: this->value.data.logptr[elem] = (val1 != val2); break;
- case GT: this->value.data.logptr[elem] = (val1 > val2); break;
- case LT: this->value.data.logptr[elem] = (val1 < val2); break;
- case LTE: this->value.data.logptr[elem] = (val1 <= val2); break;
- case GTE: this->value.data.logptr[elem] = (val1 >= val2); break;
-
- case '+': this->value.data.dblptr[elem] = (val1 + val2); break;
- case '-': this->value.data.dblptr[elem] = (val1 - val2); break;
- case '*': this->value.data.dblptr[elem] = (val1 * val2); break;
- case '%':
- if( val2 ) this->value.data.dblptr[elem] =
- val1 - val2*((int)(val1/val2));
- else {
- this->value.data.dblptr[elem] = 0.0;
- this->value.undef[elem] = 1;
- }
- break;
- case '/':
- if( val2 ) this->value.data.dblptr[elem] = (val1 / val2);
- else {
- this->value.data.dblptr[elem] = 0.0;
- this->value.undef[elem] = 1;
- }
- break;
- case POWER:
- this->value.data.dblptr[elem] = (double)pow(val1,val2);
- break;
- }
- }
- nelem = this->value.nelem;
- }
- }
- if( that1->operation>0 ) {
- free( that1->value.data.ptr );
- }
- if( that2->operation>0 ) {
- free( that2->value.data.ptr );
- }
- }
- /*
- * This Quickselect routine is based on the algorithm described in
- * "Numerical recipes in C", Second Edition,
- * Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
- * This code by Nicolas Devillard - 1998. Public domain.
- * http://ndevilla.free.fr/median/median/src/quickselect.c
- */
- #define ELEM_SWAP(a,b) { register long t=(a);(a)=(b);(b)=t; }
- /*
- * qselect_median_lng - select the median value of a long array
- *
- * This routine selects the median value of the long integer array
- * arr[]. If there are an even number of elements, the "lower median"
- * is selected.
- *
- * The array arr[] is scrambled, so users must operate on a scratch
- * array if they wish the values to be preserved.
- *
- * long arr[] - array of values
- * int n - number of elements in arr
- *
- * RETURNS: the lower median value of arr[]
- *
- */
- long qselect_median_lng(long arr[], int n)
- {
- int low, high ;
- int median;
- int middle, ll, hh;
- low = 0 ; high = n-1 ; median = (low + high) / 2;
- for (;;) {
- if (high <= low) { /* One element only */
- return arr[median];
- }
- if (high == low + 1) { /* Two elements only */
- if (arr[low] > arr[high])
- ELEM_SWAP(arr[low], arr[high]) ;
- return arr[median];
- }
- /* Find median of low, middle and high items; swap into position low */
- middle = (low + high) / 2;
- if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ;
- if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ;
- if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ;
- /* Swap low item (now in position middle) into position (low+1) */
- ELEM_SWAP(arr[middle], arr[low+1]) ;
- /* Nibble from each end towards middle, swapping items when stuck */
- ll = low + 1;
- hh = high;
- for (;;) {
- do ll++; while (arr[low] > arr[ll]) ;
- do hh--; while (arr[hh] > arr[low]) ;
- if (hh < ll)
- break;
- ELEM_SWAP(arr[ll], arr[hh]) ;
- }
- /* Swap middle item (in position low) back into correct position */
- ELEM_SWAP(arr[low], arr[hh]) ;
- /* Re-set active partition */
- if (hh <= median)
- low = ll;
- if (hh >= median)
- high = hh - 1;
- }
- }
- #undef ELEM_SWAP
- #define ELEM_SWAP(a,b) { register double t=(a);(a)=(b);(b)=t; }
- /*
- * qselect_median_dbl - select the median value of a double array
- *
- * This routine selects the median value of the double array
- * arr[]. If there are an even number of elements, the "lower median"
- * is selected.
- *
- * The array arr[] is scrambled, so users must operate on a scratch
- * array if they wish the values to be preserved.
- *
- * double arr[] - array of values
- * int n - number of elements in arr
- *
- * RETURNS: the lower median value of arr[]
- *
- */
- double qselect_median_dbl(double arr[], int n)
- {
- int low, high ;
- int median;
- int middle, ll, hh;
- low = 0 ; high = n-1 ; median = (low + high) / 2;
- for (;;) {
- if (high <= low) { /* One element only */
- return arr[median] ;
- }
- if (high == low + 1) { /* Two elements only */
- if (arr[low] > arr[high])
- ELEM_SWAP(arr[low], arr[high]) ;
- return arr[median] ;
- }
- /* Find median of low, middle and high items; swap into position low */
- middle = (low + high) / 2;
- if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ;
- if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ;
- if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ;
- /* Swap low item (now in position middle) into position (low+1) */
- ELEM_SWAP(arr[middle], arr[low+1]) ;
- /* Nibble from each end towards middle, swapping items when stuck */
- ll = low + 1;
- hh = high;
- for (;;) {
- do ll++; while (arr[low] > arr[ll]) ;
- do hh--; while (arr[hh] > arr[low]) ;
- if (hh < ll)
- break;
- ELEM_SWAP(arr[ll], arr[hh]) ;
- }
- /* Swap middle item (in position low) back into correct position */
- ELEM_SWAP(arr[low], arr[hh]) ;
- /* Re-set active partition */
- if (hh <= median)
- low = ll;
- if (hh >= median)
- high = hh - 1;
- }
- }
- #undef ELEM_SWAP
- /*
- * angsep_calc - compute angular separation between celestial coordinates
- *
- * This routine computes the angular separation between to coordinates
- * on the celestial sphere (i.e. RA and Dec). Note that all units are
- * in DEGREES, unlike the other trig functions in the calculator.
- *
- * double ra1, dec1 - RA and Dec of the first position in degrees
- * double ra2, dec2 - RA and Dec of the second position in degrees
- *
- * RETURNS: (double) angular separation in degrees
- *
- */
- double angsep_calc(double ra1, double dec1, double ra2, double dec2)
- {
- double cd;
- static double deg = 0;
- double a, sdec, sra;
-
- if (deg == 0) deg = ((double)4)*atan((double)1)/((double)180);
- /* deg = 1.0; **** UNCOMMENT IF YOU WANT RADIANS */
-
- /*
- This (commented out) algorithm uses the Low of Cosines, which becomes
- unstable for angles less than 0.1 arcsec.
-
- cd = sin(dec1*deg)*sin(dec2*deg)
- + cos(dec1*deg)*cos(dec2*deg)*cos((ra1-ra2)*deg);
- if (cd < (-1)) cd = -1;
- if (cd > (+1)) cd = +1;
- return acos(cd)/deg;
- */
- /* The algorithm is the law of Haversines. This algorithm is
- stable even when the points are close together. The normal
- Law of Cosines fails for angles around 0.1 arcsec. */
- sra = sin( (ra2 - ra1)*deg / 2 );
- sdec = sin( (dec2 - dec1)*deg / 2);
- a = sdec*sdec + cos(dec1*deg)*cos(dec2*deg)*sra*sra;
- /* Sanity checking to avoid a range error in the sqrt()'s below */
- if (a < 0) { a = 0; }
- if (a > 1) { a = 1; }
- return 2.0*atan2(sqrt(a), sqrt(1.0 - a)) / deg;
- }
- static double ran1()
- {
- static double dval = 0.0;
- double rndVal;
- if (dval == 0.0) {
- if( rand()<32768 && rand()<32768 )
- dval = 32768.0;
- else
- dval = 2147483648.0;
- }
- rndVal = (double)rand();
- while( rndVal > dval ) dval *= 2.0;
- return rndVal/dval;
- }
- /* Gaussian deviate routine from Numerical Recipes */
- static double gasdev()
- {
- static int iset = 0;
- static double gset;
- double fac, rsq, v1, v2;
- if (iset == 0) {
- do {
- v1 = 2.0*ran1()-1.0;
- v2 = 2.0*ran1()-1.0;
- rsq = v1*v1 + v2*v2;
- } while (rsq >= 1.0 || rsq == 0.0);
- fac = sqrt(-2.0*log(rsq)/rsq);
- gset = v1*fac;
- iset = 1;
- return v2*fac;
- } else {
- iset = 0;
- return gset;
- }
- }
- /* lgamma function - from Numerical Recipes */
- float gammaln(float xx)
- /* Returns the value ln Gamma[(xx)] for xx > 0. */
- {
- /*
- Internal arithmetic will be done in double precision, a nicety
- that you can omit if five-figure accuracy is good enough. */
- double x,y,tmp,ser;
- static double cof[6]={76.18009172947146,-86.50532032941677,
- 24.01409824083091,-1.231739572450155,
- 0.1208650973866179e-2,-0.5395239384953e-5};
- int j;
- y=x=xx;
- tmp=x+5.5;
- tmp -= (x+0.5)*log(tmp);
- ser=1.000000000190015;
- for (j=0;j<=5;j++) ser += cof[j]/++y;
- return (float) -tmp+log(2.5066282746310005*ser/x);
- }
- /* Poisson deviate - derived from Numerical Recipes */
- static long poidev(double xm)
- {
- static double sq, alxm, g, oldm = -1.0;
- static double pi = 0;
- double em, t, y;
- if (pi == 0) pi = ((double)4)*atan((double)1);
- if (xm < 20.0) {
- if (xm != oldm) {
- oldm = xm;
- g = exp(-xm);
- }
- em = -1;
- t = 1.0;
- do {
- em += 1;
- t *= ran1();
- } while (t > g);
- } else {
- if (xm != oldm) {
- oldm = xm;
- sq = sqrt(2.0*xm);
- alxm = log(xm);
- g = xm*alxm-gammaln( (float) (xm+1.0));
- }
- do {
- do {
- y = tan(pi*ran1());
- em = sq*y+xm;
- } while (em < 0.0);
- em = floor(em);
- t = 0.9*(1.0+y*y)*exp(em*alxm-gammaln( (float) (em+1.0) )-g);
- } while (ran1() > t);
- }
- /* Return integer version */
- return (long int) floor(em+0.5);
- }
- static void Do_Func( Node *this )
- {
- Node *theParams[MAXSUBS];
- int vector[MAXSUBS], allConst;
- lval pVals[MAXSUBS];
- char pNull[MAXSUBS];
- long ival;
- double dval;
- int i, valInit;
- long row, elem, nelem;
- i = this->nSubNodes;
- allConst = 1;
- while( i-- ) {
- theParams[i] = gParse.Nodes + this->SubNodes[i];
- vector[i] = ( theParams[i]->operation!=CONST_OP );
- if( vector[i] ) {
- allConst = 0;
- vector[i] = theParams[i]->value.nelem;
- } else {
- if( theParams[i]->type==DOUBLE ) {
- pVals[i].data.dbl = theParams[i]->value.data.dbl;
- } else if( theParams[i]->type==LONG ) {
- pVals[i].data.lng = theParams[i]->value.data.lng;
- } else if( theParams[i]->type==BOOLEAN ) {
- pVals[i].data.log = theParams[i]->value.data.log;
- } else
- strcpy(pVals[i].data.str, theParams[i]->value.data.str);
- pNull[i] = 0;
- }
- }
- if( this->nSubNodes==0 ) allConst = 0; /* These do produce scalars */
- /* Random numbers are *never* constant !! */
- if( this->operation == poirnd_fct ) allConst = 0;
- if( this->operation == gasrnd_fct ) allConst = 0;
- if( this->operation == rnd_fct ) allConst = 0;
- if( allConst ) {
- switch( this->operation ) {
- /* Non-Trig single-argument functions */
- case sum_fct:
- if( theParams[0]->type==BOOLEAN )
- this->value.data.lng = ( pVals[0].data.log ? 1 : 0 );
- else if( theParams[0]->type==LONG )
- this->value.data.lng = pVals[0].data.lng;
- else if( theParams[0]->type==DOUBLE )
- this->value.data.dbl = pVals[0].data.dbl;
- else if( theParams[0]->type==BITSTR )
- strcpy(this->value.data.str, pVals[0].data.str);
- break;
- case average_fct:
- if( theParams[0]->type==LONG )
- this->value.data.dbl = pVals[0].data.lng;
- else if( theParams[0]->type==DOUBLE )
- this->value.data.dbl = pVals[0].data.dbl;
- break;
- case stddev_fct:
- this->value.data.dbl = 0; /* Standard deviation of a constant = 0 */
- break;
- case median_fct:
- if( theParams[0]->type==BOOLEAN )
- this->value.data.lng = ( pVals[0].data.log ? 1 : 0 );
- else if( theParams[0]->type==LONG )
- this->value.data.lng = pVals[0].data.lng;
- else
- this->value.data.dbl = pVals[0].data.dbl;
- break;
- case poirnd_fct:
- if( theParams[0]->type==DOUBLE )
- this->value.data.lng = poidev(pVals[0].data.dbl);
- else
- this->value.data.lng = poidev(pVals[0].data.lng);
- break;
- case abs_fct:
- if( theParams[0]->type==DOUBLE ) {
- dval = pVals[0].data.dbl;
- this->value.data.dbl = (dval>0.0 ? dval : -dval);
- } else {
- ival = pVals[0].data.lng;
- this->value.data.lng = (ival> 0 ? ival : -ival);
- }
- break;
- /* Special Null-Handling Functions */
- case nonnull_fct:
- this->value.data.lng = 1; /* Constants are always 1-element and defined */
- break;
- case isnull_fct: /* Constants are always defined */
- this->value.data.log = 0;
- break;
- case defnull_fct:
- if( this->type==BOOLEAN )
- this->value.data.log = pVals[0].data.log;
- else if( this->type==LONG )
- this->value.data.lng = pVals[0].data.lng;
- else if( this->type==DOUBLE )
- this->value.data.dbl = pVals[0].data.dbl;
- else if( this->type==STRING )
- strcpy(this->value.data.str,pVals[0].data.str);
- break;
- /* Math functions with 1 double argument */
- case sin_fct:
- this->value.data.dbl = sin( pVals[0].data.dbl );
- break;
- case cos_fct:
- this->value.data.dbl = cos( pVals[0].data.dbl );
- break;
- case tan_fct:
- this->value.data.dbl = tan( pVals[0].data.dbl );
- break;
- case asin_fct:
- dval = pVals[0].data.dbl;
- if( dval<-1.0 || dval>1.0 )
- yyerror("Out of range argument to arcsin");
- else
- this->value.data.dbl = asin( dval );
- break;
- case acos_fct:
- dval = pVals[0].data.dbl;
- if( dval<-1.0 || dval>1.0 )
- yyerror("Out of range argument to arccos");
- else
- this->value.data.dbl = acos( dval );
- break;
- case atan_fct:
- this->value.data.dbl = atan( pVals[0].data.dbl );
- break;
- case sinh_fct:
- this->value.data.dbl = sinh( pVals[0].data.dbl );
- break;
- case cosh_fct:
- this->value.data.dbl = cosh( pVals[0].data.dbl );
- break;
- case tanh_fct:
- this->value.data.dbl = tanh( pVals[0].data.dbl );
- break;
- case exp_fct:
- this->value.data.dbl = exp( pVals[0].data.dbl );
- break;
- case log_fct:
- dval = pVals[0].data.dbl;
- if( dval<=0.0 )
- yyerror("Out of range argument to log");
- else
- this->value.data.dbl = log( dval );
- break;
- case log10_fct:
- dval = pVals[0].data.dbl;
- if( dval<=0.0 )
- yyerror("Out of range argument to log10");
- else
- this->value.data.dbl = log10( dval );
- break;
- case sqrt_fct:
- dval = pVals[0].data.dbl;
- if( dval<0.0 )
- yyerror("Out of range argument to sqrt");
- else
- this->value.data.dbl = sqrt( dval );
- break;
- case ceil_fct:
- this->value.data.dbl = ceil( pVals[0].data.dbl );
- break;
- case floor_fct:
- this->value.data.dbl = floor( pVals[0].data.dbl );
- break;
- case round_fct:
- this->value.data.dbl = floor( pVals[0].data.dbl + 0.5 );
- break;
- /* Two-argument Trig Functions */
- case atan2_fct:
- this->value.data.dbl =
- atan2( pVals[0].data.dbl, pVals[1].data.dbl );
- break;
- /* Four-argument ANGSEP function */
- case angsep_fct:
- this->value.data.dbl =
- angsep_calc(pVals[0].data.dbl, pVals[1].data.dbl,
- pVals[2].data.dbl, pVals[3].data.dbl);
- /* Min/Max functions taking 1 or 2 arguments */
- case min1_fct:
- /* No constant vectors! */
- if( this->type == DOUBLE )
- this->value.data.dbl = pVals[0].data.dbl;
- else if( this->type == LONG )
- this->value.data.lng = pVals[0].data.lng;
- else if( this->type == BITSTR )
- strcpy(this->value.data.str, pVals[0].data.str);
- break;
- case min2_fct:
- if( this->type == DOUBLE )
- this->value.data.dbl =
- minvalue( pVals[0].data.dbl, pVals[1].data.dbl );
- else if( this->type == LONG )
- this->value.data.lng =
- minvalue( pVals[0].data.lng, pVals[1].data.lng );
- break;
- case max1_fct:
- /* No constant vectors! */
- if( this->type == DOUBLE )
- this->value.data.dbl = pVals[0].data.dbl;
- else if( this->type == LONG )
- this->value.data.lng = pVals[0].data.lng;
- else if( this->type == BITSTR )
- strcpy(this->value.data.str, pVals[0].data.str);
- break;
- case max2_fct:
- if( this->type == DOUBLE )
- this->value.data.dbl =
- maxvalue( pVals[0].data.dbl, pVals[1].data.dbl );
- else if( this->type == LONG )
- this->value.data.lng =
- maxvalue( pVals[0].data.lng, pVals[1].data.lng );
- break;
- /* Boolean SAO region Functions... scalar or vector dbls */
- case near_fct:
- this->value.data.log = bnear( pVals[0].data.dbl, pVals[1].data.dbl,
- pVals[2].data.dbl );
- break;
- case circle_fct:
- this->value.data.log = circle( pVals[0].data.dbl, pVals[1].data.dbl,
- pVals[2].data.dbl, pVals[3].data.dbl,
- pVals[4].data.dbl );
- break;
- case box_fct:
- this->value.data.log = saobox( pVals[0].data.dbl, pVals[1].data.dbl,
- pVals[2].data.dbl, pVals[3].data.dbl,
- pVals[4].data.dbl, pVals[5].data.dbl,
- pVals[6].data.dbl );
- break;
- case elps_fct:
- this->value.data.log =
- ellipse( pVals[0].data.dbl, pVals[1].data.dbl,
- pVals[2].data.dbl, pVals[3].data.dbl,
- pVals[4].data.dbl, pVals[5].data.dbl,
- pVals[6].data.dbl );
- break;
- /* C Conditional expression: bool ? expr : expr */
- case ifthenelse_fct:
- switch( this->type ) {
- case BOOLEAN:
- this->value.data.log = ( pVals[2].data.log ?
- pVals[0].data.log : pVals[1].data.log );
- break;
- case LONG:
- this->value.data.lng = ( pVals[2].data.log ?
- pVals[0].data.lng : pVals[1].data.lng );
- break;
- case DOUBLE:
- this->value.data.dbl = ( pVals[2].data.log ?
- pVals[0].data.dbl : pVals[1].data.dbl );
- break;
- case STRING:
- strcpy(this->value.data.str, ( pVals[2].data.log ?
- pVals[0].data.str :
- pVals[1].data.str ) );
- break;
- }
- break;
- /* String functions */
- case strmid_fct:
- cstrmid(this->value.data.str, this->value.nelem,
- pVals[0].data.str, pVals[0].nelem,
- pVals[1].data.lng);
- break;
- case strpos_fct:
- {
- char *res = strstr(pVals[0].data.str, pVals[1].data.str);
- if (res == NULL) {
- this->value.data.lng = 0;
- } else {
- this->value.data.lng = (res - pVals[0].data.str) + 1;
- }
- break;
- }
- }
- this->operation = CONST_OP;
- } else {
- Allocate_Ptrs( this );
- row = gParse.nRows;
- elem = row * this->value.nelem;
- if( !gParse.status ) {
- switch( this->operation ) {
- /* Special functions with no arguments */
- case row_fct:
- while( row-- ) {
- this->value.data.lngptr[row] = gParse.firstRow + row;
- this->value.undef[row] = 0;
- }
- break;
- case null_fct:
- if( this->type==LONG ) {
- while( row-- ) {
- this->value.data.lngptr[row] = 0;
- this->value.undef[row] = 1;
- }
- } else if( this->type==STRING ) {
- while( row-- ) {
- this->value.data.strptr[row][0] = '\0';
- this->value.undef[row] = 1;
- }
- }
- break;
- case rnd_fct:
- while( elem-- ) {
- this->value.data.dblptr[elem] = ran1();
- this->value.undef[elem] = 0;
- }
- break;
- case gasrnd_fct:
- while( elem-- ) {
- this->value.data.dblptr[elem] = gasdev();
- this->value.undef[elem] = 0;
- }
- break;
- case poirnd_fct:
- if( theParams[0]->type==DOUBLE ) {
- if (theParams[0]->operation == CONST_OP) {
- while( elem-- ) {
- this->value.undef[elem] = (pVals[0].data.dbl < 0);
- if (! this->value.undef[elem]) {
- this->value.data.lngptr[elem] = poidev(pVals[0].data.dbl);
- }
- }
- } else {
- while( elem-- ) {
- this->value.undef[elem] = theParams[0]->value.undef[elem];
- if (theParams[0]->value.data.dblptr[elem] < 0)
- this->value.undef[elem] = 1;
- if (! this->value.undef[elem]) {
- this->value.data.lngptr[elem] =
- poidev(theParams[0]->value.data.dblptr[elem]);
- }
- } /* while */
- } /* ! CONST_OP */
- } else {
- /* LONG */
- if (theParams[0]->operation == CONST_OP) {
- while( elem-- ) {
- this->value.undef[elem] = (pVals[0].data.lng < 0);
- if (! this->value.undef[elem]) {
- this->value.data.lngptr[elem] = poidev(pVals[0].data.lng);
- }
- }
- } else {
- while( elem-- ) {
- this->value.undef[elem] = theParams[0]->value.undef[elem];
- if (theParams[0]->value.data.lngptr[elem] < 0)
- this->value.undef[elem] = 1;
- if (! this->value.undef[elem]) {
- this->value.data.lngptr[elem] =
- poidev(theParams[0]->value.data.lngptr[elem]);
- }
- } /* while */
- } /* ! CONST_OP */
- } /* END LONG */
- break;
- /* Non-Trig single-argument functions */
-
- case sum_fct:
- elem = row * theParams[0]->value.nelem;
- if( theParams[0]->type==BOOLEAN ) {
- while( row-- ) {
- this->value.data.lngptr[row] = 0;
- /* Default is UNDEF until a defined value is found */
- this->value.undef[row] = 1;
- nelem = theParams[0]->value.nelem;
- while( nelem-- ) {
- elem--;
- if ( ! theParams[0]->value.undef[elem] ) {
- this->value.data.lngptr[row] +=
- ( theParams[0]->value.data.logptr[elem] ? 1 : 0 );
- this->value.undef[row] = 0;
- }
- }
- }
- } else if( theParams[0]->type==LONG ) {
- while( row-- ) {
- this->value.data.lngptr[row] = 0;
- /* Default is UNDEF until a defined value is found */
- this->value.undef[row] = 1;
- nelem = theParams[0]->value.nelem;
- while( nelem-- ) {
- elem--;
- if ( ! theParams[0]->value.undef[elem] ) {
- this->value.data.lngptr[row] +=
- theParams[0]->value.data.lngptr[elem];
- this->value.undef[row] = 0;
- }
- }
- }
- } else if( theParams[0]->type==DOUBLE ){
- while( row-- ) {
- this->value.data.dblptr[row] = 0.0;
- /* Default is UNDEF until a defined value is found */
- this->value.undef[row] = 1;
- nelem = theParams[0]->value.nelem;
- while( nelem-- ) {
- elem--;
- if ( ! theParams[0]->value.undef[elem] ) {
- this->value.data.dblptr[row] +=
- theParams[0]->value.data.dblptr[elem];
- this->value.undef[row] = 0;
- }
- }
- }
- } else { /* BITSTR */
- nelem = theParams[0]->value.nelem;
- while( row-- ) {
- char *sptr1 = theParams[0]->value.data.strptr[row];
- this->value.data.lngptr[row] = 0;
- this->value.undef[row] = 0;
- while (*sptr1) {
- if (*sptr1 == '1') this->value.data.lngptr[row] ++;
- sptr1++;
- }
- }
- }
- break;
- case average_fct:
- elem = row * theParams[0]->value.nelem;
- if( theParams[0]->type==LONG ) {
- while( row-- ) {
- int count = 0;
- this->value.data.dblptr[row] = 0;
- nelem = theParams[0]->value.nelem;
- while( nelem-- ) {
- elem--;
- if (theParams[0]->value.undef[elem] == 0) {
- this->value.data.dblptr[row] +=
- theParams[0]->value.data.lngptr[elem];
- count ++;
- }
- }
- if (count == 0) {
- this->value.undef[row] = 1;
- } else {
- this->value.undef[row] = 0;
- this->value.data.dblptr[row] /= count;
- }
- }
- } else if( theParams[0]->type==DOUBLE ){
- while( row-- ) {
- int count = 0;
- this->value.data.dblptr[row] = 0;
- nelem = theParams[0]->value.nelem;
- while( nelem-- ) {
- elem--;
- if (theParams[0]->value.undef[elem] == 0) {
- this->value.data.dblptr[row] +=
- theParams[0]->value.data.dblptr[elem];
- count ++;
- }
- }
- if (count == 0) {
- this->value.undef[row] = 1;
- } else {
- this->value.undef[row] = 0;
- this->value.data.dblptr[row] /= count;
- }
- }
- }
- break;
- case stddev_fct:
- elem = row * theParams[0]->value.nelem;
- if( theParams[0]->type==LONG ) {
- /* Compute the mean value */
- while( row-- ) {
- int count = 0;
- double sum = 0, sum2 = 0;
- nelem = theParams[0]->value.nelem;
- while( nelem-- ) {
- elem--;
- if (theParams[0]->value.undef[elem] == 0) {
- sum += theParams[0]->value.data.lngptr[elem];
- count ++;
- }
- }
- if (count > 1) {
- sum /= count;
- /* Compute the sum of squared deviations */
- nelem = theParams[0]->value.nelem;
- elem += nelem; /* Reset elem for second pass */
- while( nelem-- ) {
- elem--;
- if (theParams[0]->value.undef[elem] == 0) {
- double dx = (theParams[0]->value.data.lngptr[elem] - sum);
- sum2 += (dx*dx);
- }
- }
- sum2 /= (double)count-1;
- this->value.undef[row] = 0;
- this->value.data.dblptr[row] = sqrt(sum2);
- } else {
- this->value.undef[row] = 0; /* STDDEV => 0 */
- this->value.data.dblptr[row] = 0;
- }
- }
- } else if( theParams[0]->type==DOUBLE ){
- /* Compute the mean value */
- while( row-- ) {
- int count = 0;
- double sum = 0, sum2 = 0;
- nelem = theParams[0]->value.nelem;
- while( nelem-- ) {
- elem--;
- if (theParams[0]->value.undef[elem] == 0) {
- sum += theParams[0]->value.data.dblptr[elem];
- count ++;
- }
- }
- if (count > 1) {
- sum /= count;
- /* Compute the sum of squared deviations */
- nelem = theParams[0]->value.nelem;
- elem += nelem; /* Reset elem for second pass */
- while( nelem-- ) {
- elem--;
- if (theParams[0]->value.undef[elem] == 0) {
- double dx = (theParams[0]->value.data.dblptr[elem] - sum);
- sum2 += (dx*dx);
- }
- }
- sum2 /= (double)count-1;
- this->value.undef[row] = 0;
- this->value.data.dblptr[row] = sqrt(sum2);
- } else {
- this->value.undef[row] = 0; /* STDDEV => 0 */
- this->value.data.dblptr[row] = 0;
- }
- }
- }
- break;
- case median_fct:
- elem = row * theParams[0]->value.nelem;
- nelem = theParams[0]->value.nelem;
- if( theParams[0]->type==LONG ) {
- long *dptr = theParams[0]->value.data.lngptr;
- char *uptr = theParams[0]->value.undef;
- long *mptr = (long *) malloc(sizeof(long)*nelem);
- int irow;
- /* Allocate temporary storage for this row, since the
- quickselect function will scramble the contents */
- if (mptr == 0) {
- yyerror("Could not allocate temporary memory in median function");
- free( this->value.data.ptr );
- break;
- }
- for (irow=0; irow<row; irow++) {
- long *p = mptr;
- int nelem1 = nelem;
- int count = 0;
- while ( nelem1-- ) {
- if (*uptr == 0) {
- *p++ = *dptr; /* Only advance the dest pointer if we copied */
- }
- dptr ++; /* Advance the source pointer ... */
- uptr ++; /* ... and source "undef" pointer */
- }
-
- nelem1 = (p - mptr); /* Number of accepted data points */
- if (nelem1 > 0) {
- this->value.undef[irow] = 0;
- this->value.data.lngptr[irow] = qselect_median_lng(mptr, nelem1);
- } else {
- this->value.undef[irow] = 1;
- this->value.data.lngptr[irow] = 0;
- }
-
- }
- free(mptr);
- } else {
- double *dptr = theParams[0]->value.data.dblptr;
- char *uptr = theParams[0]->value.undef;
- double *mptr = (double *) malloc(sizeof(double)*nelem);
- int irow;
- /* Allocate temporary storage for this row, since the
- quickselect function will scramble the contents */
- if (mptr == 0) {
- yyerror("Could not allocate temporary memory in median function");
- free( this->value.data.ptr );
- break;
- }
- for (irow=0; irow<row; irow++) {
- double *p = mptr;
- int nelem1 = nelem;
- while ( nelem1-- ) {
- if (*uptr == 0) {
- *p++ = *dptr; /* Only advance the dest pointer if we copied */
- }
- dptr ++; /* Advance the source pointer ... */
- uptr ++; /* ... and source "undef" pointer */
- }
- nelem1 = (p - mptr); /* Number of accepted data points */
- if (nelem1 > 0) {
- this->value.undef[irow] = 0;
- this->value.data.dblptr[irow] = qselect_median_dbl(mptr, nelem1);
- } else {
- this->value.undef[irow] = 1;
- this->value.data.dblptr[irow] = 0;
- }
- }
- free(mptr);
- }
- break;
- case abs_fct:
- if( theParams[0]->type==DOUBLE )
- while( elem-- ) {
- dval = theParams[0]->value.data.dblptr[elem];
- this->value.data.dblptr[elem] = (dval>0.0 ? dval : -dval);
- this->value.undef[elem] = theParams[0]->value.undef[elem];
- }
- else
- while( elem-- ) {
- ival = theParams[0]->value.data.lngptr[elem];
- this->value.data.lngptr[elem] = (ival> 0 ? ival : -ival);
- this->value.undef[elem] = theParams[0]->value.undef[elem];
- }
- break;
- /* Special Null-Handling Functions */
- case nonnull_fct:
- nelem = theParams[0]->value.nelem;
- if ( theParams[0]->type==STRING ) nelem = 1;
- elem = row * nelem;
- while( row-- ) {
- int nelem1 = nelem;
- this->value.undef[row] = 0; /* Initialize to 0 (defined) */
- this->value.data.lngptr[row] = 0;
- while( nelem1-- ) {
- elem --;
- if ( theParams[0]->value.undef[elem] == 0 ) this->value.data.lngptr[row] ++;
- }
- }
- break;
- case isnull_fct:
- if( theParams[0]->type==STRING ) elem = row;
- while( elem-- ) {
- this->value.data.logptr[elem] = theParams[0]->value.undef[elem];
- this->value.undef[elem] = 0;
- }
- break;
- case defnull_fct:
- switch( this->type ) {
- case BOOLEAN:
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- i=2; while( i-- )
- if( vector[i]>1 ) {
- pNull[i] = theParams[i]->value.undef[elem];
- pVals[i].data.log =
- theParams[i]->value.data.logptr[elem];
- } else if( vector[i] ) {
- pNull[i] = theParams[i]->value.undef[row];
- pVals[i].data.log =
- theParams[i]->value.data.logptr[row];
- }
- if( pNull[0] ) {
- this->value.undef[elem] = pNull[1];
- this->value.data.logptr[elem] = pVals[1].data.log;
- } else {
- this->value.undef[elem] = 0;
- this->value.data.logptr[elem] = pVals[0].data.log;
- }
- }
- }
- break;
- case LONG:
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- i=2; while( i-- )
- if( vector[i]>1 ) {
- pNull[i] = theParams[i]->value.undef[elem];
- pVals[i].data.lng =
- theParams[i]->value.data.lngptr[elem];
- } else if( vector[i] ) {
- pNull[i] = theParams[i]->value.undef[row];
- pVals[i].data.lng =
- theParams[i]->value.data.lngptr[row];
- }
- if( pNull[0] ) {
- this->value.undef[elem] = pNull[1];
- this->value.data.lngptr[elem] = pVals[1].data.lng;
- } else {
- this->value.undef[elem] = 0;
- this->value.data.lngptr[elem] = pVals[0].data.lng;
- }
- }
- }
- break;
- case DOUBLE:
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- i=2; while( i-- )
- if( vector[i]>1 ) {
- pNull[i] = theParams[i]->value.undef[elem];
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[elem];
- } else if( vector[i] ) {
- pNull[i] = theParams[i]->value.undef[row];
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[row];
- }
- if( pNull[0] ) {
- this->value.undef[elem] = pNull[1];
- this->value.data.dblptr[elem] = pVals[1].data.dbl;
- } else {
- this->value.undef[elem] = 0;
- this->value.data.dblptr[elem] = pVals[0].data.dbl;
- }
- }
- }
- break;
- case STRING:
- while( row-- ) {
- i=2; while( i-- )
- if( vector[i] ) {
- pNull[i] = theParams[i]->value.undef[row];
- strcpy(pVals[i].data.str,
- theParams[i]->value.data.strptr[row]);
- }
- if( pNull[0] ) {
- this->value.undef[row] = pNull[1];
- strcpy(this->value.data.strptr[row],pVals[1].data.str);
- } else {
- this->value.undef[elem] = 0;
- strcpy(this->value.data.strptr[row],pVals[0].data.str);
- }
- }
- }
- break;
- /* Math functions with 1 double argument */
- case sin_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- this->value.data.dblptr[elem] =
- sin( theParams[0]->value.data.dblptr[elem] );
- }
- break;
- case cos_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- this->value.data.dblptr[elem] =
- cos( theParams[0]->value.data.dblptr[elem] );
- }
- break;
- case tan_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- this->value.data.dblptr[elem] =
- tan( theParams[0]->value.data.dblptr[elem] );
- }
- break;
- case asin_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- dval = theParams[0]->value.data.dblptr[elem];
- if( dval<-1.0 || dval>1.0 ) {
- this->value.data.dblptr[elem] = 0.0;
- this->value.undef[elem] = 1;
- } else
- this->value.data.dblptr[elem] = asin( dval );
- }
- break;
- case acos_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- dval = theParams[0]->value.data.dblptr[elem];
- if( dval<-1.0 || dval>1.0 ) {
- this->value.data.dblptr[elem] = 0.0;
- this->value.undef[elem] = 1;
- } else
- this->value.data.dblptr[elem] = acos( dval );
- }
- break;
- case atan_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- dval = theParams[0]->value.data.dblptr[elem];
- this->value.data.dblptr[elem] = atan( dval );
- }
- break;
- case sinh_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- this->value.data.dblptr[elem] =
- sinh( theParams[0]->value.data.dblptr[elem] );
- }
- break;
- case cosh_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- this->value.data.dblptr[elem] =
- cosh( theParams[0]->value.data.dblptr[elem] );
- }
- break;
- case tanh_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- this->value.data.dblptr[elem] =
- tanh( theParams[0]->value.data.dblptr[elem] );
- }
- break;
- case exp_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- dval = theParams[0]->value.data.dblptr[elem];
- this->value.data.dblptr[elem] = exp( dval );
- }
- break;
- case log_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- dval = theParams[0]->value.data.dblptr[elem];
- if( dval<=0.0 ) {
- this->value.data.dblptr[elem] = 0.0;
- this->value.undef[elem] = 1;
- } else
- this->value.data.dblptr[elem] = log( dval );
- }
- break;
- case log10_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- dval = theParams[0]->value.data.dblptr[elem];
- if( dval<=0.0 ) {
- this->value.data.dblptr[elem] = 0.0;
- this->value.undef[elem] = 1;
- } else
- this->value.data.dblptr[elem] = log10( dval );
- }
- break;
- case sqrt_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- dval = theParams[0]->value.data.dblptr[elem];
- if( dval<0.0 ) {
- this->value.data.dblptr[elem] = 0.0;
- this->value.undef[elem] = 1;
- } else
- this->value.data.dblptr[elem] = sqrt( dval );
- }
- break;
- case ceil_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- this->value.data.dblptr[elem] =
- ceil( theParams[0]->value.data.dblptr[elem] );
- }
- break;
- case floor_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- this->value.data.dblptr[elem] =
- floor( theParams[0]->value.data.dblptr[elem] );
- }
- break;
- case round_fct:
- while( elem-- )
- if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
- this->value.data.dblptr[elem] =
- floor( theParams[0]->value.data.dblptr[elem] + 0.5);
- }
- break;
- /* Two-argument Trig Functions */
-
- case atan2_fct:
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- i=2; while( i-- )
- if( vector[i]>1 ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[elem];
- pNull[i] = theParams[i]->value.undef[elem];
- } else if( vector[i] ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[row];
- pNull[i] = theParams[i]->value.undef[row];
- }
- if( !(this->value.undef[elem] = (pNull[0] || pNull[1]) ) )
- this->value.data.dblptr[elem] =
- atan2( pVals[0].data.dbl, pVals[1].data.dbl );
- }
- }
- break;
- /* Four-argument ANGSEP Function */
-
- case angsep_fct:
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- i=4; while( i-- )
- if( vector[i]>1 ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[elem];
- pNull[i] = theParams[i]->value.undef[elem];
- } else if( vector[i] ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[row];
- pNull[i] = theParams[i]->value.undef[row];
- }
- if( !(this->value.undef[elem] = (pNull[0] || pNull[1] ||
- pNull[2] || pNull[3]) ) )
- this->value.data.dblptr[elem] =
- angsep_calc(pVals[0].data.dbl, pVals[1].data.dbl,
- pVals[2].data.dbl, pVals[3].data.dbl);
- }
- }
- break;
- /* Min/Max functions taking 1 or 2 arguments */
- case min1_fct:
- elem = row * theParams[0]->value.nelem;
- if( this->type==LONG ) {
- long minVal=0;
- while( row-- ) {
- valInit = 1;
- this->value.undef[row] = 1;
- nelem = theParams[0]->value.nelem;
- while( nelem-- ) {
- elem--;
- if ( !theParams[0]->value.undef[elem] ) {
- if ( valInit ) {
- valInit = 0;
- minVal = theParams[0]->value.data.lngptr[elem];
- } else {
- minVal = minvalue( minVal,
- theParams[0]->value.data.lngptr[elem] );
- }
- this->value.undef[row] = 0;
- }
- }
- this->value.data.lngptr[row] = minVal;
- }
- } else if( this->type==DOUBLE ) {
- double minVal=0.0;
- while( row-- ) {
- valInit = 1;
- this->value.undef[row] = 1;
- nelem = theParams[0]->value.nelem;
- while( nelem-- ) {
- elem--;
- if ( !theParams[0]->value.undef[elem] ) {
- if ( valInit ) {
- valInit = 0;
- minVal = theParams[0]->value.data.dblptr[elem];
- } else {
- minVal = minvalue( minVal,
- theParams[0]->value.data.dblptr[elem] );
- }
- this->value.undef[row] = 0;
- }
- }
- this->value.data.dblptr[row] = minVal;
- }
- } else if( this->type==BITSTR ) {
- char minVal;
- while( row-- ) {
- char *sptr1 = theParams[0]->value.data.strptr[row];
- minVal = '1';
- while (*sptr1) {
- if (*sptr1 == '0') minVal = '0';
- sptr1++;
- }
- this->value.data.strptr[row][0] = minVal;
- this->value.data.strptr[row][1] = 0; /* Null terminate */
- }
- }
- break;
- case min2_fct:
- if( this->type==LONG ) {
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- i=2; while( i-- )
- if( vector[i]>1 ) {
- pVals[i].data.lng =
- theParams[i]->value.data.lngptr[elem];
- pNull[i] = theParams[i]->value.undef[elem];
- } else if( vector[i] ) {
- pVals[i].data.lng =
- theParams[i]->value.data.lngptr[row];
- pNull[i] = theParams[i]->value.undef[row];
- }
- if( pNull[0] && pNull[1] ) {
- this->value.undef[elem] = 1;
- this->value.data.lngptr[elem] = 0;
- } else if (pNull[0]) {
- this->value.undef[elem] = 0;
- this->value.data.lngptr[elem] = pVals[1].data.lng;
- } else if (pNull[1]) {
- this->value.undef[elem] = 0;
- this->value.data.lngptr[elem] = pVals[0].data.lng;
- } else {
- this->value.undef[elem] = 0;
- this->value.data.lngptr[elem] =
- minvalue( pVals[0].data.lng, pVals[1].data.lng );
- }
- }
- }
- } else if( this->type==DOUBLE ) {
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- i=2; while( i-- )
- if( vector[i]>1 ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[elem];
- pNull[i] = theParams[i]->value.undef[elem];
- } else if( vector[i] ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[row];
- pNull[i] = theParams[i]->value.undef[row];
- }
- if( pNull[0] && pNull[1] ) {
- this->value.undef[elem] = 1;
- this->value.data.dblptr[elem] = 0;
- } else if (pNull[0]) {
- this->value.undef[elem] = 0;
- this->value.data.dblptr[elem] = pVals[1].data.dbl;
- } else if (pNull[1]) {
- this->value.undef[elem] = 0;
- this->value.data.dblptr[elem] = pVals[0].data.dbl;
- } else {
- this->value.undef[elem] = 0;
- this->value.data.dblptr[elem] =
- minvalue( pVals[0].data.dbl, pVals[1].data.dbl );
- }
- }
- }
- }
- break;
- case max1_fct:
- elem = row * theParams[0]->value.nelem;
- if( this->type==LONG ) {
- long maxVal=0;
- while( row-- ) {
- valInit = 1;
- this->value.undef[row] = 1;
- nelem = theParams[0]->value.nelem;
- while( nelem-- ) {
- elem--;
- if ( !theParams[0]->value.undef[elem] ) {
- if ( valInit ) {
- valInit = 0;
- maxVal = theParams[0]->value.data.lngptr[elem];
- } else {
- maxVal = maxvalue( maxVal,
- theParams[0]->value.data.lngptr[elem] );
- }
- this->value.undef[row] = 0;
- }
- }
- this->value.data.lngptr[row] = maxVal;
- }
- } else if( this->type==DOUBLE ) {
- double maxVal=0.0;
- while( row-- ) {
- valInit = 1;
- this->value.undef[row] = 1;
- nelem = theParams[0]->value.nelem;
- while( nelem-- ) {
- elem--;
- if ( !theParams[0]->value.undef[elem] ) {
- if ( valInit ) {
- valInit = 0;
- maxVal = theParams[0]->value.data.dblptr[elem];
- } else {
- maxVal = maxvalue( maxVal,
- theParams[0]->value.data.dblptr[elem] );
- }
- this->value.undef[row] = 0;
- }
- }
- this->value.data.dblptr[row] = maxVal;
- }
- } else if( this->type==BITSTR ) {
- char maxVal;
- while( row-- ) {
- char *sptr1 = theParams[0]->value.data.strptr[row];
- maxVal = '0';
- while (*sptr1) {
- if (*sptr1 == '1') maxVal = '1';
- sptr1++;
- }
- this->value.data.strptr[row][0] = maxVal;
- this->value.data.strptr[row][1] = 0; /* Null terminate */
- }
- }
- break;
- case max2_fct:
- if( this->type==LONG ) {
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- i=2; while( i-- )
- if( vector[i]>1 ) {
- pVals[i].data.lng =
- theParams[i]->value.data.lngptr[elem];
- pNull[i] = theParams[i]->value.undef[elem];
- } else if( vector[i] ) {
- pVals[i].data.lng =
- theParams[i]->value.data.lngptr[row];
- pNull[i] = theParams[i]->value.undef[row];
- }
- if( pNull[0] && pNull[1] ) {
- this->value.undef[elem] = 1;
- this->value.data.lngptr[elem] = 0;
- } else if (pNull[0]) {
- this->value.undef[elem] = 0;
- this->value.data.lngptr[elem] = pVals[1].data.lng;
- } else if (pNull[1]) {
- this->value.undef[elem] = 0;
- this->value.data.lngptr[elem] = pVals[0].data.lng;
- } else {
- this->value.undef[elem] = 0;
- this->value.data.lngptr[elem] =
- maxvalue( pVals[0].data.lng, pVals[1].data.lng );
- }
- }
- }
- } else if( this->type==DOUBLE ) {
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- i=2; while( i-- )
- if( vector[i]>1 ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[elem];
- pNull[i] = theParams[i]->value.undef[elem];
- } else if( vector[i] ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[row];
- pNull[i] = theParams[i]->value.undef[row];
- }
- if( pNull[0] && pNull[1] ) {
- this->value.undef[elem] = 1;
- this->value.data.dblptr[elem] = 0;
- } else if (pNull[0]) {
- this->value.undef[elem] = 0;
- this->value.data.dblptr[elem] = pVals[1].data.dbl;
- } else if (pNull[1]) {
- this->value.undef[elem] = 0;
- this->value.data.dblptr[elem] = pVals[0].data.dbl;
- } else {
- this->value.undef[elem] = 0;
- this->value.data.dblptr[elem] =
- maxvalue( pVals[0].data.dbl, pVals[1].data.dbl );
- }
- }
- }
- }
- break;
- /* Boolean SAO region Functions... scalar or vector dbls */
- case near_fct:
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- i=3; while( i-- )
- if( vector[i]>1 ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[elem];
- pNull[i] = theParams[i]->value.undef[elem];
- } else if( vector[i] ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[row];
- pNull[i] = theParams[i]->value.undef[row];
- }
- if( !(this->value.undef[elem] = (pNull[0] || pNull[1] ||
- pNull[2]) ) )
- this->value.data.logptr[elem] =
- bnear( pVals[0].data.dbl, pVals[1].data.dbl,
- pVals[2].data.dbl );
- }
- }
- break;
- case circle_fct:
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- i=5; while( i-- )
- if( vector[i]>1 ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[elem];
- pNull[i] = theParams[i]->value.undef[elem];
- } else if( vector[i] ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[row];
- pNull[i] = theParams[i]->value.undef[row];
- }
- if( !(this->value.undef[elem] = (pNull[0] || pNull[1] ||
- pNull[2] || pNull[3] ||
- pNull[4]) ) )
- this->value.data.logptr[elem] =
- circle( pVals[0].data.dbl, pVals[1].data.dbl,
- pVals[2].data.dbl, pVals[3].data.dbl,
- pVals[4].data.dbl );
- }
- }
- break;
- case box_fct:
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- i=7; while( i-- )
- if( vector[i]>1 ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[elem];
- pNull[i] = theParams[i]->value.undef[elem];
- } else if( vector[i] ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[row];
- pNull[i] = theParams[i]->value.undef[row];
- }
- if( !(this->value.undef[elem] = (pNull[0] || pNull[1] ||
- pNull[2] || pNull[3] ||
- pNull[4] || pNull[5] ||
- pNull[6] ) ) )
- this->value.data.logptr[elem] =
- saobox( pVals[0].data.dbl, pVals[1].data.dbl,
- pVals[2].data.dbl, pVals[3].data.dbl,
- pVals[4].data.dbl, pVals[5].data.dbl,
- pVals[6].data.dbl );
- }
- }
- break;
- case elps_fct:
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- i=7; while( i-- )
- if( vector[i]>1 ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[elem];
- pNull[i] = theParams[i]->value.undef[elem];
- } else if( vector[i] ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[row];
- pNull[i] = theParams[i]->value.undef[row];
- }
- if( !(this->value.undef[elem] = (pNull[0] || pNull[1] ||
- pNull[2] || pNull[3] ||
- pNull[4] || pNull[5] ||
- pNull[6] ) ) )
- this->value.data.logptr[elem] =
- ellipse( pVals[0].data.dbl, pVals[1].data.dbl,
- pVals[2].data.dbl, pVals[3].data.dbl,
- pVals[4].data.dbl, pVals[5].data.dbl,
- pVals[6].data.dbl );
- }
- }
- break;
- /* C Conditional expression: bool ? expr : expr */
- case ifthenelse_fct:
- switch( this->type ) {
- case BOOLEAN:
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- if( vector[2]>1 ) {
- pVals[2].data.log =
- theParams[2]->value.data.logptr[elem];
- pNull[2] = theParams[2]->value.undef[elem];
- } else if( vector[2] ) {
- pVals[2].data.log =
- theParams[2]->value.data.logptr[row];
- pNull[2] = theParams[2]->value.undef[row];
- }
- i=2; while( i-- )
- if( vector[i]>1 ) {
- pVals[i].data.log =
- theParams[i]->value.data.logptr[elem];
- pNull[i] = theParams[i]->value.undef[elem];
- } else if( vector[i] ) {
- pVals[i].data.log =
- theParams[i]->value.data.logptr[row];
- pNull[i] = theParams[i]->value.undef[row];
- }
- if( !(this->value.undef[elem] = pNull[2]) ) {
- if( pVals[2].data.log ) {
- this->value.data.logptr[elem] = pVals[0].data.log;
- this->value.undef[elem] = pNull[0];
- } else {
- this->value.data.logptr[elem] = pVals[1].data.log;
- this->value.undef[elem] = pNull[1];
- }
- }
- }
- }
- break;
- case LONG:
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- if( vector[2]>1 ) {
- pVals[2].data.log =
- theParams[2]->value.data.logptr[elem];
- pNull[2] = theParams[2]->value.undef[elem];
- } else if( vector[2] ) {
- pVals[2].data.log =
- theParams[2]->value.data.logptr[row];
- pNull[2] = theParams[2]->value.undef[row];
- }
- i=2; while( i-- )
- if( vector[i]>1 ) {
- pVals[i].data.lng =
- theParams[i]->value.data.lngptr[elem];
- pNull[i] = theParams[i]->value.undef[elem];
- } else if( vector[i] ) {
- pVals[i].data.lng =
- theParams[i]->value.data.lngptr[row];
- pNull[i] = theParams[i]->value.undef[row];
- }
- if( !(this->value.undef[elem] = pNull[2]) ) {
- if( pVals[2].data.log ) {
- this->value.data.lngptr[elem] = pVals[0].data.lng;
- this->value.undef[elem] = pNull[0];
- } else {
- this->value.data.lngptr[elem] = pVals[1].data.lng;
- this->value.undef[elem] = pNull[1];
- }
- }
- }
- }
- break;
- case DOUBLE:
- while( row-- ) {
- nelem = this->value.nelem;
- while( nelem-- ) {
- elem--;
- if( vector[2]>1 ) {
- pVals[2].data.log =
- theParams[2]->value.data.logptr[elem];
- pNull[2] = theParams[2]->value.undef[elem];
- } else if( vector[2] ) {
- pVals[2].data.log =
- theParams[2]->value.data.logptr[row];
- pNull[2] = theParams[2]->value.undef[row];
- }
- i=2; while( i-- )
- if( vector[i]>1 ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[elem];
- pNull[i] = theParams[i]->value.undef[elem];
- } else if( vector[i] ) {
- pVals[i].data.dbl =
- theParams[i]->value.data.dblptr[row];
- pNull[i] = theParams[i]->value.undef[row];
- }
- if( !(this->value.undef[elem] = pNull[2]) ) {
- if( pVals[2].data.log ) {
- this->value.data.dblptr[elem] = pVals[0].data.dbl;
- this->value.undef[elem] = pNull[0];
- } else {
- this->value.data.dblptr[elem] = pVals[1].data.dbl;
- this->value.undef[elem] = pNull[1];
- }
- }
- }
- }
- break;
- case STRING:
- while( row-- ) {
- if( vector[2] ) {
- pVals[2].data.log = theParams[2]->value.data.logptr[row];
- pNull[2] = theParams[2]->value.undef[row];
- }
- i=2; while( i-- )
- if( vector[i] ) {
- strcpy( pVals[i].data.str,
- theParams[i]->value.data.strptr[row] );
- pNull[i] = theParams[i]->value.undef[row];
- }
- if( !(this->value.undef[row] = pNull[2]) ) {
- if( pVals[2].data.log ) {
- strcpy( this->value.data.strptr[row],
- pVals[0].data.str );
- this->value.undef[row] = pNull[0];
- } else {
- strcpy( this->value.data.strptr[row],
- pVals[1].data.str );
- this->value.undef[row] = pNull[1];
- }
- } else {
- this->value.data.strptr[row][0] = '\0';
- }
- }
- break;
- }
- break;
- /* String functions */
- case strmid_fct:
- {
- int strconst = theParams[0]->operation == CONST_OP;
- int posconst = theParams[1]->operation == CONST_OP;
- int lenconst = theParams[2]->operation == CONST_OP;
- int dest_len = this->value.nelem;
- int src_len = theParams[0]->value.nelem;
- while (row--) {
- int pos;
- int len;
- char *str;
- int undef = 0;
- if (posconst) {
- pos = theParams[1]->value.data.lng;
- } else {
- pos = theParams[1]->value.data.lngptr[row];
- if (theParams[1]->value.undef[row]) undef = 1;
- }
- if (strconst) {
- str = theParams[0]->value.data.str;
- if (src_len == 0) src_len = strlen(str);
- } else {
- str = theParams[0]->value.data.strptr[row];
- if (theParams[0]->value.undef[row]) undef = 1;
- }
- if (lenconst) {
- len = dest_len;
- } else {
- len = theParams[2]->value.data.lngptr[row];
- if (theParams[2]->value.undef[row]) undef = 1;
- }
- this->value.data.strptr[row][0] = '\0';
- if (pos == 0) undef = 1;
- if (! undef ) {
- if (cstrmid(this->value.data.strptr[row], len,
- str, src_len, pos) < 0) break;
- }
- this->value.undef[row] = undef;
- }
- }
- break;
- /* String functions */
- case strpos_fct:
- {
- int const1 = theParams[0]->operation == CONST_OP;
- int const2 = theParams[1]->operation == CONST_OP;
- while (row--) {
- char *str1, *str2;
- int undef = 0;
- if (const1) {
- str1 = theParams[0]->value.data.str;
- } else {
- str1 = theParams[0]->value.data.strptr[row];
- if (theParams[0]->value.undef[row]) undef = 1;
- }
- if (const2) {
- str2 = theParams[1]->value.data.str;
- } else {
- str2 = theParams[1]->value.data.strptr[row];
- if (theParams[1]->value.undef[row]) undef = 1;
- }
- this->value.data.lngptr[row] = 0;
- if (! undef ) {
- char *res = strstr(str1, str2);
- if (res == NULL) {
- undef = 1;
- this->value.data.lngptr[row] = 0;
- } else {
- this->value.data.lngptr[row] = (res - str1) + 1;
- }
- }
- this->value.undef[row] = undef;
- }
- }
- break;
-
- } /* End switch(this->operation) */
- } /* End if (!gParse.status) */
- } /* End non-constant operations */
- i = this->nSubNodes;
- while( i-- ) {
- if( theParams[i]->operation>0 ) {
- /* Currently only numeric params allowed */
- free( theParams[i]->value.data.ptr );
- }
- }
- }
- static void Do_Deref( Node *this )
- {
- Node *theVar, *theDims[MAXDIMS];
- int isConst[MAXDIMS], allConst;
- long dimVals[MAXDIMS];
- int i, nDims;
- long row, elem, dsize;
- theVar = gParse.Nodes + this->SubNodes[0];
- i = nDims = this->nSubNodes-1;
- allConst = 1;
- while( i-- ) {
- theDims[i] = gParse.Nodes + this->SubNodes[i+1];
- isConst[i] = ( theDims[i]->operation==CONST_OP );
- if( isConst[i] )
- dimVals[i] = theDims[i]->value.data.lng;
- else
- allConst = 0;
- }
- if( this->type==DOUBLE ) {
- dsize = sizeof( double );
- } else if( this->type==LONG ) {
- dsize = sizeof( long );
- } else if( this->type==BOOLEAN ) {
- dsize = sizeof( char );
- } else
- dsize = 0;
- Allocate_Ptrs( this );
- if( !gParse.status ) {
- if( allConst && theVar->value.naxis==nDims ) {
- /* Dereference completely using constant indices */
- elem = 0;
- i = nDims;
- while( i-- ) {
- if( dimVals[i]<1 || dimVals[i]>theVar->value.naxes[i] ) break;
- elem = theVar->value.naxes[i]*elem + dimVals[i]-1;
- }
- if( i<0 ) {
- for( row=0; row<gParse.nRows; row++ ) {
- if( this->type==STRING )
- this->value.undef[row] = theVar->value.undef[row];
- else if( this->type==BITSTR )
- this->value.undef; /* Dummy - BITSTRs do not have undefs */
- else
- this->value.undef[row] = theVar->value.undef[elem];
- if( this->type==DOUBLE )
- this->value.data.dblptr[row] =
- theVar->value.data.dblptr[elem];
- else if( this->type==LONG )
- this->value.data.lngptr[row] =
- theVar->value.data.lngptr[elem];
- else if( this->type==BOOLEAN )
- this->value.data.logptr[row] =
- theVar->value.data.logptr[elem];
- else {
- /* XXX Note, the below expression uses knowledge of
- the layout of the string format, namely (nelem+1)
- characters per string, followed by (nelem+1)
- "undef" values. */
- this->value.data.strptr[row][0] =
- theVar->value.data.strptr[0][elem+row];
- this->value.data.strptr[row][1] = 0; /* Null terminate */
- }
- elem += theVar->value.nelem;
- }
- } else {
- yyerror("Index out of range");
- free( this->value.data.ptr );
- }
-
- } else if( allConst && nDims==1 ) {
-
- /* Reduce dimensions by 1, using a constant index */
-
- if( dimVals[0] < 1 ||
- dimVals[0] > theVar->value.naxes[ theVar->value.naxis-1 ] ) {
- yyerror("Index out of range");
- free( this->value.data.ptr );
- } else if ( this->type == BITSTR || this->type == STRING ) {
- elem = this->value.nelem * (dimVals[0]-1);
- for( row=0; row<gParse.nRows; row++ ) {
- if (this->value.undef)
- this->value.undef[row] = theVar->value.undef[row];
- memcpy( (char*)this->value.data.strptr[0]
- + row*sizeof(char)*(this->value.nelem+1),
- (char*)theVar->value.data.strptr[0] + elem*sizeof(char),
- this->value.nelem * sizeof(char) );
- /* Null terminate */
- this->value.data.strptr[row][this->value.nelem] = 0;
- elem += theVar->value.nelem+1;
- }
- } else {
- elem = this->value.nelem * (dimVals[0]-1);
- for( row=0; row<gParse.nRows; row++ ) {
- memcpy( this->value.undef + row*this->value.nelem,
- theVar->value.undef + elem,
- this->value.nelem * sizeof(char) );
- memcpy( (char*)this->value.data.ptr
- + row*dsize*this->value.nelem,
- (char*)theVar->value.data.ptr + elem*dsize,
- this->value.nelem * dsize );
- elem += theVar->value.nelem;
- }
- }
-
- } else if( theVar->value.naxis==nDims ) {
- /* Dereference completely using an expression for the indices */
- for( row=0; row<gParse.nRows; row++ ) {
- for( i=0; i<nDims; i++ ) {
- if( !isConst[i] ) {
- if( theDims[i]->value.undef[row] ) {
- yyerror("Null encountered as vector index");
- free( this->value.data.ptr );
- break;
- } else
- dimVals[i] = theDims[i]->value.data.lngptr[row];
- }
- }
- if( gParse.status ) break;
- elem = 0;
- i = nDims;
- while( i-- ) {
- if( dimVals[i]<1 || dimVals[i]>theVar->value.naxes[i] ) break;
- elem = theVar->value.naxes[i]*elem + dimVals[i]-1;
- }
- if( i<0 ) {
- elem += row*theVar->value.nelem;
- if( this->type==STRING )
- this->value.undef[row] = theVar->value.undef[row];
- else if( this->type==BITSTR )
- this->value.undef; /* Dummy - BITSTRs do not have undefs */
- else
- this->value.undef[row] = theVar->value.undef[elem];
- if( this->type==DOUBLE )
- this->value.data.dblptr[row] =
- theVar->value.data.dblptr[elem];
- else if( this->type==LONG )
- this->value.data.lngptr[row] =
- theVar->value.data.lngptr[elem];
- else if( this->type==BOOLEAN )
- this->value.data.logptr[row] =
- theVar->value.data.logptr[elem];
- else {
- /* XXX Note, the below expression uses knowledge of
- the layout of the string format, namely (nelem+1)
- characters per string, followed by (nelem+1)
- "undef" values. */
- this->value.data.strptr[row][0] =
- theVar->value.data.strptr[0][elem+row];
- this->value.data.strptr[row][1] = 0; /* Null terminate */
- }
- } else {
- yyerror("Index out of range");
- free( this->value.data.ptr );
- }
- }
- } else {
- /* Reduce dimensions by 1, using a nonconstant expression */
- for( row=0; row<gParse.nRows; row++ ) {
- /* Index cannot be a constant */
- if( theDims[0]->value.undef[row] ) {
- yyerror("Null encountered as vector index");
- free( this->value.data.ptr );
- break;
- } else
- dimVals[0] = theDims[0]->value.data.lngptr[row];
- if( dimVals[0] < 1 ||
- dimVals[0] > theVar->value.naxes[ theVar->value.naxis-1 ] ) {
- yyerror("Index out of range");
- free( this->value.data.ptr );
- } else if ( this->type == BITSTR || this->type == STRING ) {
- elem = this->value.nelem * (dimVals[0]-1);
- elem += row*(theVar->value.nelem+1);
- if (this->value.undef)
- this->value.undef[row] = theVar->value.undef[row];
- memcpy( (char*)this->value.data.strptr[0]
- + row*sizeof(char)*(this->value.nelem+1),
- (char*)theVar->value.data.strptr[0] + elem*sizeof(char),
- this->value.nelem * sizeof(char) );
- /* Null terminate */
- this->value.data.strptr[row][this->value.nelem] = 0;
- } else {
- elem = this->value.nelem * (dimVals[0]-1);
- elem += row*theVar->value.nelem;
- memcpy( this->value.undef + row*this->value.nelem,
- theVar->value.undef + elem,
- this->value.nelem * sizeof(char) );
- memcpy( (char*)this->value.data.ptr
- + row*dsize*this->value.nelem,
- (char*)theVar->value.data.ptr + elem*dsize,
- this->value.nelem * dsize );
- }
- }
- }
- }
- if( theVar->operation>0 ) {
- if (theVar->type == STRING || theVar->type == BITSTR)
- free(theVar->value.data.strptr[0] );
- else
- free( theVar->value.data.ptr );
- }
- for( i=0; i<nDims; i++ )
- if( theDims[i]->operation>0 ) {
- free( theDims[i]->value.data.ptr );
- }
- }
- static void Do_GTI( Node *this )
- {
- Node *theExpr, *theTimes;
- double *start, *stop, *times;
- long elem, nGTI, gti;
- int ordered;
- theTimes = gParse.Nodes + this->SubNodes[0];
- theExpr = gParse.Nodes + this->SubNodes[1];
- nGTI = theTimes->value.nelem;
- start = theTimes->value.data.dblptr;
- stop = theTimes->value.data.dblptr + nGTI;
- ordered = theTimes->type;
- if( theExpr->operation==CONST_OP ) {
- this->value.data.log =
- (Search_GTI( theExpr->value.data.dbl, nGTI, start, stop, ordered )>=0);
- this->operation = CONST_OP;
- } else {
- Allocate_Ptrs( this );
- times = theExpr->value.data.dblptr;
- if( !gParse.status ) {
- elem = gParse.nRows * this->value.nelem;
- if( nGTI ) {
- gti = -1;
- while( elem-- ) {
- if( (this->value.undef[elem] = theExpr->value.undef[elem]) )
- continue;
- /* Before searching entire GTI, check the GTI found last time */
- if( gti<0 || times[elem]<start[gti] || times[elem]>stop[gti] ) {
- gti = Search_GTI( times[elem], nGTI, start, stop, ordered );
- }
- this->value.data.logptr[elem] = ( gti>=0 );
- }
- } else
- while( elem-- ) {
- this->value.data.logptr[elem] = 0;
- this->value.undef[elem] = 0;
- }
- }
- }
- if( theExpr->operation>0 )
- free( theExpr->value.data.ptr );
- }
- static long Search_GTI( double evtTime, long nGTI, double *start,
- double *stop, int ordered )
- {
- long gti, step;
-
- if( ordered && nGTI>15 ) { /* If time-ordered and lots of GTIs, */
- /* use "FAST" Binary search algorithm */
- if( evtTime>=start[0] && evtTime<=stop[nGTI-1] ) {
- gti = step = (nGTI >> 1);
- while(1) {
- if( step>1L ) step >>= 1;
-
- if( evtTime>stop[gti] ) {
- if( evtTime>=start[gti+1] )
- gti += step;
- else {
- gti = -1L;
- break;
- }
- } else if( evtTime<start[gti] ) {
- if( evtTime<=stop[gti-1] )
- gti -= step;
- else {
- gti = -1L;
- break;
- }
- } else {
- break;
- }
- }
- } else
- gti = -1L;
-
- } else { /* Use "SLOW" linear search */
- gti = nGTI;
- while( gti-- )
- if( evtTime>=start[gti] && evtTime<=stop[gti] )
- break;
- }
- return( gti );
- }
- static void Do_REG( Node *this )
- {
- Node *theRegion, *theX, *theY;
- double Xval=0.0, Yval=0.0;
- char Xnull=0, Ynull=0;
- int Xvector, Yvector;
- long nelem, elem, rows;
- theRegion = gParse.Nodes + this->SubNodes[0];
- theX = gParse.Nodes + this->SubNodes[1];
- theY = gParse.Nodes + this->SubNodes[2];
- Xvector = ( theX->operation!=CONST_OP );
- if( Xvector )
- Xvector = theX->value.nelem;
- else {
- Xval = theX->value.data.dbl;
- }
- Yvector = ( theY->operation!=CONST_OP );
- if( Yvector )
- Yvector = theY->value.nelem;
- else {
- Yval = theY->value.data.dbl;
- }
- if( !Xvector && !Yvector ) {
- this->value.data.log =
- ( fits_in_region( Xval, Yval, (SAORegion *)theRegion->value.data.ptr )
- != 0 );
- this->operation = CONST_OP;
- } else {
- Allocate_Ptrs( this );
- if( !gParse.status ) {
- rows = gParse.nRows;
- nelem = this->value.nelem;
- elem = rows*nelem;
- while( rows-- ) {
- while( nelem-- ) {
- elem--;
- if( Xvector>1 ) {
- Xval = theX->value.data.dblptr[elem];
- Xnull = theX->value.undef[elem];
- } else if( Xvector ) {
- Xval = theX->value.data.dblptr[rows];
- Xnull = theX->value.undef[rows];
- }
- if( Yvector>1 ) {
- Yval = theY->value.data.dblptr[elem];
- Ynull = theY->value.undef[elem];
- } else if( Yvector ) {
- Yval = theY->value.data.dblptr[rows];
- Ynull = theY->value.undef[rows];
- }
- this->value.undef[elem] = ( Xnull || Ynull );
- if( this->value.undef[elem] )
- continue;
- this->value.data.logptr[elem] =
- ( fits_in_region( Xval, Yval,
- (SAORegion *)theRegion->value.data.ptr )
- != 0 );
- }
- nelem = this->value.nelem;
- }
- }
- }
- if( theX->operation>0 )
- free( theX->value.data.ptr );
- if( theY->operation>0 )
- free( theY->value.data.ptr );
- }
- static void Do_Vector( Node *this )
- {
- Node *that;
- long row, elem, idx, jdx, offset=0;
- int node;
- Allocate_Ptrs( this );
- if( !gParse.status ) {
- for( node=0; node<this->nSubNodes; node++ ) {
- that = gParse.Nodes + this->SubNodes[node];
- if( that->operation == CONST_OP ) {
- idx = gParse.nRows*this->value.nelem + offset;
- while( (idx-=this->value.nelem)>=0 ) {
-
- this->value.undef[idx] = 0;
- switch( this->type ) {
- case BOOLEAN:
- this->value.data.logptr[idx] = that->value.data.log;
- break;
- case LONG:
- this->value.data.lngptr[idx] = that->value.data.lng;
- break;
- case DOUBLE:
- this->value.data.dblptr[idx] = that->value.data.dbl;
- break;
- }
- }
-
- } else {
-
- row = gParse.nRows;
- idx = row * that->value.nelem;
- while( row-- ) {
- elem = that->value.nelem;
- jdx = row*this->value.nelem + offset;
- while( elem-- ) {
- this->value.undef[jdx+elem] =
- that->value.undef[--idx];
- switch( this->type ) {
- case BOOLEAN:
- this->value.data.logptr[jdx+elem] =
- that->value.data.logptr[idx];
- break;
- case LONG:
- this->value.data.lngptr[jdx+elem] =
- that->value.data.lngptr[idx];
- break;
- case DOUBLE:
- this->value.data.dblptr[jdx+elem] =
- that->value.data.dblptr[idx];
- break;
- }
- }
- }
- }
- offset += that->value.nelem;
- }
- }
- for( node=0; node < this->nSubNodes; node++ )
- if( OPER(this->SubNodes[node])>0 )
- free( gParse.Nodes[this->SubNodes[node]].value.data.ptr );
- }
- /*****************************************************************************/
- /* Utility routines which perform the calculations on bits and SAO regions */
- /*****************************************************************************/
- static char bitlgte(char *bits1, int oper, char *bits2)
- {
- int val1, val2, nextbit;
- char result;
- int i, l1, l2, length, ldiff;
- char stream[256];
- char chr1, chr2;
- l1 = strlen(bits1);
- l2 = strlen(bits2);
- if (l1 < l2)
- {
- length = l2;
- ldiff = l2 - l1;
- i=0;
- while( ldiff-- ) stream[i++] = '0';
- while( l1-- ) stream[i++] = *(bits1++);
- stream[i] = '\0';
- bits1 = stream;
- }
- else if (l2 < l1)
- {
- length = l1;
- ldiff = l1 - l2;
- i=0;
- while( ldiff-- ) stream[i++] = '0';
- while( l2-- ) stream[i++] = *(bits2++);
- stream[i] = '\0';
- bits2 = stream;
- }
- else
- length = l1;
- val1 = val2 = 0;
- nextbit = 1;
- while( length-- )
- {
- chr1 = bits1[length];
- chr2 = bits2[length];
- if ((chr1 != 'x')&&(chr1 != 'X')&&(chr2 != 'x')&&(chr2 != 'X'))
- {
- if (chr1 == '1') val1 += nextbit;
- if (chr2 == '1') val2 += nextbit;
- nextbit *= 2;
- }
- }
- result = 0;
- switch (oper)
- {
- case LT:
- if (val1 < val2) result = 1;
- break;
- case LTE:
- if (val1 <= val2) result = 1;
- break;
- case GT:
- if (val1 > val2) result = 1;
- break;
- case GTE:
- if (val1 >= val2) result = 1;
- break;
- }
- return (result);
- }
- static void bitand(char *result,char *bitstrm1,char *bitstrm2)
- {
- int i, l1, l2, ldiff;
- char stream[256];
- char chr1, chr2;
- l1 = strlen(bitstrm1);
- l2 = strlen(bitstrm2);
- if (l1 < l2)
- {
- ldiff = l2 - l1;
- i=0;
- while( ldiff-- ) stream[i++] = '0';
- while( l1-- ) stream[i++] = *(bitstrm1++);
- stream[i] = '\0';
- bitstrm1 = stream;
- }
- else if (l2 < l1)
- {
- ldiff = l1 - l2;
- i=0;
- while( ldiff-- ) stream[i++] = '0';
- while( l2-- ) stream[i++] = *(bitstrm2++);
- stream[i] = '\0';
- bitstrm2 = stream;
- }
- while ( (chr1 = *(bitstrm1++)) )
- {
- chr2 = *(bitstrm2++);
- if ((chr1 == 'x') || (chr2 == 'x'))
- *result = 'x';
- else if ((chr1 == '1') && (chr2 == '1'))
- *result = '1';
- else
- *result = '0';
- result++;
- }
- *result = '\0';
- }
- static void bitor(char *result,char *bitstrm1,char *bitstrm2)
- {
- int i, l1, l2, ldiff;
- char stream[256];
- char chr1, chr2;
- l1 = strlen(bitstrm1);
- l2 = strlen(bitstrm2);
- if (l1 < l2)
- {
- ldiff = l2 - l1;
- i=0;
- while( ldiff-- ) stream[i++] = '0';
- while( l1-- ) stream[i++] = *(bitstrm1++);
- stream[i] = '\0';
- bitstrm1 = stream;
- }
- else if (l2 < l1)
- {
- ldiff = l1 - l2;
- i=0;
- while( ldiff-- ) stream[i++] = '0';
- while( l2-- ) stream[i++] = *(bitstrm2++);
- stream[i] = '\0';
- bitstrm2 = stream;
- }
- while ( (chr1 = *(bitstrm1++)) )
- {
- chr2 = *(bitstrm2++);
- if ((chr1 == '1') || (chr2 == '1'))
- *result = '1';
- else if ((chr1 == '0') || (chr2 == '0'))
- *result = '0';
- else
- *result = 'x';
- result++;
- }
- *result = '\0';
- }
- static void bitnot(char *result,char *bits)
- {
- int length;
- char chr;
- length = strlen(bits);
- while( length-- ) {
- chr = *(bits++);
- *(result++) = ( chr=='1' ? '0' : ( chr=='0' ? '1' : chr ) );
- }
- *result = '\0';
- }
- static char bitcmp(char *bitstrm1, char *bitstrm2)
- {
- int i, l1, l2, ldiff;
- char stream[256];
- char chr1, chr2;
- l1 = strlen(bitstrm1);
- l2 = strlen(bitstrm2);
- if (l1 < l2)
- {
- ldiff = l2 - l1;
- i=0;
- while( ldiff-- ) stream[i++] = '0';
- while( l1-- ) stream[i++] = *(bitstrm1++);
- stream[i] = '\0';
- bitstrm1 = stream;
- }
- else if (l2 < l1)
- {
- ldiff = l1 - l2;
- i=0;
- while( ldiff-- ) stream[i++] = '0';
- while( l2-- ) stream[i++] = *(bitstrm2++);
- stream[i] = '\0';
- bitstrm2 = stream;
- }
- while( (chr1 = *(bitstrm1++)) )
- {
- chr2 = *(bitstrm2++);
- if ( ((chr1 == '0') && (chr2 == '1'))
- || ((chr1 == '1') && (chr2 == '0')) )
- return( 0 );
- }
- return( 1 );
- }
- static char bnear(double x, double y, double tolerance)
- {
- if (fabs(x - y) < tolerance)
- return ( 1 );
- else
- return ( 0 );
- }
- static char saobox(double xcen, double ycen, double xwid, double ywid,
- double rot, double xcol, double ycol)
- {
- double x,y,xprime,yprime,xmin,xmax,ymin,ymax,theta;
- theta = (rot / 180.0) * myPI;
- xprime = xcol - xcen;
- yprime = ycol - ycen;
- x = xprime * cos(theta) + yprime * sin(theta);
- y = -xprime * sin(theta) + yprime * cos(theta);
- xmin = - 0.5 * xwid; xmax = 0.5 * xwid;
- ymin = - 0.5 * ywid; ymax = 0.5 * ywid;
- if ((x >= xmin) && (x <= xmax) && (y >= ymin) && (y <= ymax))
- return ( 1 );
- else
- return ( 0 );
- }
- static char circle(double xcen, double ycen, double rad,
- double xcol, double ycol)
- {
- double r2,dx,dy,dlen;
- dx = xcol - xcen;
- dy = ycol - ycen;
- dx *= dx; dy *= dy;
- dlen = dx + dy;
- r2 = rad * rad;
- if (dlen <= r2)
- return ( 1 );
- else
- return ( 0 );
- }
- static char ellipse(double xcen, double ycen, double xrad, double yrad,
- double rot, double xcol, double ycol)
- {
- double x,y,xprime,yprime,dx,dy,dlen,theta;
- theta = (rot / 180.0) * myPI;
- xprime = xcol - xcen;
- yprime = ycol - ycen;
- x = xprime * cos(theta) + yprime * sin(theta);
- y = -xprime * sin(theta) + yprime * cos(theta);
- dx = x / xrad; dy = y / yrad;
- dx *= dx; dy *= dy;
- dlen = dx + dy;
- if (dlen <= 1.0)
- return ( 1 );
- else
- return ( 0 );
- }
- /*
- * Extract substring
- */
- int cstrmid(char *dest_str, int dest_len,
- char *src_str, int src_len,
- int pos)
- {
- /* char fill_char = ' '; */
- char fill_char = '\0';
- if (src_len == 0) { src_len = strlen(src_str); } /* .. if constant */
- /* Fill destination with blanks */
- if (pos < 0) {
- yyerror("STRMID(S,P,N) P must be 0 or greater");
- return -1;
- }
- if (pos > src_len || pos == 0) {
- /* pos==0: blank string requested */
- memset(dest_str, fill_char, dest_len);
- } else if (pos+dest_len > src_len) {
- /* Copy a subset */
- int nsub = src_len-pos+1;
- int npad = dest_len - nsub;
- memcpy(dest_str, src_str+pos-1, nsub);
- /* Fill remaining string with blanks */
- memset(dest_str+nsub, fill_char, npad);
- } else {
- /* Full string copy */
- memcpy(dest_str, src_str+pos-1, dest_len);
- }
- dest_str[dest_len] = '\0'; /* Null-terminate */
- return 0;
- }
- static void yyerror(char *s)
- {
- char msg[80];
- if( !gParse.status ) gParse.status = PARSE_SYNTAX_ERR;
- strncpy(msg, s, 80);
- msg[79] = '\0';
- ffpmsg(msg);
- }