"Fossies" - the Fresh Open Source Software Archive

Member "cfitsio-4.0.0/eval_f.c" (13 Jul 2021, 99081 Bytes) of package /linux/misc/cfitsio-4.0.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file.

    1 /************************************************************************/
    2 /*                                                                      */
    3 /*                       CFITSIO Lexical Parser                         */
    4 /*                                                                      */
    5 /* This file is one of 3 files containing code which parses an          */
    6 /* arithmetic expression and evaluates it in the context of an input    */
    7 /* FITS file table extension.  The CFITSIO lexical parser is divided    */
    8 /* into the following 3 parts/files: the CFITSIO "front-end",           */
    9 /* eval_f.c, contains the interface between the user/CFITSIO and the    */
   10 /* real core of the parser; the FLEX interpreter, eval_l.c, takes the   */
   11 /* input string and parses it into tokens and identifies the FITS       */
   12 /* information required to evaluate the expression (ie, keywords and    */
   13 /* columns); and, the BISON grammar and evaluation routines, eval_y.c,  */
   14 /* receives the FLEX output and determines and performs the actual      */
   15 /* operations.  The files eval_l.c and eval_y.c are produced from       */
   16 /* running flex and bison on the files eval.l and eval.y, respectively. */
   17 /* (flex and bison are available from any GNU archive: see www.gnu.org) */
   18 /*                                                                      */
   19 /* The grammar rules, rather than evaluating the expression in situ,    */
   20 /* builds a tree, or Nodal, structure mapping out the order of          */
   21 /* operations and expression dependencies.  This "compilation" process  */
   22 /* allows for much faster processing of multiple rows.  This technique  */
   23 /* was developed by Uwe Lammers of the XMM Science Analysis System,     */
   24 /* although the CFITSIO implementation is entirely code original.       */
   25 /*                                                                      */
   26 /*                                                                      */
   27 /* Modification History:                                                */
   28 /*                                                                      */
   29 /*   Kent Blackburn      c1992  Original parser code developed for the  */
   30 /*                              FTOOLS software package, in particular, */
   31 /*                              the fselect task.                       */
   32 /*   Kent Blackburn      c1995  BIT column support added                */
   33 /*   Peter D Wilson   Feb 1998  Vector column support added             */
   34 /*   Peter D Wilson   May 1998  Ported to CFITSIO library.  User        */
   35 /*                              interface routines written, in essence  */
   36 /*                              making fselect, fcalc, and maketime     */
   37 /*                              capabilities available to all tools     */
   38 /*                              via single function calls.              */
   39 /*   Peter D Wilson   Jun 1998  Major rewrite of parser core, so as to  */
   40 /*                              create a run-time evaluation tree,      */
   41 /*                              inspired by the work of Uwe Lammers,    */
   42 /*                              resulting in a speed increase of        */
   43 /*                              10-100 times.                           */
   44 /*   Peter D Wilson   Jul 1998  gtifilter(a,b,c,d) function added       */
   45 /*   Peter D Wilson   Aug 1998  regfilter(a,b,c,d) function added       */
   46 /*   Peter D Wilson   Jul 1999  Make parser fitsfile-independent,       */
   47 /*                              allowing a purely vector-based usage    */
   48 /*   Peter D Wilson   Aug 1999  Add row-offset capability               */
   49 /*   Peter D Wilson   Sep 1999  Add row-range capability to ffcalc_rng  */
   50 /*                                                                      */
   51 /************************************************************************/
   52 
   53 #include <limits.h>
   54 #include <ctype.h>
   55 #include "eval_defs.h"
   56 #include "region.h"
   57 
   58 typedef struct {
   59      int  datatype;   /* Data type to cast parse results into for user       */
   60      void *dataPtr;   /* Pointer to array of results, NULL if to use iterCol */
   61      void *nullPtr;   /* Pointer to nulval, use zero if NULL                 */
   62      long maxRows;    /* Max No. of rows to process, -1=all, 0=1 iteration   */
   63      int  anyNull;    /* Flag indicating at least 1 undef value encountered  */
   64 } parseInfo;
   65 
   66 /*  Internal routines needed to allow the evaluator to operate on FITS data  */
   67 
   68 static void Setup_DataArrays( int nCols, iteratorCol *cols,
   69                               long fRow, long nRows );
   70 static int  find_column( char *colName, void *itslval );
   71 static int  find_keywd ( char *key,     void *itslval );
   72 static int  allocateCol( int nCol, int *status );
   73 static int  load_column( int varNum, long fRow, long nRows,
   74                          void *data, char *undef );
   75 
   76 static int DEBUG_PIXFILTER;
   77 
   78 #define FREE(x) { if (x) free(x); else printf("invalid free(" #x ") at %s:%d\n", __FILE__, __LINE__); }
   79 
   80 /*---------------------------------------------------------------------------*/
   81 int fffrow( fitsfile *fptr,         /* I - Input FITS file                   */
   82             char     *expr,         /* I - Boolean expression                */
   83             long     firstrow,      /* I - First row of table to eval        */
   84             long     nrows,         /* I - Number of rows to evaluate        */
   85             long     *n_good_rows,  /* O - Number of rows eval to True       */
   86             char     *row_status,   /* O - Array of boolean results          */
   87             int      *status )      /* O - Error status                      */
   88 /*                                                                           */
   89 /* Evaluate a boolean expression using the indicated rows, returning an      */
   90 /* array of flags indicating which rows evaluated to TRUE/FALSE              */
   91 /*---------------------------------------------------------------------------*/
   92 {
   93    parseInfo Info;
   94    int naxis, constant;
   95    long nelem, naxes[MAXDIMS], elem;
   96    char result;
   97 
   98    if( *status ) return( *status );
   99 
  100    FFLOCK;
  101    if( ffiprs( fptr, 0, expr, MAXDIMS, &Info.datatype, &nelem, &naxis,
  102                naxes, status ) ) {
  103       ffcprs();
  104       FFUNLOCK;
  105       return( *status );
  106    }
  107    if( nelem<0 ) {
  108       constant = 1;
  109       nelem = -nelem;
  110    } else
  111       constant = 0;
  112 
  113    if( Info.datatype!=TLOGICAL || nelem!=1 ) {
  114       ffcprs();
  115       ffpmsg("Expression does not evaluate to a logical scalar.");
  116       FFUNLOCK;
  117       return( *status = PARSE_BAD_TYPE );
  118    }
  119 
  120    if( constant ) { /* No need to call parser... have result from ffiprs */
  121       result = gParse.Nodes[gParse.resultNode].value.data.log;
  122       *n_good_rows = nrows;
  123       for( elem=0; elem<nrows; elem++ )
  124          row_status[elem] = result;
  125    } else {
  126       firstrow     = (firstrow>1 ? firstrow : 1);
  127       Info.dataPtr = row_status;
  128       Info.nullPtr = NULL;
  129       Info.maxRows = nrows;
  130 
  131       if( ffiter( gParse.nCols, gParse.colData, firstrow-1, 0,
  132                   parse_data, (void*)&Info, status ) == -1 )
  133          *status = 0;  /* -1 indicates exitted without error before end... OK */
  134 
  135       if( *status ) {
  136 
  137          /***********************/
  138          /* Error... Do nothing */
  139          /***********************/
  140 
  141       } else {
  142 
  143          /***********************************/
  144          /* Count number of good rows found */
  145          /***********************************/
  146 
  147          *n_good_rows = 0L;
  148          for( elem=0; elem<Info.maxRows; elem++ ) {
  149             if( row_status[elem]==1 ) ++*n_good_rows;
  150          }
  151       }
  152    }
  153 
  154    ffcprs();
  155    FFUNLOCK;
  156    return(*status);
  157 }
  158 
  159 /*--------------------------------------------------------------------------*/
  160 int ffsrow( fitsfile *infptr,   /* I - Input FITS file                      */
  161             fitsfile *outfptr,  /* I - Output FITS file                     */
  162             char     *expr,     /* I - Boolean expression                   */
  163             int      *status )  /* O - Error status                         */
  164 /*                                                                          */
  165 /* Evaluate an expression on all rows of a table.  If the input and output  */
  166 /* files are not the same, copy the TRUE rows to the output file.  If the   */
  167 /* files are the same, delete the FALSE rows (preserve the TRUE rows).      */
  168 /* Can copy rows between extensions of the same file, *BUT* if output       */
  169 /* extension is before the input extension, the second extension *MUST* be  */
  170 /* opened using ffreopen, so that CFITSIO can handle changing file lengths. */
  171 /*--------------------------------------------------------------------------*/
  172 {
  173    parseInfo Info;
  174    int naxis, constant;
  175    long nelem, rdlen, naxes[MAXDIMS], maxrows, nbuff, nGood, inloc, outloc;
  176    LONGLONG ntodo, inbyteloc, outbyteloc, hsize;
  177    long freespace;
  178    unsigned char *buffer, result;
  179    struct {
  180       LONGLONG rowLength, numRows, heapSize;
  181       LONGLONG dataStart, heapStart;
  182    } inExt, outExt;
  183 
  184    if( *status ) return( *status );
  185 
  186    FFLOCK;
  187    if( ffiprs( infptr, 0, expr, MAXDIMS, &Info.datatype, &nelem, &naxis,
  188                naxes, status ) ) {
  189       ffcprs();
  190       FFUNLOCK;
  191       return( *status );
  192    }
  193 
  194    if( nelem<0 ) {
  195       constant = 1;
  196       nelem = -nelem;
  197    } else
  198       constant = 0;
  199 
  200    /**********************************************************************/
  201    /* Make sure expression evaluates to the right type... logical scalar */
  202    /**********************************************************************/
  203 
  204    if( Info.datatype!=TLOGICAL || nelem!=1 ) {
  205       ffcprs();
  206       ffpmsg("Expression does not evaluate to a logical scalar.");
  207       FFUNLOCK;
  208       return( *status = PARSE_BAD_TYPE );
  209    }
  210 
  211    /***********************************************************/
  212    /*  Extract various table information from each extension  */
  213    /***********************************************************/
  214 
  215    if( infptr->HDUposition != (infptr->Fptr)->curhdu )
  216       ffmahd( infptr, (infptr->HDUposition) + 1, NULL, status );
  217    if( *status ) {
  218       ffcprs();
  219       FFUNLOCK;
  220       return( *status );
  221    }
  222    inExt.rowLength = (long) (infptr->Fptr)->rowlength;
  223    inExt.numRows   = (infptr->Fptr)->numrows;
  224    inExt.heapSize  = (infptr->Fptr)->heapsize;
  225    if( inExt.numRows == 0 ) { /* Nothing to copy */
  226       ffcprs();
  227       FFUNLOCK;
  228       return( *status );
  229    }
  230 
  231    if( outfptr->HDUposition != (outfptr->Fptr)->curhdu )
  232       ffmahd( outfptr, (outfptr->HDUposition) + 1, NULL, status );
  233    if( (outfptr->Fptr)->datastart < 0 )
  234       ffrdef( outfptr, status );
  235    if( *status ) {
  236       ffcprs();
  237       FFUNLOCK;
  238       return( *status );
  239    }
  240    outExt.rowLength = (long) (outfptr->Fptr)->rowlength;
  241    outExt.numRows   = (outfptr->Fptr)->numrows;
  242    if( !outExt.numRows )
  243       (outfptr->Fptr)->heapsize = 0L;
  244    outExt.heapSize  = (outfptr->Fptr)->heapsize;
  245 
  246    if( inExt.rowLength != outExt.rowLength ) {
  247       ffpmsg("Output table has different row length from input");
  248       ffcprs();
  249       FFUNLOCK;
  250       return( *status = PARSE_BAD_OUTPUT );
  251    }
  252 
  253    /***********************************/
  254    /*  Fill out Info data for parser  */
  255    /***********************************/
  256 
  257    Info.dataPtr = (char *)malloc( (size_t) ((inExt.numRows + 1) * sizeof(char)) );
  258    Info.nullPtr = NULL;
  259    Info.maxRows = (long) inExt.numRows;
  260    if( !Info.dataPtr ) {
  261       ffpmsg("Unable to allocate memory for row selection");
  262       ffcprs();
  263       FFUNLOCK;
  264       return( *status = MEMORY_ALLOCATION );
  265    }
  266    
  267    /* make sure array is zero terminated */
  268    ((char*)Info.dataPtr)[inExt.numRows] = 0;
  269 
  270    if( constant ) { /*  Set all rows to the same value from constant result  */
  271 
  272       result = gParse.Nodes[gParse.resultNode].value.data.log;
  273       for( ntodo = 0; ntodo<inExt.numRows; ntodo++ )
  274          ((char*)Info.dataPtr)[ntodo] = result;
  275       nGood = (long) (result ? inExt.numRows : 0);
  276 
  277    } else {
  278 
  279       ffiter( gParse.nCols, gParse.colData, 0L, 0L,
  280               parse_data, (void*)&Info, status );
  281 
  282       nGood = 0;
  283       for( ntodo = 0; ntodo<inExt.numRows; ntodo++ )
  284          if( ((char*)Info.dataPtr)[ntodo] ) nGood++;
  285    }
  286 
  287    if( *status ) {
  288       /* Error... Do nothing */
  289    } else {
  290       rdlen  = (long) inExt.rowLength;
  291       buffer = (unsigned char *)malloc(maxvalue(500000,rdlen) * sizeof(char) );
  292       if( buffer==NULL ) {
  293          ffcprs();
  294          FFUNLOCK;
  295          return( *status=MEMORY_ALLOCATION );
  296       }
  297       maxrows = maxvalue( (500000L/rdlen), 1);
  298       nbuff = 0;
  299       inloc = 1;
  300       if( infptr==outfptr ) { /* Skip initial good rows if input==output file */
  301          while( ((char*)Info.dataPtr)[inloc-1] ) inloc++;
  302          outloc = inloc;
  303       } else {
  304          outloc = (long) (outExt.numRows + 1);
  305          if (outloc > 1) 
  306             ffirow( outfptr, outExt.numRows, nGood, status );
  307       }
  308 
  309       do {
  310          if( ((char*)Info.dataPtr)[inloc-1] ) {
  311             ffgtbb( infptr, inloc, 1L, rdlen, buffer+rdlen*nbuff, status );
  312             nbuff++;
  313             if( nbuff==maxrows ) {
  314                ffptbb( outfptr, outloc, 1L, rdlen*nbuff, buffer,  status );
  315                outloc += nbuff;
  316                nbuff = 0;
  317             }
  318          }
  319          inloc++;
  320       } while( !*status && inloc<=inExt.numRows );
  321 
  322       if( nbuff ) {
  323          ffptbb( outfptr, outloc, 1L, rdlen*nbuff, buffer,  status );
  324          outloc += nbuff;
  325       }
  326 
  327       if( infptr==outfptr ) {
  328 
  329          if( outloc<=inExt.numRows )
  330             ffdrow( infptr, outloc, inExt.numRows-outloc+1, status );
  331 
  332       } else if( inExt.heapSize && nGood ) {
  333 
  334          /* Copy heap, if it exists and at least one row copied */
  335 
  336          /********************************************************/
  337          /*  Get location information from the output extension  */
  338          /********************************************************/
  339 
  340          if( outfptr->HDUposition != (outfptr->Fptr)->curhdu )
  341             ffmahd( outfptr, (outfptr->HDUposition) + 1, NULL, status );
  342          outExt.dataStart = (outfptr->Fptr)->datastart;
  343          outExt.heapStart = (outfptr->Fptr)->heapstart;
  344 
  345          /*************************************************/
  346          /*  Insert more space into outfptr if necessary  */
  347          /*************************************************/
  348 
  349          hsize     = outExt.heapStart + outExt.heapSize;
  350          freespace = (long) (( ( (hsize + 2879) / 2880) * 2880) - hsize);
  351          ntodo     = inExt.heapSize;
  352 
  353          if ( (freespace - ntodo) < 0) {       /* not enough existing space? */
  354             ntodo = (ntodo - freespace + 2879) / 2880;  /* number of blocks  */
  355             ffiblk(outfptr, (long) ntodo, 1, status);   /* insert the blocks */
  356          }
  357          ffukyj( outfptr, "PCOUNT", inExt.heapSize+outExt.heapSize,
  358                  NULL, status );
  359 
  360          /*******************************************************/
  361          /*  Get location information from the input extension  */
  362          /*******************************************************/
  363 
  364          if( infptr->HDUposition != (infptr->Fptr)->curhdu )
  365             ffmahd( infptr, (infptr->HDUposition) + 1, NULL, status );
  366          inExt.dataStart = (infptr->Fptr)->datastart;
  367          inExt.heapStart = (infptr->Fptr)->heapstart;
  368 
  369          /**********************************/
  370          /*  Finally copy heap to outfptr  */
  371          /**********************************/
  372 
  373          ntodo  =  inExt.heapSize;
  374          inbyteloc  =  inExt.heapStart +  inExt.dataStart;
  375          outbyteloc = outExt.heapStart + outExt.dataStart + outExt.heapSize;
  376 
  377          while ( ntodo && !*status ) {
  378             rdlen = (long) minvalue(ntodo,500000);
  379             ffmbyt( infptr,  inbyteloc,  REPORT_EOF, status );
  380             ffgbyt( infptr,  rdlen,  buffer,     status );
  381             ffmbyt( outfptr, outbyteloc, IGNORE_EOF, status );
  382             ffpbyt( outfptr, rdlen,  buffer,     status );
  383             inbyteloc  += rdlen;
  384             outbyteloc += rdlen;
  385             ntodo  -= rdlen;
  386          }
  387 
  388          /***********************************************************/
  389          /*  But must update DES if data is being appended to a     */
  390          /*  pre-existing heap space.  Edit each new entry in file  */
  391          /***********************************************************/
  392 
  393          if( outExt.heapSize ) {
  394             LONGLONG repeat, offset, j;
  395             int i;
  396             for( i=1; i<=(outfptr->Fptr)->tfield; i++ ) {
  397                if( (outfptr->Fptr)->tableptr[i-1].tdatatype<0 ) {
  398                   for( j=outExt.numRows+1; j<=outExt.numRows+nGood; j++ ) {
  399                      ffgdesll( outfptr, i, j, &repeat, &offset, status );
  400                      offset += outExt.heapSize;
  401                      ffpdes( outfptr, i, j, repeat, offset, status );
  402                   }
  403                }
  404             }
  405          }
  406 
  407       } /*  End of HEAP copy  */
  408 
  409       FREE(buffer);
  410    }
  411 
  412    FREE(Info.dataPtr);
  413    ffcprs();
  414 
  415    ffcmph(outfptr, status);  /* compress heap, deleting any orphaned data */
  416    FFUNLOCK;
  417    return(*status);
  418 }
  419 
  420 /*---------------------------------------------------------------------------*/
  421 int ffcrow( fitsfile *fptr,      /* I - Input FITS file                      */
  422             int      datatype,   /* I - Datatype to return results as        */
  423             char     *expr,      /* I - Arithmetic expression                */
  424             long     firstrow,   /* I - First row to evaluate                */
  425             long     nelements,  /* I - Number of elements to return         */
  426             void     *nulval,    /* I - Ptr to value to use as UNDEF         */
  427             void     *array,     /* O - Array of results                     */
  428             int      *anynul,    /* O - Were any UNDEFs encountered?         */
  429             int      *status )   /* O - Error status                         */
  430 /*                                                                           */
  431 /* Calculate an expression for the indicated rows of a table, returning      */
  432 /* the results, cast as datatype (TSHORT, TDOUBLE, etc), in array.  If       */
  433 /* nulval==NULL, UNDEFs will be zeroed out.  For vector results, the number  */
  434 /* of elements returned may be less than nelements if nelements is not an    */
  435 /* even multiple of the result dimension.  Call fftexp to obtain the         */
  436 /* dimensions of the results.                                                */
  437 /*---------------------------------------------------------------------------*/
  438 {
  439    parseInfo Info;
  440    int naxis;
  441    long nelem1, naxes[MAXDIMS];
  442 
  443    if( *status ) return( *status );
  444 
  445    FFLOCK;
  446    if( ffiprs( fptr, 0, expr, MAXDIMS, &Info.datatype, &nelem1, &naxis,
  447                naxes, status ) ) {
  448       ffcprs();
  449       FFUNLOCK;
  450       return( *status );
  451    }
  452    if( nelem1<0 ) nelem1 = - nelem1;
  453 
  454    if( nelements<nelem1 ) {
  455       ffcprs();
  456       ffpmsg("Array not large enough to hold at least one row of data.");
  457       FFUNLOCK;
  458       return( *status = PARSE_LRG_VECTOR );
  459    }
  460 
  461    firstrow = (firstrow>1 ? firstrow : 1);
  462 
  463    if( datatype ) Info.datatype = datatype;
  464 
  465    Info.dataPtr = array;
  466    Info.nullPtr = nulval;
  467    Info.maxRows = nelements / nelem1;
  468    
  469    if( ffiter( gParse.nCols, gParse.colData, firstrow-1, 0,
  470                parse_data, (void*)&Info, status ) == -1 )
  471       *status=0;  /* -1 indicates exitted without error before end... OK */
  472 
  473    *anynul = Info.anyNull;
  474    ffcprs();
  475    FFUNLOCK;
  476    return( *status );
  477 }
  478 
  479 /*--------------------------------------------------------------------------*/
  480 int ffcalc( fitsfile *infptr,   /* I - Input FITS file                      */
  481             char     *expr,     /* I - Arithmetic expression                */
  482             fitsfile *outfptr,  /* I - Output fits file                     */
  483             char     *parName,  /* I - Name of output parameter             */
  484             char     *parInfo,  /* I - Extra information on parameter       */
  485             int      *status )  /* O - Error status                         */
  486 /*                                                                          */
  487 /* Evaluate an expression for all rows of a table.  Call ffcalc_rng with    */
  488 /* a row range of 1-MAX.                                                    */
  489 {
  490    long start=1, end=LONG_MAX;
  491 
  492    return ffcalc_rng( infptr, expr, outfptr, parName, parInfo,
  493                       1, &start, &end, status );
  494 }
  495 
  496 /*--------------------------------------------------------------------------*/
  497 int ffcalc_rng( fitsfile *infptr,   /* I - Input FITS file                  */
  498                 char     *expr,     /* I - Arithmetic expression            */
  499                 fitsfile *outfptr,  /* I - Output fits file                 */
  500                 char     *parName,  /* I - Name of output parameter         */
  501                 char     *parInfo,  /* I - Extra information on parameter   */
  502                 int      nRngs,     /* I - Row range info                   */
  503                 long     *start,    /* I - Row range info                   */
  504                 long     *end,      /* I - Row range info                   */
  505                 int      *status )  /* O - Error status                     */
  506 /*                                                                          */
  507 /* Evaluate an expression using the data in the input FITS file and place   */
  508 /* the results into either a column or keyword in the output fits file,     */
  509 /* depending on the value of parName (keywords normally prefixed with '#')  */
  510 /* and whether the expression evaluates to a constant or a table column.    */
  511 /* The logic is as follows:                                                 */
  512 /*    (1) If a column exists with name, parName, put results there.         */
  513 /*    (2) If parName starts with '#', as in #NAXIS, put result there,       */
  514 /*        with parInfo used as the comment. If expression does not evaluate */
  515 /*        to a constant, flag an error.                                     */
  516 /*    (3) If a keyword exists with name, parName, and expression is a       */
  517 /*        constant, put result there, using parInfo as the new comment.     */
  518 /*    (4) Else, create a new column with name parName and TFORM parInfo.    */
  519 /*        If parInfo is NULL, use a default data type for the column.       */
  520 /*--------------------------------------------------------------------------*/
  521 {
  522    parseInfo Info;
  523    int naxis, constant, typecode, newNullKwd=0;
  524    long nelem, naxes[MAXDIMS], repeat, width;
  525    int col_cnt, colNo;
  526    Node *result;
  527    char card[81], tform[16], nullKwd[9], tdimKwd[9];
  528 
  529    if( *status ) return( *status );
  530 
  531    FFLOCK;
  532    if( ffiprs( infptr, 0, expr, MAXDIMS, &Info.datatype, &nelem, &naxis,
  533                naxes, status ) ) {
  534 
  535       ffcprs();
  536       FFUNLOCK;
  537       return( *status );
  538    }
  539    if( nelem<0 ) {
  540       constant = 1;
  541       nelem = -nelem;
  542    } else
  543       constant = 0;
  544 
  545    /*  Case (1): If column exists put it there  */
  546 
  547    colNo = 0;
  548    if( ffgcno( outfptr, CASEINSEN, parName, &colNo, status )==COL_NOT_FOUND ) {
  549 
  550       /*  Output column doesn't exist.  Test for keyword. */
  551 
  552       /* Case (2): Does parName indicate result should be put into keyword */
  553 
  554       *status = 0;
  555       if( parName[0]=='#' ) {
  556          if( ! constant ) {
  557             ffcprs();
  558             ffpmsg( "Cannot put tabular result into keyword (ffcalc)" );
  559             FFUNLOCK;
  560             return( *status = PARSE_BAD_TYPE );
  561          }
  562          parName++;  /* Advance past '#' */
  563      if ( (fits_strcasecmp(parName,"HISTORY") == 0 || fits_strcasecmp(parName,"COMMENT") == 0) &&
  564           Info.datatype != TSTRING ) {
  565             ffcprs();
  566             ffpmsg( "HISTORY and COMMENT values must be strings (ffcalc)" );
  567         FFUNLOCK;
  568         return( *status = PARSE_BAD_TYPE );
  569      }
  570 
  571       } else if( constant ) {
  572 
  573          /* Case (3): Does a keyword named parName already exist */
  574 
  575          if( ffgcrd( outfptr, parName, card, status )==KEY_NO_EXIST ) {
  576             colNo = -1;
  577          } else if( *status ) {
  578             ffcprs();
  579             FFUNLOCK;
  580             return( *status );
  581          }
  582 
  583       } else
  584          colNo = -1;
  585 
  586       if( colNo<0 ) {
  587 
  588          /* Case (4): Create new column */
  589 
  590          *status = 0;
  591          ffgncl( outfptr, &colNo, status );
  592          colNo++;
  593          if( parInfo==NULL || *parInfo=='\0' ) {
  594             /*  Figure out best default column type  */
  595             if( gParse.hdutype==BINARY_TBL ) {
  596                snprintf(tform,15,"%ld",nelem);
  597                switch( Info.datatype ) {
  598                case TLOGICAL:  strcat(tform,"L");  break;
  599                case TLONG:     strcat(tform,"J");  break;
  600                case TDOUBLE:   strcat(tform,"D");  break;
  601                case TSTRING:   strcat(tform,"A");  break;
  602                case TBIT:      strcat(tform,"X");  break;
  603                case TLONGLONG: strcat(tform,"K");  break;
  604                }
  605             } else {
  606                switch( Info.datatype ) {
  607                case TLOGICAL:
  608                   ffcprs();
  609                   ffpmsg("Cannot create LOGICAL column in ASCII table");
  610                   FFUNLOCK;
  611                   return( *status = NOT_BTABLE );
  612                case TLONG:     strcpy(tform,"I11");     break;
  613                case TDOUBLE:   strcpy(tform,"D23.15");  break;
  614                case TSTRING:   
  615                case TBIT:      snprintf(tform,16,"A%ld",nelem);  break;
  616                }
  617             }
  618             parInfo = tform;
  619          } else if( !(isdigit((int) *parInfo)) && gParse.hdutype==BINARY_TBL ) {
  620             if( Info.datatype==TBIT && *parInfo=='B' )
  621                nelem = (nelem+7)/8;
  622             snprintf(tform,16,"%ld%s",nelem,parInfo);
  623             parInfo = tform;
  624          }
  625          fficol( outfptr, colNo, parName, parInfo, status );
  626          if( naxis>1 )
  627             ffptdm( outfptr, colNo, naxis, naxes, status );
  628 
  629          /*  Setup TNULLn keyword in case NULLs are encountered  */
  630 
  631          ffkeyn("TNULL", colNo, nullKwd, status);
  632          if( ffgcrd( outfptr, nullKwd, card, status )==KEY_NO_EXIST ) {
  633             *status = 0;
  634             if( gParse.hdutype==BINARY_TBL ) {
  635            LONGLONG nullVal=0;
  636                fits_binary_tform( parInfo, &typecode, &repeat, &width, status );
  637                if( typecode==TBYTE )
  638                   nullVal = UCHAR_MAX;
  639                else if( typecode==TSHORT )
  640                   nullVal = SHRT_MIN;
  641                else if( typecode==TINT )
  642                   nullVal = INT_MIN;
  643                else if( typecode==TLONG ) {
  644                   if (sizeof(long) == 8 && sizeof(int) == 4)
  645                      nullVal = INT_MIN;
  646                   else
  647                      nullVal = LONG_MIN;
  648                }
  649                else if( typecode==TLONGLONG )
  650                   nullVal = LONGLONG_MIN;
  651           
  652                if( nullVal ) {
  653                   ffpkyj( outfptr, nullKwd, nullVal, "Null value", status );
  654                   fits_set_btblnull( outfptr, colNo, nullVal, status );
  655                   newNullKwd = 1;
  656                }
  657             } else if( gParse.hdutype==ASCII_TBL ) {
  658                ffpkys( outfptr, nullKwd, "NULL", "Null value string", status );
  659                fits_set_atblnull( outfptr, colNo, "NULL", status );
  660                newNullKwd = 1;
  661             }
  662          }
  663 
  664       }
  665 
  666    } else if( *status ) {
  667       ffcprs();
  668       FFUNLOCK;
  669       return( *status );
  670    } else {
  671 
  672       /********************************************************/
  673       /*  Check if a TDIM keyword should be written/updated.  */
  674       /********************************************************/
  675 
  676       ffkeyn("TDIM", colNo, tdimKwd, status);
  677       ffgcrd( outfptr, tdimKwd, card, status );
  678       if( *status==0 ) {
  679          /*  TDIM exists, so update it with result's dimension  */
  680          ffptdm( outfptr, colNo, naxis, naxes, status );
  681       } else if( *status==KEY_NO_EXIST ) {
  682          /*  TDIM does not exist, so clear error stack and     */
  683          /*  write a TDIM only if result is multi-dimensional  */
  684          *status = 0;
  685          ffcmsg();
  686          if( naxis>1 )
  687             ffptdm( outfptr, colNo, naxis, naxes, status );
  688       }
  689       if( *status ) {
  690          /*  Either some other error happened in ffgcrd   */
  691          /*  or one happened in ffptdm                    */
  692          ffcprs();
  693          FFUNLOCK;
  694          return( *status );
  695       }
  696 
  697    }
  698 
  699    if( colNo>0 ) {
  700 
  701       /*  Output column exists (now)... put results into it  */
  702 
  703       int anyNull = 0;
  704       int nPerLp, i;
  705       long totaln;
  706 
  707       ffgkyj(infptr, "NAXIS2", &totaln, 0, status);
  708 
  709       /*************************************/
  710       /* Create new iterator Output Column */
  711       /*************************************/
  712 
  713       col_cnt = gParse.nCols;
  714       if( allocateCol( col_cnt, status ) ) {
  715          ffcprs();
  716          FFUNLOCK;
  717          return( *status );
  718       }
  719 
  720       fits_iter_set_by_num( gParse.colData+col_cnt, outfptr,
  721                             colNo, 0, OutputCol );
  722       gParse.nCols++;
  723 
  724       for( i=0; i<nRngs; i++ ) {
  725          Info.dataPtr = NULL;
  726          Info.maxRows = end[i]-start[i]+1;
  727 
  728           /*
  729             If there is only 1 range, and it includes all the rows,
  730             and there are 10 or more rows, then set nPerLp = 0 so
  731             that the iterator function will dynamically choose the
  732             most efficient number of rows to process in each loop.
  733             Otherwise, set nPerLp to the number of rows in this range.
  734          */
  735 
  736          if( (Info.maxRows >= 10) && (nRngs == 1) &&
  737              (start[0] == 1) && (end[0] == totaln))
  738               nPerLp = 0;
  739          else
  740               nPerLp = Info.maxRows;
  741 
  742          if( ffiter( gParse.nCols, gParse.colData, start[i]-1,
  743                      nPerLp, parse_data, (void*)&Info, status ) == -1 )
  744             *status = 0;
  745          else if( *status ) {
  746             ffcprs();
  747             FFUNLOCK;
  748             return( *status );
  749          }
  750          if( Info.anyNull ) anyNull = 1;
  751       }
  752 
  753       if( newNullKwd && !anyNull ) {
  754          ffdkey( outfptr, nullKwd, status );
  755       }
  756 
  757    } else {
  758 
  759       /* Put constant result into keyword */
  760 
  761       result  = gParse.Nodes + gParse.resultNode;
  762       switch( Info.datatype ) {
  763       case TDOUBLE:
  764          ffukyd( outfptr, parName, result->value.data.dbl, 15,
  765                  parInfo, status );
  766          break;
  767       case TLONG:
  768          ffukyj( outfptr, parName, result->value.data.lng, parInfo, status );
  769          break;
  770       case TLOGICAL:
  771          ffukyl( outfptr, parName, result->value.data.log, parInfo, status );
  772          break;
  773       case TBIT:
  774       case TSTRING:
  775      if (fits_strcasecmp(parName,"HISTORY") == 0) {
  776        ffphis( outfptr, result->value.data.str, status);
  777      } else if (fits_strcasecmp(parName,"COMMENT") == 0) {
  778        ffpcom( outfptr, result->value.data.str, status);
  779      } else {
  780        ffukys( outfptr, parName, result->value.data.str, parInfo, status );
  781      }
  782          break;
  783       }
  784    }
  785 
  786    ffcprs();
  787    FFUNLOCK;
  788    return( *status );
  789 }
  790 
  791 /*--------------------------------------------------------------------------*/
  792 int fftexp( fitsfile *fptr,      /* I - Input FITS file                     */
  793             char     *expr,      /* I - Arithmetic expression               */
  794             int      maxdim,     /* I - Max Dimension of naxes              */
  795             int      *datatype,  /* O - Data type of result                 */
  796             long     *nelem,     /* O - Vector length of result             */
  797             int      *naxis,     /* O - # of dimensions of result           */
  798             long     *naxes,     /* O - Size of each dimension              */
  799             int      *status )   /* O - Error status                        */
  800 /*                                                                          */
  801 /* Evaluate the given expression and return information on the result.      */
  802 /*--------------------------------------------------------------------------*/
  803 {
  804    FFLOCK;
  805    ffiprs( fptr, 0, expr, maxdim, datatype, nelem, naxis, naxes, status );
  806    ffcprs();
  807    FFUNLOCK;
  808    return( *status );
  809 }
  810 
  811 /*--------------------------------------------------------------------------*/
  812 int ffiprs( fitsfile *fptr,      /* I - Input FITS file                     */
  813             int      compressed, /* I - Is FITS file hkunexpanded?          */
  814             char     *expr,      /* I - Arithmetic expression               */
  815             int      maxdim,     /* I - Max Dimension of naxes              */
  816             int      *datatype,  /* O - Data type of result                 */
  817             long     *nelem,     /* O - Vector length of result             */
  818             int      *naxis,     /* O - # of dimensions of result           */
  819             long     *naxes,     /* O - Size of each dimension              */
  820             int      *status )   /* O - Error status                        */
  821 /*                                                                          */
  822 /* Initialize the parser and determine what type of result the expression   */
  823 /* produces.                                                                */
  824 /*--------------------------------------------------------------------------*/
  825 {
  826    Node *result;
  827    int  i,lexpr, tstatus = 0;
  828    int xaxis, bitpix;
  829    long xaxes[9];
  830    static iteratorCol dmyCol;
  831 
  832    if( *status ) return( *status );
  833 
  834    /* make sure all internal structures for this HDU are current */
  835    if ( ffrdef(fptr, status) ) return(*status);
  836 
  837    /*  Initialize the Parser structure  */
  838 
  839    gParse.def_fptr   = fptr;
  840    gParse.compressed = compressed;
  841    gParse.nCols      = 0;
  842    gParse.colData    = NULL;
  843    gParse.varData    = NULL;
  844    gParse.getData    = find_column;
  845    gParse.loadData   = load_column;
  846    gParse.Nodes      = NULL;
  847    gParse.nNodesAlloc= 0;
  848    gParse.nNodes     = 0;
  849    gParse.hdutype    = 0;
  850    gParse.status     = 0;
  851 
  852    fits_get_hdu_type(fptr, &gParse.hdutype, status );
  853 
  854    if (gParse.hdutype == IMAGE_HDU) {
  855 
  856       fits_get_img_param(fptr, 9, &bitpix, &xaxis, xaxes, status);
  857       if (*status) {
  858          ffpmsg("ffiprs: unable to get image dimensions");
  859          return( *status );
  860       }
  861       gParse.totalRows = xaxis > 0 ? 1 : 0;
  862       for (i = 0; i < xaxis; ++i)
  863          gParse.totalRows *= xaxes[i];
  864       if (DEBUG_PIXFILTER)
  865          printf("naxis=%d, gParse.totalRows=%ld\n", xaxis, gParse.totalRows);
  866    }
  867    else if( ffgkyj(fptr, "NAXIS2", &gParse.totalRows, 0, &tstatus) )
  868    {
  869       /* this might be a 1D or null image with no NAXIS2 keyword */
  870       gParse.totalRows = 0;
  871    } 
  872    
  873 
  874    /*  Copy expression into parser... read from file if necessary  */
  875 
  876 
  877    if( expr[0]=='@' ) {
  878       if( ffimport_file( expr+1, &gParse.expr, status ) ) return( *status );
  879       lexpr = strlen(gParse.expr);
  880    } else {
  881       lexpr = strlen(expr);
  882       gParse.expr = (char*)malloc( (2+lexpr)*sizeof(char));
  883       strcpy(gParse.expr,expr);
  884    }
  885    strcat(gParse.expr + lexpr,"\n");
  886    gParse.index    = 0;
  887    gParse.is_eobuf = 0;
  888 
  889    /*  Parse the expression, building the Nodes and determing  */
  890    /*  which columns are needed and what data type is returned  */
  891 
  892    ffrestart(NULL);
  893    if( ffparse() ) {
  894       return( *status = PARSE_SYNTAX_ERR );
  895    }
  896    /*  Check results  */
  897 
  898    *status = gParse.status;
  899    if( *status ) return(*status);
  900 
  901    if( !gParse.nNodes ) {
  902       ffpmsg("Blank expression");
  903       return( *status = PARSE_SYNTAX_ERR );
  904    }
  905    if( !gParse.nCols ) {
  906       dmyCol.fptr = fptr;         /* This allows iterator to know value of */
  907       gParse.colData = &dmyCol;   /* fptr when no columns are referenced   */
  908    }
  909 
  910    result = gParse.Nodes + gParse.resultNode;
  911 
  912    *naxis = result->value.naxis;
  913    *nelem = result->value.nelem;
  914    for( i=0; i<*naxis && i<maxdim; i++ )
  915       naxes[i] = result->value.naxes[i];
  916 
  917    switch( result->type ) {
  918    case BOOLEAN:
  919       *datatype = TLOGICAL;
  920       break;
  921    case LONG:
  922       *datatype = TLONG;
  923       break;
  924    case DOUBLE:
  925       *datatype = TDOUBLE;
  926       break;
  927    case BITSTR:
  928       *datatype = TBIT;
  929       break;
  930    case STRING:
  931       *datatype = TSTRING;
  932       break;
  933    default:
  934       *datatype = 0;
  935       ffpmsg("Bad return data type");
  936       *status = gParse.status = PARSE_BAD_TYPE;
  937       break;
  938    }
  939    gParse.datatype = *datatype;
  940    FREE(gParse.expr);
  941 
  942    if( result->operation==CONST_OP ) *nelem = - *nelem;
  943    return(*status);
  944 }
  945 
  946 /*--------------------------------------------------------------------------*/
  947 void ffcprs( void )  /*  No parameters                                      */
  948 /*                                                                          */
  949 /* Clear the parser, making it ready to accept a new expression.            */
  950 /*--------------------------------------------------------------------------*/
  951 {
  952    int col, node, i;
  953 
  954    if( gParse.nCols > 0 ) {
  955       FREE( gParse.colData  );
  956       for( col=0; col<gParse.nCols; col++ ) {
  957          if( gParse.varData[col].undef == NULL ) continue;
  958          if( gParse.varData[col].type  == BITSTR )
  959            FREE( ((char**)gParse.varData[col].data)[0] );
  960          free( gParse.varData[col].undef );
  961       }
  962       FREE( gParse.varData );
  963       gParse.nCols = 0;
  964    }
  965 
  966    if( gParse.nNodes > 0 ) {
  967       node = gParse.nNodes;
  968       while( node-- ) {
  969          if( gParse.Nodes[node].operation==gtifilt_fct ) {
  970             i = gParse.Nodes[node].SubNodes[0];
  971             if (gParse.Nodes[ i ].value.data.ptr)
  972             FREE( gParse.Nodes[ i ].value.data.ptr );
  973          }
  974          else if( gParse.Nodes[node].operation==regfilt_fct ) {
  975             i = gParse.Nodes[node].SubNodes[0];
  976             fits_free_region( (SAORegion *)gParse.Nodes[ i ].value.data.ptr );
  977          }
  978       }
  979       gParse.nNodes = 0;
  980    }
  981    if( gParse.Nodes ) free( gParse.Nodes );
  982    gParse.Nodes = NULL;
  983 
  984    gParse.hdutype = ANY_HDU;
  985    gParse.pixFilter = 0;
  986 }
  987 
  988 /*---------------------------------------------------------------------------*/
  989 int parse_data( long    totalrows,     /* I - Total rows to be processed     */
  990                 long    offset,        /* I - Number of rows skipped at start*/
  991                 long    firstrow,      /* I - First row of this iteration    */
  992                 long    nrows,         /* I - Number of rows in this iter    */
  993                 int      nCols,        /* I - Number of columns in use       */
  994                 iteratorCol *colData,  /* IO- Column information/data        */
  995                 void    *userPtr )     /* I - Data handling instructions     */
  996 /*                                                                           */
  997 /* Iterator work function which calls the parser and copies the results      */
  998 /* into either an OutputCol or a data pointer supplied in the userPtr        */
  999 /* structure.                                                                */
 1000 /*---------------------------------------------------------------------------*/
 1001 {
 1002     int status, constant=0, anyNullThisTime=0;
 1003     long jj, kk, idx, remain, ntodo;
 1004     Node *result;
 1005     iteratorCol * outcol;
 1006 
 1007     /* declare variables static to preserve their values between calls */
 1008     static void *Data, *Null;
 1009     static int  datasize;
 1010     static long lastRow, repeat, resDataSize;
 1011     static LONGLONG jnull;
 1012     static parseInfo *userInfo;
 1013     static long zeros[4] = {0,0,0,0};
 1014 
 1015     if (DEBUG_PIXFILTER)
 1016        printf("parse_data(total=%ld, offset=%ld, first=%ld, rows=%ld, cols=%d)\n",
 1017                 totalrows, offset, firstrow, nrows, nCols);
 1018     /*--------------------------------------------------------*/
 1019     /*  Initialization procedures: execute on the first call  */
 1020     /*--------------------------------------------------------*/
 1021     outcol = colData + (nCols - 1);
 1022     if (firstrow == offset+1)
 1023     {
 1024        userInfo = (parseInfo*)userPtr;
 1025        userInfo->anyNull = 0;
 1026 
 1027        if( userInfo->maxRows>0 )
 1028           userInfo->maxRows = minvalue(totalrows,userInfo->maxRows);
 1029        else if( userInfo->maxRows<0 )
 1030           userInfo->maxRows = totalrows;
 1031        else
 1032           userInfo->maxRows = nrows;
 1033 
 1034        lastRow = firstrow + userInfo->maxRows - 1;
 1035 
 1036        if( userInfo->dataPtr==NULL ) {
 1037 
 1038           if( outcol->iotype == InputCol ) {
 1039              ffpmsg("Output column for parser results not found!");
 1040              return( PARSE_NO_OUTPUT );
 1041           }
 1042           /* Data gets set later */
 1043           Null = outcol->array;
 1044           userInfo->datatype = outcol->datatype;
 1045 
 1046           /* Check for a TNULL/BLANK keyword for output column/image */
 1047 
 1048           status = 0;
 1049           jnull = 0;
 1050           if (gParse.hdutype == IMAGE_HDU) {
 1051              if (gParse.pixFilter->blank)
 1052                 jnull = (LONGLONG) gParse.pixFilter->blank;
 1053           }
 1054           else {
 1055              ffgknjj( outcol->fptr, "TNULL", outcol->colnum,
 1056                         1, &jnull, (int*)&jj, &status );
 1057 
 1058              if( status==BAD_INTKEY ) {
 1059                 /*  Probably ASCII table with text TNULL keyword  */
 1060                 switch( userInfo->datatype ) {
 1061                    case TSHORT:  jnull = (LONGLONG) SHRT_MIN;      break;
 1062                    case TINT:    jnull = (LONGLONG) INT_MIN;       break;
 1063                    case TLONG:   jnull = (LONGLONG) LONG_MIN;      break;
 1064                 }
 1065              }
 1066           }
 1067           repeat = outcol->repeat;
 1068 /*
 1069           if (DEBUG_PIXFILTER)
 1070             printf("parse_data: using null value %ld\n", jnull);
 1071 */
 1072        } else {
 1073 
 1074           Data = userInfo->dataPtr;
 1075           Null = (userInfo->nullPtr ? userInfo->nullPtr : zeros);
 1076           repeat = gParse.Nodes[gParse.resultNode].value.nelem;
 1077 
 1078        }
 1079 
 1080        /* Determine the size of each element of the returned result */
 1081 
 1082        switch( userInfo->datatype ) {
 1083        case TBIT:       /*  Fall through to TBYTE  */
 1084        case TLOGICAL:   /*  Fall through to TBYTE  */
 1085        case TBYTE:     datasize = sizeof(char);     break;
 1086        case TSHORT:    datasize = sizeof(short);    break;
 1087        case TINT:      datasize = sizeof(int);      break;
 1088        case TLONG:     datasize = sizeof(long);     break;
 1089        case TLONGLONG: datasize = sizeof(LONGLONG); break;
 1090        case TFLOAT:    datasize = sizeof(float);    break;
 1091        case TDOUBLE:   datasize = sizeof(double);   break;
 1092        case TSTRING:   datasize = sizeof(char*);    break;
 1093        }
 1094 
 1095        /* Determine the size of each element of the calculated result */
 1096        /*   (only matters for numeric/logical data)                   */
 1097 
 1098        switch( gParse.Nodes[gParse.resultNode].type ) {
 1099        case BOOLEAN:   resDataSize = sizeof(char);    break;
 1100        case LONG:      resDataSize = sizeof(long);    break;
 1101        case DOUBLE:    resDataSize = sizeof(double);  break;
 1102        }
 1103     }
 1104 
 1105     /*-------------------------------------------*/
 1106     /*  Main loop: process all the rows of data  */
 1107     /*-------------------------------------------*/
 1108 
 1109     /*  If writing to output column, set first element to appropriate  */
 1110     /*  null value.  If no NULLs encounter, zero out before returning. */
 1111 /*
 1112           if (DEBUG_PIXFILTER)
 1113             printf("parse_data: using null value %ld\n", jnull);
 1114 */
 1115 
 1116     if( userInfo->dataPtr == NULL ) {
 1117        /* First, reset Data pointer to start of output array */
 1118        Data = (char*) outcol->array + datasize;
 1119 
 1120        switch( userInfo->datatype ) {
 1121        case TLOGICAL: *(char  *)Null = 'U';             break;
 1122        case TBYTE:    *(char  *)Null = (char )jnull;    break;
 1123        case TSHORT:   *(short *)Null = (short)jnull;    break;
 1124        case TINT:     *(int   *)Null = (int  )jnull;    break;
 1125        case TLONG:    *(long  *)Null = (long )jnull;    break;
 1126        case TLONGLONG: *(LONGLONG  *)Null = (LONGLONG )jnull;    break;
 1127        case TFLOAT:   *(float *)Null = FLOATNULLVALUE;  break;
 1128        case TDOUBLE:  *(double*)Null = DOUBLENULLVALUE; break;
 1129        case TSTRING: (*(char **)Null)[0] = '\1';
 1130                      (*(char **)Null)[1] = '\0';        break;
 1131        }
 1132     }
 1133 
 1134     /* Alter nrows in case calling routine didn't want to do all rows */
 1135 
 1136     nrows = minvalue(nrows,lastRow-firstrow+1);
 1137 
 1138     Setup_DataArrays( nCols, colData, firstrow, nrows );
 1139 
 1140     /* Parser allocates arrays for each column and calculation it performs. */
 1141     /* Limit number of rows processed during each pass to reduce memory     */
 1142     /* requirements... In most cases, iterator will limit rows to less      */
 1143     /* than 2500 rows per iteration, so this is really only relevant for    */
 1144     /* hk-compressed files which must be decompressed in memory and sent    */
 1145     /* whole to parse_data in a single iteration.                           */
 1146 
 1147     remain = nrows;
 1148     while( remain ) {
 1149        ntodo = minvalue(remain,2500);
 1150        Evaluate_Parser ( firstrow, ntodo );
 1151        if( gParse.status ) break;
 1152 
 1153        firstrow += ntodo;
 1154        remain   -= ntodo;
 1155 
 1156        /*  Copy results into data array  */
 1157 
 1158        result = gParse.Nodes + gParse.resultNode;
 1159        if( result->operation==CONST_OP ) constant = 1;
 1160 
 1161        switch( result->type ) {
 1162 
 1163        case BOOLEAN:
 1164        case LONG:
 1165        case DOUBLE:
 1166           if( constant ) {
 1167              char undef=0;
 1168              for( kk=0; kk<ntodo; kk++ )
 1169                 for( jj=0; jj<repeat; jj++ )
 1170                    ffcvtn( gParse.datatype,
 1171                            &(result->value.data),
 1172                            &undef, result->value.nelem /* 1 */,
 1173                            userInfo->datatype, Null,
 1174                            (char*)Data + (kk*repeat+jj)*datasize,
 1175                            &anyNullThisTime, &gParse.status );
 1176           } else {
 1177              if ( repeat == result->value.nelem ) {
 1178                 ffcvtn( gParse.datatype,
 1179                         result->value.data.ptr,
 1180                         result->value.undef,
 1181                         result->value.nelem*ntodo,
 1182                         userInfo->datatype, Null, Data,
 1183                         &anyNullThisTime, &gParse.status );
 1184              } else if( result->value.nelem == 1 ) {
 1185                 for( kk=0; kk<ntodo; kk++ )
 1186                    for( jj=0; jj<repeat; jj++ ) {
 1187                       ffcvtn( gParse.datatype,
 1188                               (char*)result->value.data.ptr + kk*resDataSize,
 1189                               (char*)result->value.undef + kk,
 1190                               1, userInfo->datatype, Null,
 1191                               (char*)Data + (kk*repeat+jj)*datasize,
 1192                               &anyNullThisTime, &gParse.status );
 1193                    }
 1194              } else {
 1195                 int nCopy;
 1196                 nCopy = minvalue( repeat, result->value.nelem );
 1197                 for( kk=0; kk<ntodo; kk++ ) {
 1198                    ffcvtn( gParse.datatype,
 1199                            (char*)result->value.data.ptr
 1200                                   + kk*result->value.nelem*resDataSize,
 1201                            (char*)result->value.undef
 1202                                   + kk*result->value.nelem,
 1203                            nCopy, userInfo->datatype, Null,
 1204                            (char*)Data + (kk*repeat)*datasize,
 1205                            &anyNullThisTime, &gParse.status );
 1206                    if( nCopy < repeat ) {
 1207                       memset( (char*)Data + (kk*repeat+nCopy)*datasize,
 1208                               0, (repeat-nCopy)*datasize);
 1209                    }
 1210                 }
 1211 
 1212              }
 1213              if( result->operation>0 ) {
 1214                 FREE( result->value.data.ptr );
 1215              }
 1216           }
 1217           if( gParse.status==OVERFLOW_ERR ) {
 1218              gParse.status = NUM_OVERFLOW;
 1219              ffpmsg("Numerical overflow while converting expression to necessary datatype");
 1220           }
 1221           break;
 1222 
 1223        case BITSTR:
 1224           switch( userInfo->datatype ) {
 1225           case TBYTE:
 1226              idx = -1;
 1227              for( kk=0; kk<ntodo; kk++ ) {
 1228                 for( jj=0; jj<result->value.nelem; jj++ ) {
 1229                    if( jj%8 == 0 )
 1230                       ((char*)Data)[++idx] = 0;
 1231                    if( constant ) {
 1232                       if( result->value.data.str[jj]=='1' )
 1233                          ((char*)Data)[idx] |= 128>>(jj%8);
 1234                    } else {
 1235                       if( result->value.data.strptr[kk][jj]=='1' )
 1236                          ((char*)Data)[idx] |= 128>>(jj%8);
 1237                    }
 1238                 }
 1239              }
 1240              break;
 1241           case TBIT:
 1242           case TLOGICAL:
 1243              if( constant ) {
 1244                 for( kk=0; kk<ntodo; kk++ )
 1245                    for( jj=0; jj<result->value.nelem; jj++ ) {
 1246                       ((char*)Data)[ jj+kk*result->value.nelem ] =
 1247                          ( result->value.data.str[jj]=='1' );
 1248                    }
 1249              } else {
 1250                 for( kk=0; kk<ntodo; kk++ )
 1251                    for( jj=0; jj<result->value.nelem; jj++ ) {
 1252                       ((char*)Data)[ jj+kk*result->value.nelem ] =
 1253                          ( result->value.data.strptr[kk][jj]=='1' );
 1254                    }
 1255              }
 1256              break; 
 1257           case TSTRING:
 1258              if( constant ) {
 1259                 for( jj=0; jj<ntodo; jj++ ) {
 1260                    strcpy( ((char**)Data)[jj], result->value.data.str );
 1261                 }
 1262              } else {
 1263                 for( jj=0; jj<ntodo; jj++ ) {
 1264                    strcpy( ((char**)Data)[jj], result->value.data.strptr[jj] );
 1265                 }
 1266              }
 1267              break;
 1268           default:
 1269              ffpmsg("Cannot convert bit expression to desired type.");
 1270              gParse.status = PARSE_BAD_TYPE;
 1271              break;
 1272           }
 1273           if( result->operation>0 ) {
 1274              FREE( result->value.data.strptr[0] );
 1275              FREE( result->value.data.strptr );
 1276           }
 1277           break;
 1278 
 1279        case STRING:
 1280           if( userInfo->datatype==TSTRING ) {
 1281              if( constant ) {
 1282                 for( jj=0; jj<ntodo; jj++ )
 1283                    strcpy( ((char**)Data)[jj], result->value.data.str );
 1284              } else {
 1285                 for( jj=0; jj<ntodo; jj++ )
 1286                    if( result->value.undef[jj] ) {
 1287                       anyNullThisTime = 1;
 1288                       strcpy( ((char**)Data)[jj],
 1289                               *(char **)Null );
 1290                    } else {
 1291                       strcpy( ((char**)Data)[jj],
 1292                               result->value.data.strptr[jj] );
 1293                    }
 1294              }
 1295           } else {
 1296              ffpmsg("Cannot convert string expression to desired type.");
 1297              gParse.status = PARSE_BAD_TYPE;
 1298           }
 1299           if( result->operation>0 ) {
 1300              FREE( result->value.data.strptr[0] );
 1301              FREE( result->value.data.strptr );
 1302           }
 1303           break;
 1304        }
 1305 
 1306        if( gParse.status ) break;
 1307 
 1308        /*  Increment Data to point to where the next block should go  */
 1309 
 1310        if( result->type==BITSTR && userInfo->datatype==TBYTE )
 1311           Data = (char*)Data
 1312                     + datasize * ( (result->value.nelem+7)/8 ) * ntodo;
 1313        else if( result->type==STRING )
 1314           Data = (char*)Data + datasize * ntodo;
 1315        else
 1316           Data = (char*)Data + datasize * ntodo * repeat;
 1317     }
 1318 
 1319     /* If no NULLs encountered during this pass, set Null value to */
 1320     /* zero to make the writing of the output column data faster   */
 1321 
 1322     if( anyNullThisTime )
 1323        userInfo->anyNull = 1;
 1324     else if( userInfo->dataPtr == NULL ) {
 1325        if( userInfo->datatype == TSTRING )
 1326           memcpy( *(char **)Null, zeros, 2 );
 1327        else 
 1328           memcpy( Null, zeros, datasize );
 1329     }
 1330 
 1331     /*-------------------------------------------------------*/
 1332     /*  Clean up procedures:  after processing all the rows  */
 1333     /*-------------------------------------------------------*/
 1334 
 1335     /*  if the calling routine specified that only a limited number    */
 1336     /*  of rows in the table should be processed, return a value of -1 */
 1337     /*  once all the rows have been done, if no other error occurred.  */
 1338 
 1339     if (gParse.hdutype != IMAGE_HDU && firstrow - 1 == lastRow) {
 1340            if (!gParse.status && userInfo->maxRows<totalrows) {
 1341                   return (-1);
 1342            }
 1343     }
 1344 
 1345     return(gParse.status);  /* return successful status */
 1346 }
 1347 
 1348 static void Setup_DataArrays( int nCols, iteratorCol *cols,
 1349                               long fRow, long nRows )
 1350     /***********************************************************************/
 1351     /*  Setup the varData array in gParse to contain the fits column data. */
 1352     /*  Then, allocate and initialize the necessary UNDEF arrays for each  */
 1353     /*  column used by the parser.                                         */
 1354     /***********************************************************************/
 1355 {
 1356    int     i;
 1357    long    nelem, len, row, idx;
 1358    char  **bitStrs;
 1359    char  **sptr;
 1360    char   *barray;
 1361    long   *iarray;
 1362    double *rarray;
 1363    char msg[80];
 1364 
 1365    gParse.firstDataRow = fRow;
 1366    gParse.nDataRows    = nRows;
 1367 
 1368    /*  Resize and fill in UNDEF arrays for each column  */
 1369 
 1370    for( i=0; i<nCols; i++ ) {
 1371 
 1372       iteratorCol *icol = cols + i;
 1373       DataInfo *varData = gParse.varData + i;
 1374 
 1375       if( icol->iotype == OutputCol ) continue;
 1376 
 1377       nelem  = varData->nelem;
 1378       len    = nelem * nRows;
 1379 
 1380       switch ( varData->type ) {
 1381 
 1382       case BITSTR:
 1383       /* No need for UNDEF array, but must make string DATA array */
 1384          len = (nelem+1)*nRows;   /* Count '\0' */
 1385          bitStrs = (char**)varData->data;
 1386          if( bitStrs ) FREE( bitStrs[0] );
 1387          free( bitStrs );
 1388          bitStrs = (char**)malloc( nRows*sizeof(char*) );
 1389          if( bitStrs==NULL ) {
 1390             varData->data = varData->undef = NULL;
 1391             gParse.status = MEMORY_ALLOCATION;
 1392             break;
 1393          }
 1394          bitStrs[0] = (char*)malloc( len*sizeof(char) );
 1395          if( bitStrs[0]==NULL ) {
 1396             free( bitStrs );
 1397             varData->data = varData->undef = NULL;
 1398             gParse.status = MEMORY_ALLOCATION;
 1399             break;
 1400          }
 1401 
 1402          for( row=0; row<nRows; row++ ) {
 1403             bitStrs[row] = bitStrs[0] + row*(nelem+1);
 1404             idx = (row)*( (nelem+7)/8 ) + 1;
 1405             for(len=0; len<nelem; len++) {
 1406                if( ((char*)icol->array)[idx] & (1<<(7-len%8)) )
 1407                   bitStrs[row][len] = '1';
 1408                else
 1409                   bitStrs[row][len] = '0';
 1410                if( len%8==7 ) idx++;
 1411             }
 1412             bitStrs[row][len] = '\0';
 1413          }
 1414          varData->undef = (char*)bitStrs;
 1415          varData->data  = (char*)bitStrs;
 1416          break;
 1417 
 1418       case STRING:
 1419          sptr = (char**)icol->array;
 1420          if (varData->undef)
 1421             free( varData->undef );
 1422          varData->undef = (char*)malloc( nRows*sizeof(char) );
 1423          if( varData->undef==NULL ) {
 1424             gParse.status = MEMORY_ALLOCATION;
 1425             break;
 1426          }
 1427          row = nRows;
 1428          while( row-- )
 1429             varData->undef[row] =
 1430                ( **sptr != '\0' && FSTRCMP( sptr[0], sptr[row+1] )==0 );
 1431          varData->data  = sptr + 1;
 1432          break;
 1433 
 1434       case BOOLEAN:
 1435          barray = (char*)icol->array;
 1436          if (varData->undef)
 1437             free( varData->undef );
 1438          varData->undef = (char*)malloc( len*sizeof(char) );
 1439          if( varData->undef==NULL ) {
 1440             gParse.status = MEMORY_ALLOCATION;
 1441             break;
 1442          }
 1443          while( len-- ) {
 1444             varData->undef[len] = 
 1445                ( barray[0]!=0 && barray[0]==barray[len+1] );
 1446          }
 1447          varData->data  = barray + 1;
 1448          break;
 1449 
 1450       case LONG:
 1451          iarray = (long*)icol->array;
 1452          if (varData->undef)
 1453             free( varData->undef );
 1454          varData->undef = (char*)malloc( len*sizeof(char) );
 1455          if( varData->undef==NULL ) {
 1456             gParse.status = MEMORY_ALLOCATION;
 1457             break;
 1458          }
 1459          while( len-- ) {
 1460             varData->undef[len] = 
 1461                ( iarray[0]!=0L && iarray[0]==iarray[len+1] );
 1462          }
 1463          varData->data  = iarray + 1;
 1464          break;
 1465 
 1466       case DOUBLE:
 1467          rarray = (double*)icol->array;
 1468          if (varData->undef)
 1469             free( varData->undef );
 1470          varData->undef = (char*)malloc( len*sizeof(char) );
 1471          if( varData->undef==NULL ) {
 1472             gParse.status = MEMORY_ALLOCATION;
 1473             break;
 1474          }
 1475          while( len-- ) {
 1476             varData->undef[len] = 
 1477                ( rarray[0]!=0.0 && rarray[0]==rarray[len+1]);
 1478          }
 1479          varData->data  = rarray + 1;
 1480          break;
 1481 
 1482       default:
 1483          snprintf(msg, 80, "SetupDataArrays, unhandled type %d\n",
 1484                 varData->type);
 1485          ffpmsg(msg);
 1486       }
 1487 
 1488       if( gParse.status ) {  /*  Deallocate NULL arrays of previous columns */
 1489          while( i-- ) {
 1490             varData = gParse.varData + i;
 1491             if( varData->type==BITSTR )
 1492                FREE( ((char**)varData->data)[0] );
 1493             FREE( varData->undef );
 1494             varData->undef = NULL;
 1495          }
 1496          return;
 1497       }
 1498    }
 1499 }
 1500 
 1501 /*--------------------------------------------------------------------------*/
 1502 int ffcvtn( int   inputType,  /* I - Data type of input array               */
 1503             void  *input,     /* I - Input array of type inputType          */
 1504             char  *undef,     /* I - Array of flags indicating UNDEF elems  */
 1505             long  ntodo,      /* I - Number of elements to process          */
 1506             int   outputType, /* I - Data type of output array              */
 1507             void  *nulval,    /* I - Ptr to value to use for UNDEF elements */
 1508             void  *output,    /* O - Output array of type outputType        */
 1509             int   *anynull,   /* O - Any nulls flagged?                     */
 1510             int   *status )   /* O - Error status                           */
 1511 /*                                                                          */
 1512 /* Convert an array of any input data type to an array of any output        */
 1513 /* data type, using an array of UNDEF flags to assign nulvals to            */
 1514 /*--------------------------------------------------------------------------*/
 1515 {
 1516    long i;
 1517 
 1518    switch( outputType ) {
 1519 
 1520    case TLOGICAL:
 1521       switch( inputType ) {
 1522       case TLOGICAL:
 1523       case TBYTE:
 1524          for( i=0; i<ntodo; i++ )
 1525             if( ((unsigned char*)input)[i] )
 1526                 ((unsigned char*)output)[i] = 1;
 1527             else
 1528                 ((unsigned char*)output)[i] = 0;
 1529          break;
 1530       case TSHORT:
 1531          for( i=0; i<ntodo; i++ )
 1532             if( ((short*)input)[i] )
 1533                 ((unsigned char*)output)[i] = 1;
 1534             else
 1535                 ((unsigned char*)output)[i] = 0;
 1536          break;
 1537       case TLONG:
 1538          for( i=0; i<ntodo; i++ )
 1539             if( ((long*)input)[i] )
 1540                 ((unsigned char*)output)[i] = 1;
 1541             else
 1542                 ((unsigned char*)output)[i] = 0;
 1543          break;
 1544       case TFLOAT:
 1545          for( i=0; i<ntodo; i++ )
 1546             if( ((float*)input)[i] )
 1547                 ((unsigned char*)output)[i] = 1;
 1548             else
 1549                 ((unsigned char*)output)[i] = 0;
 1550          break;
 1551       case TDOUBLE:
 1552          for( i=0; i<ntodo; i++ )
 1553             if( ((double*)input)[i] )
 1554                 ((unsigned char*)output)[i] = 1;
 1555             else
 1556                 ((unsigned char*)output)[i] = 0;
 1557          break;
 1558       default:
 1559          *status = BAD_DATATYPE;
 1560          break;
 1561       }
 1562       for(i=0;i<ntodo;i++) {
 1563          if( undef[i] ) {
 1564             ((unsigned char*)output)[i] = *(unsigned char*)nulval;
 1565             *anynull = 1;
 1566          }
 1567       }
 1568       break;
 1569 
 1570    case TBYTE:
 1571       switch( inputType ) {
 1572       case TLOGICAL:
 1573       case TBYTE:
 1574          for( i=0; i<ntodo; i++ )
 1575             ((unsigned char*)output)[i] = ((unsigned char*)input)[i];
 1576          break;
 1577       case TSHORT:
 1578          fffi2i1((short*)input,ntodo,1.,0.,0,0,0,NULL,NULL,(unsigned char*)output,status);
 1579          break;
 1580       case TLONG:
 1581          for (i = 0; i < ntodo; i++) {
 1582             if( undef[i] ) {
 1583                ((unsigned char*)output)[i] = *(unsigned char*)nulval;
 1584                *anynull = 1;
 1585             } else {
 1586                if( ((long*)input)[i] < 0 ) {
 1587                   *status = OVERFLOW_ERR;
 1588                   ((unsigned char*)output)[i] = 0;
 1589                } else if( ((long*)input)[i] > UCHAR_MAX ) {
 1590                   *status = OVERFLOW_ERR;
 1591                   ((unsigned char*)output)[i] = UCHAR_MAX;
 1592                } else
 1593                   ((unsigned char*)output)[i] = 
 1594                      (unsigned char) ((long*)input)[i];
 1595             }
 1596          }
 1597          return( *status );
 1598       case TFLOAT:
 1599          fffr4i1((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
 1600                  (unsigned char*)output,status);
 1601          break;
 1602       case TDOUBLE:
 1603          fffr8i1((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
 1604                  (unsigned char*)output,status);
 1605          break;
 1606       default:
 1607          *status = BAD_DATATYPE;
 1608          break;
 1609       }
 1610       for(i=0;i<ntodo;i++) {
 1611          if( undef[i] ) {
 1612             ((unsigned char*)output)[i] = *(unsigned char*)nulval;
 1613             *anynull = 1;
 1614          }
 1615       }
 1616       break;
 1617 
 1618    case TSHORT:
 1619       switch( inputType ) {
 1620       case TLOGICAL:
 1621       case TBYTE:
 1622          for( i=0; i<ntodo; i++ )
 1623             ((short*)output)[i] = ((unsigned char*)input)[i];
 1624          break;
 1625       case TSHORT:
 1626          for( i=0; i<ntodo; i++ )
 1627             ((short*)output)[i] = ((short*)input)[i];
 1628          break;
 1629       case TLONG:
 1630          for (i = 0; i < ntodo; i++) {
 1631             if( undef[i] ) {
 1632                ((short*)output)[i] = *(short*)nulval;
 1633                *anynull = 1;
 1634             } else {
 1635                if( ((long*)input)[i] < SHRT_MIN ) {
 1636                   *status = OVERFLOW_ERR;
 1637                   ((short*)output)[i] = SHRT_MIN;
 1638                } else if ( ((long*)input)[i] > SHRT_MAX ) {
 1639                   *status = OVERFLOW_ERR;
 1640                   ((short*)output)[i] = SHRT_MAX;
 1641                } else
 1642                   ((short*)output)[i] = (short) ((long*)input)[i];
 1643             }
 1644          }
 1645          return( *status );
 1646       case TFLOAT:
 1647          fffr4i2((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
 1648                  (short*)output,status);
 1649          break;
 1650       case TDOUBLE:
 1651          fffr8i2((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
 1652                  (short*)output,status);
 1653          break;
 1654       default:
 1655          *status = BAD_DATATYPE;
 1656          break;
 1657       }
 1658       for(i=0;i<ntodo;i++) {
 1659          if( undef[i] ) {
 1660             ((short*)output)[i] = *(short*)nulval;
 1661             *anynull = 1;
 1662          }
 1663       }
 1664       break;
 1665 
 1666    case TINT:
 1667       switch( inputType ) {
 1668       case TLOGICAL:
 1669       case TBYTE:
 1670          for( i=0; i<ntodo; i++ )
 1671             ((int*)output)[i] = ((unsigned char*)input)[i];
 1672          break;
 1673       case TSHORT:
 1674          for( i=0; i<ntodo; i++ )
 1675             ((int*)output)[i] = ((short*)input)[i];
 1676          break;
 1677       case TLONG:
 1678          for( i=0; i<ntodo; i++ )
 1679             ((int*)output)[i] = ((long*)input)[i];
 1680          break;
 1681       case TFLOAT:
 1682          fffr4int((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
 1683                   (int*)output,status);
 1684          break;
 1685       case TDOUBLE:
 1686          fffr8int((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
 1687                   (int*)output,status);
 1688          break;
 1689       default:
 1690          *status = BAD_DATATYPE;
 1691          break;
 1692       }
 1693       for(i=0;i<ntodo;i++) {
 1694          if( undef[i] ) {
 1695             ((int*)output)[i] = *(int*)nulval;
 1696             *anynull = 1;
 1697          }
 1698       }
 1699       break;
 1700 
 1701    case TLONG:
 1702       switch( inputType ) {
 1703       case TLOGICAL:
 1704       case TBYTE:
 1705          for( i=0; i<ntodo; i++ )
 1706             ((long*)output)[i] = ((unsigned char*)input)[i];
 1707          break;
 1708       case TSHORT:
 1709          for( i=0; i<ntodo; i++ )
 1710             ((long*)output)[i] = ((short*)input)[i];
 1711          break;
 1712       case TLONG:
 1713          for( i=0; i<ntodo; i++ )
 1714             ((long*)output)[i] = ((long*)input)[i];
 1715          break;
 1716       case TFLOAT:
 1717          fffr4i4((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
 1718                  (long*)output,status);
 1719          break;
 1720       case TDOUBLE:
 1721          fffr8i4((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
 1722                  (long*)output,status);
 1723          break;
 1724       default:
 1725          *status = BAD_DATATYPE;
 1726          break;
 1727       }
 1728       for(i=0;i<ntodo;i++) {
 1729          if( undef[i] ) {
 1730             ((long*)output)[i] = *(long*)nulval;
 1731             *anynull = 1;
 1732          }
 1733       }
 1734       break;
 1735 
 1736    case TLONGLONG:
 1737       switch( inputType ) {
 1738       case TLOGICAL:
 1739       case TBYTE:
 1740          for( i=0; i<ntodo; i++ )
 1741             ((LONGLONG*)output)[i] = ((unsigned char*)input)[i];
 1742          break;
 1743       case TSHORT:
 1744          for( i=0; i<ntodo; i++ )
 1745             ((LONGLONG*)output)[i] = ((short*)input)[i];
 1746          break;
 1747       case TLONG:
 1748          for( i=0; i<ntodo; i++ )
 1749             ((LONGLONG*)output)[i] = ((long*)input)[i];
 1750          break;
 1751       case TFLOAT:
 1752          fffr4i8((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
 1753                  (LONGLONG*)output,status);
 1754          break;
 1755       case TDOUBLE:
 1756          fffr8i8((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
 1757                  (LONGLONG*)output,status);
 1758 
 1759          break;
 1760       default:
 1761          *status = BAD_DATATYPE;
 1762          break;
 1763       }
 1764       for(i=0;i<ntodo;i++) {
 1765          if( undef[i] ) {
 1766             ((LONGLONG*)output)[i] = *(LONGLONG*)nulval;
 1767             *anynull = 1;
 1768          }
 1769       }
 1770       break;
 1771 
 1772    case TFLOAT:
 1773       switch( inputType ) {
 1774       case TLOGICAL:
 1775       case TBYTE:
 1776          for( i=0; i<ntodo; i++ )
 1777             ((float*)output)[i] = ((unsigned char*)input)[i];
 1778          break;
 1779       case TSHORT:
 1780          for( i=0; i<ntodo; i++ )
 1781             ((float*)output)[i] = ((short*)input)[i];
 1782          break;
 1783       case TLONG:
 1784          for( i=0; i<ntodo; i++ )
 1785             ((float*)output)[i] = (float) ((long*)input)[i];
 1786          break;
 1787       case TFLOAT:
 1788          for( i=0; i<ntodo; i++ )
 1789             ((float*)output)[i] = ((float*)input)[i];
 1790          break;
 1791       case TDOUBLE:
 1792          fffr8r4((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
 1793                  (float*)output,status);
 1794          break;
 1795       default:
 1796          *status = BAD_DATATYPE;
 1797          break;
 1798       }
 1799       for(i=0;i<ntodo;i++) {
 1800          if( undef[i] ) {
 1801             ((float*)output)[i] = *(float*)nulval;
 1802             *anynull = 1;
 1803          }
 1804       }
 1805       break;
 1806 
 1807    case TDOUBLE:
 1808       switch( inputType ) {
 1809       case TLOGICAL:
 1810       case TBYTE:
 1811          for( i=0; i<ntodo; i++ )
 1812             ((double*)output)[i] = ((unsigned char*)input)[i];
 1813          break;
 1814       case TSHORT:
 1815          for( i=0; i<ntodo; i++ )
 1816             ((double*)output)[i] = ((short*)input)[i];
 1817          break;
 1818       case TLONG:
 1819          for( i=0; i<ntodo; i++ )
 1820             ((double*)output)[i] = ((long*)input)[i];
 1821          break;
 1822       case TFLOAT:
 1823          for( i=0; i<ntodo; i++ )
 1824             ((double*)output)[i] = ((float*)input)[i];
 1825          break;
 1826       case TDOUBLE:
 1827          for( i=0; i<ntodo; i++ )
 1828             ((double*)output)[i] = ((double*)input)[i];
 1829          break;
 1830       default:
 1831          *status = BAD_DATATYPE;
 1832          break;
 1833       }
 1834       for(i=0;i<ntodo;i++) {
 1835          if( undef[i] ) {
 1836             ((double*)output)[i] = *(double*)nulval;
 1837             *anynull = 1;
 1838          }
 1839       }
 1840       break;
 1841 
 1842    default:
 1843       *status = BAD_DATATYPE;
 1844       break;
 1845    }
 1846 
 1847    return ( *status );
 1848 }
 1849 
 1850 /*---------------------------------------------------------------------------*/
 1851 int fffrwc( fitsfile *fptr,        /* I - Input FITS file                    */
 1852             char     *expr,        /* I - Boolean expression                 */
 1853             char     *timeCol,     /* I - Name of time column                */
 1854             char     *parCol,      /* I - Name of parameter column           */
 1855             char     *valCol,      /* I - Name of value column               */
 1856             long     ntimes,       /* I - Number of distinct times in file   */
 1857             double   *times,       /* O - Array of times in file             */
 1858             char     *time_status, /* O - Array of boolean results           */
 1859             int      *status )     /* O - Error status                       */
 1860 /*                                                                           */
 1861 /* Evaluate a boolean expression for each time in a compressed file,         */
 1862 /* returning an array of flags indicating which times evaluated to TRUE/FALSE*/
 1863 /*---------------------------------------------------------------------------*/
 1864 {
 1865    parseInfo Info;
 1866    long alen, width;
 1867    int parNo, typecode;
 1868    int naxis, constant, nCol=0;
 1869    long nelem, naxes[MAXDIMS], elem;
 1870    char result;
 1871 
 1872    if( *status ) return( *status );
 1873 
 1874    fits_get_colnum( fptr, CASEINSEN, timeCol, &gParse.timeCol, status );
 1875    fits_get_colnum( fptr, CASEINSEN, parCol,  &gParse.parCol , status );
 1876    fits_get_colnum( fptr, CASEINSEN, valCol,  &gParse.valCol, status );
 1877    if( *status ) return( *status );
 1878    
 1879    if( ffiprs( fptr, 1, expr, MAXDIMS, &Info.datatype, &nelem,
 1880                &naxis, naxes, status ) ) {
 1881       ffcprs();
 1882       return( *status );
 1883    }
 1884    if( nelem<0 ) {
 1885       constant = 1;
 1886       nelem = -nelem;
 1887       nCol = gParse.nCols;
 1888       gParse.nCols = 0;    /*  Ignore all column references  */
 1889    } else
 1890       constant = 0;
 1891 
 1892    if( Info.datatype!=TLOGICAL || nelem!=1 ) {
 1893       ffcprs();
 1894       ffpmsg("Expression does not evaluate to a logical scalar.");
 1895       return( *status = PARSE_BAD_TYPE );
 1896    }
 1897 
 1898    /*******************************************/
 1899    /* Allocate data arrays for each parameter */
 1900    /*******************************************/
 1901    
 1902    parNo = gParse.nCols;
 1903    while( parNo-- ) {
 1904       switch( gParse.colData[parNo].datatype ) {
 1905       case TLONG:
 1906          if( (gParse.colData[parNo].array =
 1907               (long *)malloc( (ntimes+1)*sizeof(long) )) )
 1908             ((long*)gParse.colData[parNo].array)[0] = 1234554321;
 1909          else
 1910             *status = MEMORY_ALLOCATION;
 1911          break;
 1912       case TDOUBLE:
 1913          if( (gParse.colData[parNo].array =
 1914               (double *)malloc( (ntimes+1)*sizeof(double) )) )
 1915             ((double*)gParse.colData[parNo].array)[0] = DOUBLENULLVALUE;
 1916          else
 1917             *status = MEMORY_ALLOCATION;
 1918          break;
 1919       case TSTRING:
 1920          if( !fits_get_coltype( fptr, gParse.valCol, &typecode,
 1921                                 &alen, &width, status ) ) {
 1922             alen++;
 1923             if( (gParse.colData[parNo].array =
 1924                  (char **)malloc( (ntimes+1)*sizeof(char*) )) ) {
 1925                if( (((char **)gParse.colData[parNo].array)[0] =
 1926                     (char *)malloc( (ntimes+1)*sizeof(char)*alen )) ) {
 1927                   for( elem=1; elem<=ntimes; elem++ )
 1928                      ((char **)gParse.colData[parNo].array)[elem] =
 1929                         ((char **)gParse.colData[parNo].array)[elem-1]+alen;
 1930                   ((char **)gParse.colData[parNo].array)[0][0] = '\0';
 1931                } else {
 1932                   free( gParse.colData[parNo].array );
 1933                   *status = MEMORY_ALLOCATION;
 1934                }
 1935             } else {
 1936                *status = MEMORY_ALLOCATION;
 1937             }
 1938          }
 1939          break;
 1940       }
 1941       if( *status ) {
 1942          while( parNo-- ) {
 1943             if( gParse.colData[parNo].datatype==TSTRING )
 1944                FREE( ((char **)gParse.colData[parNo].array)[0] );
 1945             FREE( gParse.colData[parNo].array );
 1946          }
 1947          return( *status );
 1948       }
 1949    }
 1950    
 1951    /**********************************************************************/
 1952    /* Read data from columns needed for the expression and then parse it */
 1953    /**********************************************************************/
 1954    
 1955    if( !uncompress_hkdata( fptr, ntimes, times, status ) ) {
 1956       if( constant ) {
 1957          result = gParse.Nodes[gParse.resultNode].value.data.log;
 1958          elem = ntimes;
 1959          while( elem-- ) time_status[elem] = result;
 1960       } else {
 1961          Info.dataPtr  = time_status;
 1962          Info.nullPtr  = NULL;
 1963          Info.maxRows  = ntimes;
 1964          *status       = parse_data( ntimes, 0, 1, ntimes, gParse.nCols,
 1965                                      gParse.colData, (void*)&Info );
 1966       }
 1967    }
 1968    
 1969    /************/
 1970    /* Clean up */
 1971    /************/
 1972    
 1973    parNo = gParse.nCols;
 1974    while ( parNo-- ) {
 1975       if( gParse.colData[parNo].datatype==TSTRING )
 1976          FREE( ((char **)gParse.colData[parNo].array)[0] );
 1977       FREE( gParse.colData[parNo].array );
 1978    }
 1979    
 1980    if( constant ) gParse.nCols = nCol;
 1981 
 1982    ffcprs();
 1983    return(*status);
 1984 }
 1985 
 1986 /*---------------------------------------------------------------------------*/
 1987 int uncompress_hkdata( fitsfile *fptr,
 1988                        long     ntimes,
 1989                        double   *times,
 1990                        int      *status )
 1991 /*                                                                           */
 1992 /* description                                                               */
 1993 /*---------------------------------------------------------------------------*/
 1994 {
 1995    char parName[256], *sPtr[1], found[1000];
 1996    int parNo, anynul;
 1997    long naxis2, row, currelem;
 1998    double currtime, newtime;
 1999 
 2000    sPtr[0] = parName;
 2001    currelem = 0;
 2002    currtime = -1e38;
 2003 
 2004    parNo=gParse.nCols;
 2005    while( parNo-- ) found[parNo] = 0;
 2006 
 2007    if( ffgkyj( fptr, "NAXIS2", &naxis2, NULL, status ) ) return( *status );
 2008 
 2009    for( row=1; row<=naxis2; row++ ) {
 2010       if( ffgcvd( fptr, gParse.timeCol, row, 1L, 1L, 0.0,
 2011                   &newtime, &anynul, status ) ) return( *status );
 2012       if( newtime != currtime ) {
 2013          /*  New time encountered... propogate parameters to next row  */
 2014          if( currelem==ntimes ) {
 2015             ffpmsg("Found more unique time stamps than caller indicated");
 2016             return( *status = PARSE_BAD_COL );
 2017          }
 2018          times[currelem++] = currtime = newtime;
 2019          parNo = gParse.nCols;
 2020          while( parNo-- ) {
 2021             switch( gParse.colData[parNo].datatype ) {
 2022             case TLONG:
 2023                ((long*)gParse.colData[parNo].array)[currelem] =
 2024                   ((long*)gParse.colData[parNo].array)[currelem-1];
 2025                break;
 2026             case TDOUBLE:
 2027                ((double*)gParse.colData[parNo].array)[currelem] =
 2028                   ((double*)gParse.colData[parNo].array)[currelem-1];
 2029                break;
 2030             case TSTRING:
 2031                strcpy( ((char **)gParse.colData[parNo].array)[currelem],
 2032                        ((char **)gParse.colData[parNo].array)[currelem-1] );
 2033                break;
 2034             }
 2035          }
 2036       }
 2037 
 2038       if( ffgcvs( fptr, gParse.parCol, row, 1L, 1L, "",
 2039                   sPtr, &anynul, status ) ) return( *status );
 2040       parNo = gParse.nCols;
 2041       while( parNo-- )
 2042          if( !fits_strcasecmp( parName, gParse.varData[parNo].name ) ) break;
 2043 
 2044       if( parNo>=0 ) {
 2045          found[parNo] = 1; /* Flag this parameter as found */
 2046          switch( gParse.colData[parNo].datatype ) {
 2047          case TLONG:
 2048             ffgcvj( fptr, gParse.valCol, row, 1L, 1L,
 2049                     ((long*)gParse.colData[parNo].array)[0],
 2050                     ((long*)gParse.colData[parNo].array)+currelem,
 2051                     &anynul, status );
 2052             break;
 2053          case TDOUBLE:
 2054             ffgcvd( fptr, gParse.valCol, row, 1L, 1L,
 2055                     ((double*)gParse.colData[parNo].array)[0],
 2056                     ((double*)gParse.colData[parNo].array)+currelem,
 2057                     &anynul, status );
 2058             break;
 2059          case TSTRING:
 2060             ffgcvs( fptr, gParse.valCol, row, 1L, 1L,
 2061                     ((char**)gParse.colData[parNo].array)[0],
 2062                     ((char**)gParse.colData[parNo].array)+currelem,
 2063                     &anynul, status );
 2064             break;
 2065          }
 2066          if( *status ) return( *status );
 2067       }
 2068    }
 2069 
 2070    if( currelem<ntimes ) {
 2071       ffpmsg("Found fewer unique time stamps than caller indicated");
 2072       return( *status = PARSE_BAD_COL );
 2073    }
 2074 
 2075    /*  Check for any parameters which were not located in the table  */
 2076    parNo = gParse.nCols;
 2077    while( parNo-- )
 2078       if( !found[parNo] ) {
 2079          snprintf( parName, 256, "Parameter not found: %-30s", 
 2080                   gParse.varData[parNo].name );
 2081          ffpmsg( parName );
 2082          *status = PARSE_SYNTAX_ERR;
 2083       }
 2084    return( *status );
 2085 }
 2086 
 2087 /*---------------------------------------------------------------------------*/
 2088 int ffffrw( fitsfile *fptr,         /* I - Input FITS file                   */
 2089             char     *expr,         /* I - Boolean expression                */
 2090             long     *rownum,       /* O - First row of table to eval to T   */
 2091             int      *status )      /* O - Error status                      */
 2092 /*                                                                           */
 2093 /* Evaluate a boolean expression, returning the row number of the first      */
 2094 /* row which evaluates to TRUE                                               */
 2095 /*---------------------------------------------------------------------------*/
 2096 {
 2097    int naxis, constant, dtype;
 2098    long nelem, naxes[MAXDIMS];
 2099    char result;
 2100 
 2101    if( *status ) return( *status );
 2102 
 2103    FFLOCK;
 2104    if( ffiprs( fptr, 0, expr, MAXDIMS, &dtype, &nelem, &naxis,
 2105                naxes, status ) ) {
 2106       ffcprs();
 2107       FFUNLOCK;
 2108       return( *status );
 2109    }
 2110    if( nelem<0 ) {
 2111       constant = 1;
 2112       nelem = -nelem;
 2113    } else
 2114       constant = 0;
 2115 
 2116    if( dtype!=TLOGICAL || nelem!=1 ) {
 2117       ffcprs();
 2118       ffpmsg("Expression does not evaluate to a logical scalar.");
 2119       FFUNLOCK;
 2120       return( *status = PARSE_BAD_TYPE );
 2121    }
 2122 
 2123    *rownum = 0;
 2124    if( constant ) { /* No need to call parser... have result from ffiprs */
 2125       result = gParse.Nodes[gParse.resultNode].value.data.log;
 2126       if( result ) {
 2127          /*  Make sure there is at least 1 row in table  */
 2128          ffgnrw( fptr, &nelem, status );
 2129          if( nelem )
 2130             *rownum = 1;
 2131       }
 2132    } else {
 2133       if( ffiter( gParse.nCols, gParse.colData, 0, 0,
 2134                   ffffrw_work, (void*)rownum, status ) == -1 )
 2135          *status = 0;  /* -1 indicates exitted without error before end... OK */
 2136    }
 2137 
 2138    ffcprs();
 2139    FFUNLOCK;
 2140    return(*status);
 2141 }
 2142 
 2143 /*---------------------------------------------------------------------------*/
 2144 int ffffrw_work(long        totalrows, /* I - Total rows to be processed     */
 2145                 long        offset,    /* I - Number of rows skipped at start*/
 2146                 long        firstrow,  /* I - First row of this iteration    */
 2147                 long        nrows,     /* I - Number of rows in this iter    */
 2148                 int         nCols,     /* I - Number of columns in use       */
 2149                 iteratorCol *colData,  /* IO- Column information/data        */
 2150                 void        *userPtr ) /* I - Data handling instructions     */
 2151 /*                                                                           */
 2152 /* Iterator work function which calls the parser and searches for the        */
 2153 /* first row which evaluates to TRUE.                                        */
 2154 /*---------------------------------------------------------------------------*/
 2155 {
 2156     long idx;
 2157     Node *result;
 2158 
 2159     Evaluate_Parser( firstrow, nrows );
 2160 
 2161     if( !gParse.status ) {
 2162 
 2163        result = gParse.Nodes + gParse.resultNode;
 2164        if( result->operation==CONST_OP ) {
 2165 
 2166           if( result->value.data.log ) {
 2167              *(long*)userPtr = firstrow;
 2168              return( -1 );
 2169           }
 2170 
 2171        } else {
 2172 
 2173           for( idx=0; idx<nrows; idx++ )
 2174              if( result->value.data.logptr[idx] && !result->value.undef[idx] ) {
 2175                 *(long*)userPtr = firstrow + idx;
 2176                 return( -1 );
 2177              }
 2178        }
 2179     }
 2180 
 2181     return( gParse.status );
 2182 }
 2183 
 2184 
 2185 static int set_image_col_types (fitsfile * fptr, const char * name, int bitpix,
 2186                 DataInfo * varInfo, iteratorCol *colIter) {
 2187 
 2188    int istatus;
 2189    double tscale, tzero;
 2190    char temp[80];
 2191 
 2192    switch (bitpix) {
 2193       case BYTE_IMG:
 2194       case SHORT_IMG:
 2195       case LONG_IMG:
 2196          istatus = 0;
 2197          if (fits_read_key(fptr, TDOUBLE, "BZERO", &tzero, NULL, &istatus))
 2198             tzero = 0.0;
 2199 
 2200          istatus = 0;
 2201          if (fits_read_key(fptr, TDOUBLE, "BSCALE", &tscale, NULL, &istatus))
 2202             tscale = 1.0;
 2203 
 2204          if (tscale == 1.0 && (tzero == 0.0 || tzero == 32768.0 )) {
 2205             varInfo->type     = LONG;
 2206             colIter->datatype = TLONG;
 2207          }
 2208          else {
 2209             varInfo->type     = DOUBLE;
 2210             colIter->datatype = TDOUBLE;
 2211             if (DEBUG_PIXFILTER)
 2212                 printf("use DOUBLE for %s with BSCALE=%g/BZERO=%g\n",
 2213                         name, tscale, tzero);
 2214          }
 2215          break;
 2216 
 2217       case LONGLONG_IMG:
 2218       case FLOAT_IMG:
 2219       case DOUBLE_IMG:
 2220          varInfo->type     = DOUBLE;
 2221          colIter->datatype = TDOUBLE;
 2222          break;
 2223       default:
 2224          snprintf(temp, 80,"set_image_col_types: unrecognized image bitpix [%d]\n",
 2225                 bitpix);
 2226          ffpmsg(temp);
 2227          return gParse.status = PARSE_BAD_TYPE;
 2228    }
 2229    return 0;
 2230 }
 2231 
 2232 
 2233 /*************************************************************************
 2234 
 2235         Functions used by the evaluator to access FITS data
 2236             (find_column, find_keywd, allocateCol, load_column)
 2237 
 2238  *************************************************************************/
 2239 
 2240 static int find_column( char *colName, void *itslval )
 2241 {
 2242    FFSTYPE *thelval = (FFSTYPE*)itslval;
 2243    int col_cnt, status;
 2244    int colnum, typecode, type;
 2245    long repeat, width;
 2246    fitsfile *fptr;
 2247    char temp[80];
 2248    double tzero,tscale;
 2249    int istatus;
 2250    DataInfo *varInfo;
 2251    iteratorCol *colIter;
 2252 
 2253 if (DEBUG_PIXFILTER)
 2254    printf("find_column(%s)\n", colName);
 2255 
 2256    if( *colName == '#' )
 2257       return( find_keywd( colName + 1, itslval ) );
 2258 
 2259    fptr = gParse.def_fptr;
 2260 
 2261    status = 0;
 2262    col_cnt = gParse.nCols;
 2263 
 2264 if (gParse.hdutype == IMAGE_HDU) {
 2265    int i;
 2266    if (!gParse.pixFilter) {
 2267       gParse.status = COL_NOT_FOUND;
 2268       ffpmsg("find_column: IMAGE_HDU but no PixelFilter");
 2269       return pERROR;
 2270    }
 2271 
 2272    colnum = -1;
 2273    for (i = 0; i < gParse.pixFilter->count; ++i) {
 2274       if (!fits_strcasecmp(colName, gParse.pixFilter->tag[i]))
 2275          colnum = i;
 2276    }
 2277    if (colnum < 0) {
 2278       snprintf(temp, 80, "find_column: PixelFilter tag %s not found", colName);
 2279       ffpmsg(temp);
 2280       gParse.status = COL_NOT_FOUND;
 2281       return pERROR;
 2282    }
 2283 
 2284    if( allocateCol( col_cnt, &gParse.status ) ) return pERROR;
 2285 
 2286    varInfo = gParse.varData + col_cnt;
 2287    colIter = gParse.colData + col_cnt;
 2288 
 2289    fptr = gParse.pixFilter->ifptr[colnum];
 2290    fits_get_img_param(fptr,
 2291                 MAXDIMS,
 2292                 &typecode, /* actually bitpix */
 2293                 &varInfo->naxis,
 2294                 &varInfo->naxes[0],
 2295                 &status);
 2296    varInfo->nelem = 1;
 2297    type = COLUMN;
 2298    if (set_image_col_types(fptr, colName, typecode, varInfo, colIter))
 2299       return pERROR;
 2300    colIter->fptr = fptr;
 2301    colIter->iotype = InputCol;
 2302 }
 2303 else { /* HDU holds a table */
 2304    if( gParse.compressed )
 2305       colnum = gParse.valCol;
 2306    else
 2307       if( fits_get_colnum( fptr, CASEINSEN, colName, &colnum, &status ) ) {
 2308          if( status == COL_NOT_FOUND ) {
 2309             type = find_keywd( colName, itslval );
 2310             if( type != pERROR ) ffcmsg();
 2311             return( type );
 2312          }
 2313          gParse.status = status;
 2314          return pERROR;
 2315       }
 2316    
 2317    if( fits_get_coltype( fptr, colnum, &typecode,
 2318                          &repeat, &width, &status ) ) {
 2319       gParse.status = status;
 2320       return pERROR;
 2321    }
 2322 
 2323    if( allocateCol( col_cnt, &gParse.status ) ) return pERROR;
 2324 
 2325    varInfo = gParse.varData + col_cnt;
 2326    colIter = gParse.colData + col_cnt;
 2327 
 2328    fits_iter_set_by_num( colIter, fptr, colnum, 0, InputCol );
 2329 }
 2330 
 2331    /*  Make sure we don't overflow variable name array  */
 2332    strncpy(varInfo->name,colName,MAXVARNAME);
 2333    varInfo->name[MAXVARNAME] = '\0';
 2334 
 2335 if (gParse.hdutype != IMAGE_HDU) {
 2336    switch( typecode ) {
 2337    case TBIT:
 2338       varInfo->type     = BITSTR;
 2339       colIter->datatype = TBYTE;
 2340       type = BITCOL;
 2341       break;
 2342    case TBYTE:
 2343    case TSHORT:
 2344    case TLONG:
 2345       /* The datatype of column with TZERO and TSCALE keywords might be 
 2346          float or double. 
 2347       */
 2348       snprintf(temp,80,"TZERO%d",colnum);
 2349       istatus = 0;
 2350       if(fits_read_key(fptr,TDOUBLE,temp,&tzero,NULL,&istatus)) {
 2351           tzero = 0.0;
 2352       } 
 2353       snprintf(temp,80,"TSCAL%d",colnum);
 2354       istatus = 0;
 2355       if(fits_read_key(fptr,TDOUBLE,temp,&tscale,NULL,&istatus)) {
 2356           tscale = 1.0;
 2357       } 
 2358       if (tscale == 1.0 && (tzero == 0.0 || tzero == 32768.0 )) {
 2359           varInfo->type     = LONG;
 2360           colIter->datatype = TLONG;
 2361 /*    Reading an unsigned long column as a long can cause overflow errors.
 2362       Treat the column as a double instead.
 2363       } else if (tscale == 1.0 &&  tzero == 2147483648.0 ) {
 2364           varInfo->type     = LONG;
 2365           colIter->datatype = TULONG;
 2366  */
 2367 
 2368       }
 2369       else {
 2370           varInfo->type     = DOUBLE;
 2371           colIter->datatype = TDOUBLE;
 2372       }
 2373       type = COLUMN;
 2374       break;
 2375 /* 
 2376   For now, treat 8-byte integer columns as type double.
 2377   This can lose precision, so the better long term solution
 2378   will be to add support for TLONGLONG as a separate datatype.
 2379 */
 2380    case TLONGLONG:
 2381    case TFLOAT:
 2382    case TDOUBLE:
 2383       varInfo->type     = DOUBLE;
 2384       colIter->datatype = TDOUBLE;
 2385       type = COLUMN;
 2386       break;
 2387    case TLOGICAL:
 2388       varInfo->type     = BOOLEAN;
 2389       colIter->datatype = TLOGICAL;
 2390       type = BCOLUMN;
 2391       break;
 2392    case TSTRING:
 2393       varInfo->type     = STRING;
 2394       colIter->datatype = TSTRING;
 2395       type = SCOLUMN;
 2396       if ( width >= MAX_STRLEN ) {
 2397     snprintf(temp, 80, "column %d is wider than maximum %d characters",
 2398         colnum, MAX_STRLEN-1);
 2399         ffpmsg(temp);
 2400     gParse.status = PARSE_LRG_VECTOR;
 2401     return pERROR;
 2402       }
 2403       if( gParse.hdutype == ASCII_TBL ) repeat = width;
 2404       break;
 2405    default:
 2406       if (typecode < 0) {
 2407         snprintf(temp, 80,"variable-length array columns are not supported. typecode = %d", typecode);
 2408         ffpmsg(temp);
 2409       }
 2410       gParse.status = PARSE_BAD_TYPE;
 2411       return pERROR;
 2412    }
 2413    varInfo->nelem = repeat;
 2414    if( repeat>1 && typecode!=TSTRING ) {
 2415       if( fits_read_tdim( fptr, colnum, MAXDIMS,
 2416                           &varInfo->naxis,
 2417                           &varInfo->naxes[0], &status )
 2418           ) {
 2419          gParse.status = status;
 2420          return pERROR;
 2421       }
 2422    } else {
 2423       varInfo->naxis = 1;
 2424       varInfo->naxes[0] = 1;
 2425    }
 2426 }
 2427    gParse.nCols++;
 2428    thelval->lng = col_cnt;
 2429 
 2430    return( type );
 2431 }
 2432 
 2433 static int find_keywd(char *keyname, void *itslval )
 2434 {
 2435    FFSTYPE *thelval = (FFSTYPE*)itslval;
 2436    int status, type;
 2437    char keyvalue[FLEN_VALUE], dtype;
 2438    fitsfile *fptr;
 2439    double rval;
 2440    int bval;
 2441    long ival;
 2442 
 2443    status = 0;
 2444    fptr = gParse.def_fptr;
 2445    if( fits_read_keyword( fptr, keyname, keyvalue, NULL, &status ) ) {
 2446       if( status == KEY_NO_EXIST ) {
 2447          /*  Do this since ffgkey doesn't put an error message on stack  */
 2448          snprintf(keyvalue,FLEN_VALUE, "ffgkey could not find keyword: %s",keyname);
 2449          ffpmsg(keyvalue);
 2450       }
 2451       gParse.status = status;
 2452       return( pERROR );
 2453    }
 2454       
 2455    if( fits_get_keytype( keyvalue, &dtype, &status ) ) {
 2456       gParse.status = status;
 2457       return( pERROR );
 2458    }
 2459       
 2460    switch( dtype ) {
 2461    case 'C':
 2462       fits_read_key_str( fptr, keyname, keyvalue, NULL, &status );
 2463       type = STRING;
 2464       strcpy( thelval->str , keyvalue );
 2465       break;
 2466    case 'L':
 2467       fits_read_key_log( fptr, keyname, &bval, NULL, &status );
 2468       type = BOOLEAN;
 2469       thelval->log = bval;
 2470       break;
 2471    case 'I':
 2472       fits_read_key_lng( fptr, keyname, &ival, NULL, &status );
 2473       type = LONG;
 2474       thelval->lng = ival;
 2475       break;
 2476    case 'F':
 2477       fits_read_key_dbl( fptr, keyname, &rval, NULL, &status );
 2478       type = DOUBLE;
 2479       thelval->dbl = rval;
 2480       break;
 2481    default:
 2482       type = pERROR;
 2483       break;
 2484    }
 2485 
 2486    if( status ) {
 2487       gParse.status=status;
 2488       return pERROR;
 2489    }
 2490 
 2491    return( type );
 2492 }
 2493 
 2494 static int allocateCol( int nCol, int *status )
 2495 {
 2496    if( (nCol%25)==0 ) {
 2497       if( nCol ) {
 2498          gParse.colData  = (iteratorCol*) realloc( gParse.colData,
 2499                                               (nCol+25)*sizeof(iteratorCol) );
 2500          gParse.varData  = (DataInfo   *) realloc( gParse.varData,
 2501                                               (nCol+25)*sizeof(DataInfo)    );
 2502       } else {
 2503          gParse.colData  = (iteratorCol*) malloc( 25*sizeof(iteratorCol) );
 2504          gParse.varData  = (DataInfo   *) malloc( 25*sizeof(DataInfo)    );
 2505       }
 2506       if(    gParse.colData  == NULL
 2507           || gParse.varData  == NULL    ) {
 2508          if( gParse.colData  ) free(gParse.colData);
 2509          if( gParse.varData  ) free(gParse.varData);
 2510          gParse.colData = NULL;
 2511          gParse.varData = NULL;
 2512          return( *status = MEMORY_ALLOCATION );
 2513       }
 2514    }
 2515    gParse.varData[nCol].data  = NULL;
 2516    gParse.varData[nCol].undef = NULL;
 2517    return 0;
 2518 }
 2519 
 2520 static int load_column( int varNum, long fRow, long nRows,
 2521                         void *data, char *undef )
 2522 {
 2523    iteratorCol *var = gParse.colData+varNum;
 2524    long nelem,nbytes,row,len,idx;
 2525    char **bitStrs, msg[80];
 2526    unsigned char *bytes;
 2527    int status = 0, anynul;
 2528 
 2529   if (gParse.hdutype == IMAGE_HDU) {
 2530     /* This test would need to be on a per varNum basis to support
 2531      * cross HDU operations */
 2532     fits_read_imgnull(var->fptr, var->datatype, fRow, nRows,
 2533                 data, undef, &anynul, &status);
 2534     if (DEBUG_PIXFILTER)
 2535         printf("load_column: IMAGE_HDU fRow=%ld, nRows=%ld => %d\n",
 2536                         fRow, nRows, status);
 2537   } else { 
 2538 
 2539    nelem = nRows * var->repeat;
 2540 
 2541    switch( var->datatype ) {
 2542    case TBYTE:
 2543       nbytes = ((var->repeat+7)/8) * nRows;
 2544       bytes = (unsigned char *)malloc( nbytes * sizeof(char) );
 2545 
 2546       ffgcvb(var->fptr, var->colnum, fRow, 1L, nbytes,
 2547              0, bytes, &anynul, &status);
 2548 
 2549       nelem = var->repeat;
 2550       bitStrs = (char **)data;
 2551       for( row=0; row<nRows; row++ ) {
 2552          idx = (row)*( (nelem+7)/8 ) + 1;
 2553          for(len=0; len<nelem; len++) {
 2554             if( bytes[idx] & (1<<(7-len%8)) )
 2555                bitStrs[row][len] = '1';
 2556             else
 2557                bitStrs[row][len] = '0';
 2558             if( len%8==7 ) idx++;
 2559          }
 2560          bitStrs[row][len] = '\0';
 2561       }
 2562 
 2563       FREE( (char *)bytes );
 2564       break;
 2565    case TSTRING:
 2566       ffgcfs(var->fptr, var->colnum, fRow, 1L, nRows,
 2567              (char **)data, undef, &anynul, &status);
 2568       break;
 2569    case TLOGICAL:
 2570       ffgcfl(var->fptr, var->colnum, fRow, 1L, nelem,
 2571              (char *)data, undef, &anynul, &status);
 2572       break;
 2573    case TLONG:
 2574       ffgcfj(var->fptr, var->colnum, fRow, 1L, nelem,
 2575              (long *)data, undef, &anynul, &status);
 2576       break;
 2577    case TDOUBLE:
 2578       ffgcfd(var->fptr, var->colnum, fRow, 1L, nelem,
 2579              (double *)data, undef, &anynul, &status);
 2580       break;
 2581    default:
 2582       snprintf(msg,80,"load_column: unexpected datatype %d", var->datatype);
 2583       ffpmsg(msg);
 2584    }
 2585   }
 2586    if( status ) {
 2587       gParse.status = status;
 2588       return pERROR;
 2589    }
 2590 
 2591    return 0;
 2592 }
 2593 
 2594 
 2595 /*--------------------------------------------------------------------------*/
 2596 int fits_pixel_filter (PixelFilter * filter, int * status)
 2597 /* Evaluate an expression using the data in the input FITS file(s)          */
 2598 /*--------------------------------------------------------------------------*/
 2599 {
 2600    parseInfo Info = { 0 };
 2601    int naxis, bitpix;
 2602    long nelem, naxes[MAXDIMS];
 2603    int col_cnt;
 2604    Node *result;
 2605    int datatype;
 2606    fitsfile * infptr;
 2607    fitsfile * outfptr;
 2608    char * DEFAULT_TAGS[] = { "X" };
 2609    char msg[256];
 2610    int writeBlankKwd = 0;   /* write BLANK if any output nulls? */
 2611 
 2612    DEBUG_PIXFILTER = getenv("DEBUG_PIXFILTER") ? 1 : 0;
 2613 
 2614    if (*status)
 2615       return (*status);
 2616 
 2617    FFLOCK;
 2618    if (!filter->tag || !filter->tag[0] || !filter->tag[0][0]) {
 2619       filter->tag = DEFAULT_TAGS;
 2620       if (DEBUG_PIXFILTER)
 2621          printf("using default tag '%s'\n", filter->tag[0]);
 2622    }
 2623 
 2624    infptr = filter->ifptr[0];
 2625    outfptr = filter->ofptr;
 2626    gParse.pixFilter = filter;
 2627 
 2628    if (ffiprs(infptr, 0, filter->expression, MAXDIMS,
 2629             &Info.datatype, &nelem, &naxis, naxes, status)) {
 2630       goto CLEANUP;
 2631    }
 2632 
 2633    if (nelem < 0) {
 2634       nelem = -nelem;
 2635    }
 2636 
 2637    {
 2638       /* validate result type */
 2639       const char * type = 0;
 2640       switch (Info.datatype) {
 2641          case TLOGICAL:  type = "LOGICAL"; break;
 2642          case TLONG:     type = "LONG"; break;
 2643          case TDOUBLE:   type = "DOUBLE"; break;
 2644          case TSTRING:   type = "STRING";
 2645                          *status = pERROR;
 2646                          ffpmsg("pixel_filter: cannot have string image");
 2647          case TBIT:      type = "BIT";
 2648                          if (DEBUG_PIXFILTER)
 2649                             printf("hmm, image from bits?\n");
 2650                          break;
 2651          default:       type = "UNKNOWN?!";
 2652                         *status = pERROR;
 2653                         ffpmsg("pixel_filter: unexpected result datatype");
 2654       }
 2655       if (DEBUG_PIXFILTER)
 2656          printf("result type is %s [%d]\n", type, Info.datatype);
 2657       if (*status)
 2658          goto CLEANUP;
 2659    }
 2660 
 2661    if (fits_get_img_param(infptr, MAXDIMS,
 2662             &bitpix, &naxis, &naxes[0], status)) {
 2663       ffpmsg("pixel_filter: unable to read input image parameters");
 2664       goto CLEANUP;
 2665    }
 2666 
 2667    if (DEBUG_PIXFILTER)
 2668       printf("input bitpix %d\n", bitpix);
 2669 
 2670    if (Info.datatype == TDOUBLE) {
 2671        /*  for floating point expressions, set the default output image to
 2672            bitpix = -32 (float) unless the default is already a double */
 2673        if (bitpix != DOUBLE_IMG)
 2674            bitpix = FLOAT_IMG;
 2675    }
 2676 
 2677    /* override output image bitpix if specified by caller */
 2678    if (filter->bitpix)
 2679       bitpix = filter->bitpix;
 2680    if (DEBUG_PIXFILTER)
 2681       printf("output bitpix %d\n", bitpix);
 2682 
 2683    if (fits_create_img(outfptr, bitpix, naxis, naxes, status)) {
 2684       ffpmsg("pixel_filter: unable to create output image");
 2685       goto CLEANUP;
 2686    }
 2687 
 2688    /* transfer keycards */
 2689    {
 2690       int i, ncards, more;
 2691       if (fits_get_hdrspace(infptr, &ncards, &more, status)) {
 2692          ffpmsg("pixel_filter: unable to determine number of keycards");
 2693          goto CLEANUP;
 2694       }
 2695 
 2696       for (i = 1; i <= ncards; ++i) {
 2697 
 2698          int keyclass;
 2699          char card[FLEN_CARD];
 2700 
 2701          if (fits_read_record(infptr, i, card, status)) {
 2702             snprintf(msg, 256,"pixel_filter: unable to read keycard %d", i);
 2703             ffpmsg(msg);
 2704             goto CLEANUP;
 2705          }
 2706 
 2707          keyclass = fits_get_keyclass(card);
 2708          if (keyclass == TYP_STRUC_KEY) {
 2709             /* output structure defined by fits_create_img */
 2710          }
 2711          else if (keyclass == TYP_COMM_KEY && i < 12) {
 2712             /* assume this is one of the FITS standard comments */
 2713          }
 2714          else if (keyclass == TYP_NULL_KEY && bitpix < 0) {
 2715             /* do not transfer BLANK to real output image */
 2716          }
 2717          else if (keyclass == TYP_SCAL_KEY && bitpix < 0) {
 2718             /* do not transfer BZERO, BSCALE to real output image */
 2719          }
 2720          else if (fits_write_record(outfptr, card, status)) {
 2721             snprintf(msg,256, "pixel_filter: unable to write keycard '%s' [%d]\n",
 2722                         card, *status);
 2723             ffpmsg(msg);
 2724             goto CLEANUP;
 2725          }
 2726       }
 2727    }
 2728 
 2729    switch (bitpix) {
 2730       case BYTE_IMG: datatype = TLONG; Info.datatype = TBYTE; break;
 2731       case SHORT_IMG: datatype = TLONG; Info.datatype = TSHORT; break;
 2732       case LONG_IMG: datatype = TLONG; Info.datatype = TLONG; break;
 2733       case FLOAT_IMG: datatype = TDOUBLE; Info.datatype = TFLOAT; break;
 2734       case DOUBLE_IMG: datatype = TDOUBLE; Info.datatype = TDOUBLE; break;
 2735 
 2736       default:
 2737            snprintf(msg, 256,"pixel_filter: unexpected output bitpix %d\n", bitpix);
 2738            ffpmsg(msg);
 2739            *status = pERROR;
 2740            goto CLEANUP;
 2741    }
 2742 
 2743    if (bitpix > 0) { /* arrange for NULLs in output */
 2744       long nullVal = filter->blank;
 2745       if (!filter->blank) {
 2746          int tstatus = 0;
 2747          if (fits_read_key_lng(infptr, "BLANK", &nullVal, 0, &tstatus)) {
 2748 
 2749             writeBlankKwd = 1;
 2750 
 2751             if (bitpix == BYTE_IMG)
 2752                 nullVal = UCHAR_MAX;
 2753             else if (bitpix == SHORT_IMG)
 2754                 nullVal = SHRT_MIN;
 2755             else if (bitpix == LONG_IMG) {
 2756                if (sizeof(long) == 8 && sizeof(int) == 4)
 2757                   nullVal = INT_MIN;
 2758                else
 2759                   nullVal = LONG_MIN;
 2760             }
 2761             else
 2762                 printf("unhandled positive output BITPIX %d\n", bitpix);
 2763          }
 2764 
 2765          filter->blank = nullVal;
 2766       }
 2767 
 2768       fits_set_imgnull(outfptr, filter->blank, status);
 2769       if (DEBUG_PIXFILTER)
 2770          printf("using blank %ld\n", nullVal);
 2771 
 2772    }
 2773 
 2774    if (!filter->keyword[0]) {
 2775       iteratorCol * colIter;
 2776       DataInfo * varInfo;
 2777 
 2778       /*************************************/
 2779       /* Create new iterator Output Column */
 2780       /*************************************/
 2781       col_cnt = gParse.nCols;
 2782       if (allocateCol(col_cnt, status))
 2783          goto CLEANUP;
 2784       gParse.nCols++;
 2785 
 2786       colIter = &gParse.colData[col_cnt];
 2787       colIter->fptr = filter->ofptr;
 2788       colIter->iotype = OutputCol;
 2789       varInfo = &gParse.varData[col_cnt];
 2790       set_image_col_types(colIter->fptr, "CREATED", bitpix, varInfo, colIter);
 2791 
 2792       Info.maxRows = -1;
 2793 
 2794       if (ffiter(gParse.nCols, gParse.colData, 0,
 2795                      0, parse_data, &Info, status) == -1)
 2796             *status = 0;
 2797       else if (*status)
 2798          goto CLEANUP;
 2799 
 2800       if (Info.anyNull) {
 2801          if (writeBlankKwd) {
 2802             fits_update_key_lng(outfptr, "BLANK", filter->blank, "NULL pixel value", status);
 2803             if (*status)
 2804                 ffpmsg("pixel_filter: unable to write BLANK keyword");
 2805             if (DEBUG_PIXFILTER) {
 2806                 printf("output has NULLs\n");
 2807                 printf("wrote blank [%d]\n", *status);
 2808             }
 2809          }
 2810       }
 2811       else if (bitpix > 0) /* never used a null */
 2812          if (fits_set_imgnull(outfptr, -1234554321, status))
 2813             ffpmsg("pixel_filter: unable to reset imgnull");
 2814    }
 2815    else {
 2816 
 2817       /* Put constant result into keyword */
 2818       char * parName = filter->keyword;
 2819       char * parInfo = filter->comment;
 2820 
 2821       result  = gParse.Nodes + gParse.resultNode;
 2822       switch (Info.datatype) {
 2823       case TDOUBLE:
 2824          ffukyd(outfptr, parName, result->value.data.dbl, 15, parInfo, status);
 2825          break;
 2826       case TLONG:
 2827          ffukyj(outfptr, parName, result->value.data.lng, parInfo, status);
 2828          break;
 2829       case TLOGICAL:
 2830          ffukyl(outfptr, parName, result->value.data.log, parInfo, status);
 2831          break;
 2832       case TBIT:
 2833       case TSTRING:
 2834          ffukys(outfptr, parName, result->value.data.str, parInfo, status);
 2835          break;
 2836       default:
 2837          snprintf(msg, 256,"pixel_filter: unexpected constant result type [%d]\n",
 2838                 Info.datatype);
 2839          ffpmsg(msg);
 2840       }
 2841    }
 2842 
 2843 CLEANUP:
 2844    ffcprs();
 2845    FFUNLOCK;
 2846    return (*status);
 2847 }