00013 {
00014    PARSER_code * pcode = NULL ;
00015    char sym[4] ;
00016    double atoz[26] , value , del=1.0 ;
00017    int ii,jj , kvar , nopt , qvar , num=-1 , verbose=0 ;
00018    MRI_IMAGE * inim[26] ;
00019    float     * inar[26] ;
00020    MRI_IMAGE *dindex_im = NULL;
00021    float     *dindex = NULL;
00022    char abet[] = "abcdefghijklmnopqrstuvwxyz" ;
00023 
00024    /*-- help? --*/
00025 
00026    if( argc < 3 ){
00027       printf("Usage: 1deval [options] -expr 'expression'\n"
00028              "Evaluates the expression and writes the result to stdout.\n"
00029              "Any single letter from a-z can be used as the independent\n"
00030              "variable in the expression.\n"
00031              "\n"
00032              "Options:\n"
00033              "  -del d   = Use 'd' as the step for the variable in the\n"
00034              "               expression [default = 1.0]\n"
00035              "  -num n   = Evaluate the expression 'n' times.\n"
00036              "               If -num is not used, then the length of an\n"
00037              "               input time series is used.  If there are no\n"
00038              "               time series input, then -num is required.\n"
00039              "  -a q.1D  = Read time series file q.1D and assign it\n"
00040              "               to the symbol 'a' (as in 3dcalc).\n"
00041              "  -index i.1D = Read index column from file i.1D and\n"
00042              "                 write it out as 1st column of output.\n"
00043              "                 This option is useful when working with\n"
00044              "                 surface data.\n" 
00045              "Examples:\n"
00046              "  1deval -expr 'sin(2*PI*t)' -del 0.01 -num 101 > sin.1D\n"
00047              "  1deval -expr 'a*b*x' -a fred.1D -b ethel.1D > x.1D\n"
00048 
00049              "\n"
00050              TS_HELP_STRING
00051             ) ;
00052       exit(0) ;
00053    }
00054 
00055    machdep() ;
00056 
00057    /*-- initialize --*/
00058 
00059    for( ii=0 ; ii < 26 ; ii++ ){
00060       atoz[ii] = 0.0 ; inim[ii] = NULL ; inar[ii] = NULL ;
00061    }
00062 
00063    /*-- read options --*/
00064 
00065    nopt = 1 ;
00066    while( nopt < argc ){
00067 
00068       if( strcmp(argv[nopt],"-verb") == 0 ){
00069          verbose++ ;
00070          nopt++ ; continue ;
00071       }
00072 
00073       if( strlen(argv[nopt]) == 2   && argv[nopt][0] == '-' &&
00074           argv[nopt][1]      >= 'a' && argv[nopt][1] <= 'z'   ){
00075 
00076          int ival = argv[nopt][1] - 'a' ;
00077 
00078          if( inim[ival] != NULL ){
00079             fprintf(stderr,"** Can't define symbol %c twice!\n",argv[nopt][1]);
00080             exit(1) ;
00081          }
00082 
00083          nopt++ ;
00084          if( nopt >= argc ){
00085             fprintf(stderr,"** -%c needs an argument!\n",argv[nopt][1]); exit(1);
00086          }
00087 
00088          inim[ival] = mri_read_1D( argv[nopt] ) ;
00089          if( inim[ival] == NULL ){
00090             fprintf(stderr,"** Can't read time series file %s\n",argv[nopt]);
00091             exit(1) ;
00092          }
00093 
00094          inar[ival] = MRI_FLOAT_PTR(inim[ival]) ;
00095          nopt++ ; continue ;
00096       }
00097 
00098       if( strcmp(argv[nopt],"-expr") == 0 ){
00099          if( pcode != NULL ){
00100             fprintf(stderr,"** Can't have 2 -expr options!\n") ;
00101             exit(1) ;
00102          }
00103          nopt++ ;
00104          if( nopt >= argc ){
00105             fprintf(stderr,"** -expr needs an argument!\n") ; exit(1) ;
00106          }
00107          pcode = PARSER_generate_code( argv[nopt] ) ;  /* compile */
00108          if( pcode == NULL ){
00109             fprintf(stderr,"** Illegal expression!\n") ;
00110             exit(1) ;
00111          }
00112          nopt++ ; continue ;
00113       }
00114 
00115       if( strcmp(argv[nopt],"-del") == 0 ){
00116          nopt++ ;
00117          if( nopt >= argc ){
00118             fprintf(stderr,"** -del needs an argument!\n") ;
00119             exit(1) ;
00120          }
00121          del = strtod( argv[nopt] , NULL ) ;
00122          if( del == 0 ){
00123             fprintf(stderr,"** -del value must not be zero!\n") ;
00124             exit(1) ;
00125          }
00126          if( verbose ) fprintf(stderr,"del set to %g\n",del) ;
00127          nopt++ ; continue ;
00128       }
00129       
00130       if( strcmp(argv[nopt],"-index") == 0 ){
00131          nopt++ ;
00132          if( nopt >= argc ){
00133             fprintf(stderr,"** -index needs an argument!\n") ;
00134             exit(1) ;
00135          }
00136          
00137          dindex_im = mri_read_1D( argv[nopt] ) ;
00138          if( dindex_im == NULL ){
00139             fprintf(stderr,"** Can't read time series file %s\n",argv[nopt]);
00140             exit(1) ;
00141          }
00142          if (dindex_im->ny != 1) {
00143             fprintf(stderr,"** Only one column allowed for indexing.\n   Found %d columns\n", dindex_im->ny);
00144             exit(1) ;
00145          }
00146          dindex = MRI_FLOAT_PTR(dindex_im) ;
00147          nopt++ ; continue ;
00148          
00149       }
00150 
00151 
00152       if( strcmp(argv[nopt],"-num") == 0 ){
00153          nopt++ ;
00154          if( nopt >= argc ){
00155             fprintf(stderr,"** -num needs an argument!\n") ;
00156             exit(1) ;
00157          }
00158          num = strtol( argv[nopt] , NULL , 10 ) ;
00159          if( num <= 0 ){
00160             fprintf(stderr,"** -num value must be positive!\n") ;
00161             exit(1) ;
00162          }
00163          nopt++ ; continue ;
00164       }
00165 
00166       fprintf(stderr,"** %s = unknown command line option!\n",argv[nopt]) ;
00167       exit(1) ;
00168    }
00169 
00170    if( num <= 0 ){
00171       for( ii=0 ; ii < 26 ; ii++ ){
00172          if( inim[ii] != NULL ){ num = inim[ii]->nx ; break ; }
00173       }
00174       if( num > 0 ){
00175          if( verbose )
00176             fprintf(stderr,"++ Set num = %d from input time series\n",num);
00177       } else {
00178          fprintf(stderr,"** Need to supply -num on command line!\n"); exit(1);
00179       }
00180    }
00181 
00182    for( qvar=ii=0 ; ii < 26 ; ii++ ){
00183       if( inim[ii] != NULL && inim[ii]->nx < num ){
00184          fprintf(stderr,"** Time series file %s is too short!\n",inim[ii]->name);
00185          qvar++ ;
00186       }
00187    }
00188    
00189    if (dindex) {
00190       if ( dindex_im->nx != num) {
00191          fprintf(stderr,"** Number of values in index column (%d) \n   not equal to number of values in data (%d)\n", dindex_im->nx, num );
00192          exit(1) ;
00193       }
00194    }
00195    
00196    if( qvar > 0 ) exit(1) ;
00197 
00198    if( pcode == NULL ){
00199       fprintf(stderr,"** -expr is missing!\n") ; exit(1) ;
00200    }
00201 
00202    qvar = 0 ; kvar = -1 ;                       /* find symbol */
00203    for( ii=0 ; ii < 26 ; ii++ ){
00204       sym[0] = 'A' + ii ; sym[1] = '\0' ;
00205       if( PARSER_has_symbol(sym,pcode) ){
00206          if( inim[ii] == NULL ){
00207             qvar++ ; if( kvar < 0 ) kvar = ii ;
00208          }
00209       } else if( inim[ii] != NULL ){
00210          fprintf(stderr,"++ Symbol %c defined but not used!\n",abet[ii]) ;
00211       }
00212    }
00213    if( qvar > 1 ){
00214       fprintf(stderr,"++ Found %d indeterminate symbols in expression,\n"
00215                      "++   but there should only be 0 or 1.\n"
00216                      "++   Will use symbol %c as the variable.\n" ,
00217               qvar , abet[kvar] ) ;
00218    }
00219 
00220    /*-- evaluate --*/
00221 
00222    for( ii=0 ; ii < num ; ii++ ){
00223       for( jj=0 ; jj < 26 ; jj++ )
00224          if( inar[jj] != NULL ) atoz[jj] = inar[jj][ii] ;
00225       if( kvar >= 0 ) atoz[kvar] = ii * del ;
00226       value = PARSER_evaluate_one( pcode , atoz ) ;
00227       if (dindex) printf(" %d\t%g\n", (int)dindex[ii], value) ;
00228       else printf(" %g\n",value) ;
00229    }
00230    exit(0) ;
00231 }