c
c     Copyright (c) 1999,2000, Compaq Computer Corporation
c
c     See the file "license.terms" for information on usage and redistribution
c     of this file, and for a DISCLAIMER OF ALL WARRANTIES.
c
c     Subroutine package to take an expression 'ala' FORTRAN
c     (input as a character string), a set of data values and
c     output a set of result values, after substituting each data value
c     in the expression.
c
c     The expression can have variables in it (e.g. "x" "y")
c     that is a placeholder for the input data value.
c     The variable can be any string (except for predefined functions
c     see below), and defined through the input parameter "var".
c
c     There are many possible errors, the user should check return codes
c     carefully. Iptr points to the place in the input expression
c     where the error occurred. If ierr was not SUCCESS and iptr=zero,
c     then the error occurred at some point in the expression evaluation
c     where it was impossible to backtrack and find where the error occurred.
c
c     Operators:
c     +  == add, also unary plus
c     -  == subtract, also unary minus
c     *  == multiply
c     /  == divide
c     ** == exponentiation
c     >  greater than
c     <  less than
c
c     Functions:
c     sqrt == sqrt
c     log  == log base e
c     exp  == exponential function
c     log10 == log base 10
c     abs  == absolute value
c     sin  == sine function (radian argument)
c     cos  == cosine function (radian argument)
c     tan  == tangent function (radian argument)
c     asin == arcsin function (returns radians)
c     acos == arccos function (returns radians)
c     atan == arctan function (returns radians)
c     atan2== arctan function with two arguments (returns radians)
c     max(...) == maximum function (arb. number of arguments)
c     min(...) == minimum function (arb. number of arguments)
c     sind == sine function (degree argument)
c     cosd == cosine function (degree argument)
c     tand == tangent function (degree argument)
c     merge(a,b,c) == like C's syntax: c ? a : b
c
c     The following functions only make sense when calling evaluate
c     repeatedly on a vector of inputs (ival > 1)
c
c     ran  == random number generator [0..1.0]
c             Caveat: takes an argument to satisfy parser, but does not
c                     use it
c     integrate(y,x) == Integrate y(x) from initial x to current x value
c
c     deriv(y,x)
c     differentiate(y,x) == Differentiate y(x) from initial x to
c             current x value.
c             Caveats: One differentiation call per expression
c
c     step(init,step) == provide stepping function = init + step*(ival-1)
c
      subroutine ev_expression(expr,var,nvars,ival,values,isize,
     1                         result,ierr,ind)
      implicit none
c
      character*(*) expr
      integer nvars, ival, ierr, ind, isize
c     Actual second dimension of values is "nvars", but we don't code
c     this because nvars could be zero.
      double precision    values(isize,*), result(ival)
c     Actual dimension of var is "nvars", but we don't code this because
c     nvars could be zero.
      character*(*) var(*)
c
      include 'evaluate.cmn'
c
      logical ltest
      integer ilen, i, j, k, string_len
      double precision rdum1, rdum2, inputs(EVAL_MAXVARS)
c
      ierr = EVAL_SUCCESS
      ind  = 0
c
      call ev_init_common
      call ev_initialize
c
c     Copy expression and variable names into local storage, upcase.
      ilen = string_len(expr)
      if( ilen .gt. EVAL_MAXEXPR ) then
         ierr = EVAL_LARGE
         goto 10000
      endif
      ev_expr = expr(1:ilen)
      call upcase(ev_expr)
c
      ev_num_vars = nvars
      do 100 i=1,nvars
         call ev_check_varname(var(i),i-1,ierr)
         if( ierr .ne. EVAL_SUCCESS ) goto 10000
         ev_var(i)  = var(i)
         call upcase(ev_var(i))
100   continue
c
c     Break expression into tokens
      call ev_tokenize(ierr,ind)
      if( ierr .ne. EVAL_SUCCESS ) goto 10000
c
c     Convert to reverse polish notation
      call ev_rpn(ierr)
      if( ierr .ne. EVAL_SUCCESS ) goto 10000
c
c     Perform a "dry run" of the expresion evaluator, in order to catch
c     syntax errors.
      ltest = .TRUE.
      call ev_compute_expr(rdum1,rdum2,ltest,ierr)
c
c     Evaluate expression for each input value
      ltest = .FALSE.
c
c     Check for too many integration or differentiation functions
      if( ev_integrate_id .gt. EVAL_STACKSIZE ) then
         ierr = EVAL_TOOMANYINTS
         goto 10000
      else if( ev_differentiate_id .gt. 1 ) then
         ierr = EVAL_TOOMANYDIFFS
         goto 10000
      else if( ev_differentiate_id .eq. 1 ) then
c        Special case to support 2-pass derivatives
         ev_differentiate_npts  = ival
         allocate(ev_differentiate_x(1:ival),stat=ierr)
         if( ierr .ne. 0 ) then
            ierr = EVAL_ALLOCATE
            goto 10000
         endif
         allocate(ev_differentiate_y(1:ival),stat=ierr)
         if( ierr .ne. 0 ) then
            deallocate(ev_differentiate_x)
            ierr = EVAL_ALLOCATE
            goto 10000
         endif
         do k=1,2
            do i=1,ival
               ev_count = i
               do j=1,nvars
                  inputs(j) = values(i,j)
               end do
               call ev_compute_expr(inputs,result(i),ltest,ierr)
               if( ierr .ne. EVAL_SUCCESS ) then
                  deallocate(ev_differentiate_x)
                  deallocate(ev_differentiate_y)
                  goto 10000
               endif
            end do
            ev_first_pass          = .FALSE.
            ev_differentiate_count = 0
            ev_integrate_count     = 0
            ev_integrate_sum       = 0
            ev_integrate_lastx     = 0
            ev_integrate_lasty     = 0
         end do
         deallocate(ev_differentiate_x)
         deallocate(ev_differentiate_y)
      else
         do i=1,ival
            ev_count = i
            do j=1,nvars
               inputs(j) = values(i,j)
            end do
            call ev_compute_expr(inputs,result(i),ltest,ierr)
            if( ierr .ne. EVAL_SUCCESS ) goto 10000
         end do
      endif
c
10000 continue
      return
      end
c
      subroutine ev_expression_f(expr,var,nvars,ival,values,isize,
     1                           result,ierr,ind)
      implicit none
c
c     Inputs:
c
c     Version of ev_expression with single precision input variables
c
c     expr    -- Expression to evaluate
c     var     -- Array of variable names
c     nvar    -- Number of variables
c     ival    -- Number of times to evaluate the expression using different
c                input values
c     values  -- Input values
c                value(1,1) is 1st data value for variable "1"
c                value(2,1) is 2nd data value for variable "1"
c                value(3,1) is 3rd data value for variable "1"
c                ...
c                value(1,2) is 1st data value for variable "2"
c                value(2,2) is 2nd data value for variable "2"
c                value(3,2) is 3rd data value for variable "2"
c                ...
c     isize   -- As declared size of first dimension of array 'values'
c
c     Outputs:
c
c     result  -- Output vector of size ival
c     ierr    -- Error indicator, see evaluate.cmn for error codes.
c                Successful completion is EVAL_SUCCESS (zero).
c     ind     -- Index into expr() where syntax error occurred.
c
      character*(*) expr
      integer nvars, ival, ierr, ind, isize
c     Actual second dimension of values is "nvars", but we don't code
c     this because nvars could be zero.
      real    values(isize,*), result(ival)
c     Actual dimension of var is "nvars", but we don't code this because
c     nvars could be zero.
      character*(*) var(*)
c
      include 'evaluate.cmn'
c
      logical ltest
      integer ilen, i, j, k, string_len
      double precision rdum1, rdum2, inputs(EVAL_MAXVARS), dresult
c
      ierr = EVAL_SUCCESS
      ind  = 0
c
      call ev_init_common
      call ev_initialize
c
c     Copy expression and variable names into local storage, upcase.
      ilen = string_len(expr)
      if( ilen .gt. EVAL_MAXEXPR ) then
         ierr = EVAL_LARGE
         goto 10000
      endif
      ev_expr = expr(1:ilen)
      call upcase(ev_expr)
c
      ev_num_vars = nvars
      do 100 i=1,nvars
         call ev_check_varname(var(i),i-1,ierr)
         if( ierr .ne. EVAL_SUCCESS ) goto 10000
         ev_var(i)  = var(i)
         call upcase(ev_var(i))
100   continue
c
c     Break expression into tokens
      call ev_tokenize(ierr,ind)
      if( ierr .ne. EVAL_SUCCESS ) goto 10000
c
c     Convert to reverse polish notation
      call ev_rpn(ierr)
      if( ierr .ne. EVAL_SUCCESS ) goto 10000
c
c     Perform a "dry run" of the expresion evaluator, in order to catch
c     syntax errors.
      ltest = .TRUE.
      call ev_compute_expr(rdum1,rdum2,ltest,ierr)
c
c     Evaluate expression for each input value
      ltest = .FALSE.
c
c     Check for too many integration or differentiation functions
      if( ev_integrate_id .gt. EVAL_STACKSIZE ) then
         ierr = EVAL_TOOMANYINTS
         goto 10000
      else if( ev_differentiate_id .gt. 1 ) then
         ierr = EVAL_TOOMANYDIFFS
         goto 10000
      else if( ev_differentiate_id .eq. 1 ) then
c        Special case to support 2-pass derivatives
         ev_differentiate_npts  = ival
         allocate(ev_differentiate_x(1:ival),stat=ierr)
         if( ierr .ne. 0 ) then
            ierr = EVAL_ALLOCATE
            goto 10000
         endif
         allocate(ev_differentiate_y(1:ival),stat=ierr)
         if( ierr .ne. 0 ) then
            deallocate(ev_differentiate_x)
            ierr = EVAL_ALLOCATE
            goto 10000
         endif
         do k=1,2
            do i=1,ival
               ev_count = i
               do j=1,nvars
                  inputs(j) = values(i,j)
               end do
               call ev_compute_expr(inputs,dresult,ltest,ierr)
               if( ierr .ne. EVAL_SUCCESS ) then
                  deallocate(ev_differentiate_x)
                  deallocate(ev_differentiate_y)
                  goto 10000
               endif
               result(i)=MIN(EVAL_MAXFLOAT,MAX(-EVAL_MAXFLOAT,dresult))
            end do
            ev_first_pass          = .FALSE.
            ev_differentiate_count = 0
            ev_integrate_count     = 0
            ev_integrate_sum       = 0
            ev_integrate_lastx     = 0
            ev_integrate_lasty     = 0
         end do
         deallocate(ev_differentiate_x)
         deallocate(ev_differentiate_y)
      else
         do i=1,ival
            ev_count = i
            do j=1,nvars
               inputs(j) = values(i,j)
            end do
            call ev_compute_expr(inputs,dresult,ltest,ierr)
            if( ierr .ne. EVAL_SUCCESS ) goto 10000
            result(i) = MIN(EVAL_MAXFLOAT,MAX(-EVAL_MAXFLOAT,dresult))
         end do
      endif
c
10000 continue
      return
      end
c
      subroutine ev_tokenize(ierr,ind)
      implicit none
c
      integer ierr, ind
c
      include 'evaluate.cmn'
c
      integer i, ilen, itok, string_len
      character*1 ch, prev, next
      logical intoken, ev_is_var
c
      ind = 0
      intoken = .FALSE.
      itok    = 0
c
      ilen = string_len(ev_expr)
      i = 1
c     DO WHILE i <= ilen
1000  if( i .gt. ilen ) goto 2000
c
c        Get previous,current,next character
         ch = ev_expr(i:i)
         if( i .eq. 1 ) then
           prev = ' '
         else
           prev = ev_expr(i-1:i-1)
         endif
         if( i .lt. ilen ) then
            next = ev_expr(i+1:i+1)
         else
            next = ' '
         endif
c
c        Check for unrecognized character in input string
         if(
     1   INDEX('0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ().+-*/_ ><=|&,',
     1         ch) .eq. 0 ) goto 10002
c
c        Large switch statement
c        White space
         if( ch .eq. ' ' ) then
            if( intoken ) then
               ev_token_end(itok) = i-1
               call ev_eval_token(itok,ierr,ind)
               if( ierr .ne. EVAL_SUCCESS ) goto 10000
               intoken = .FALSE.
            endif
c        Left paren
         else if( ch .eq. '(' ) then
            if( intoken ) then
               ev_token_end(itok) = i-1
               call ev_eval_token(itok,ierr,ind)
               if( ierr .ne. EVAL_SUCCESS ) goto 10000
               intoken = .FALSE.
            endif
            itok = itok + 1
            if( itok .gt. EVAL_MAXTOKENS ) goto 10001
            ev_token_start(itok) = i
            ev_token_end(itok) = i
            ev_token_precedence(itok) = EVAL_PREC_NUL
            ev_token_type(itok) = EVAL_TOKEN_OPERATOR
c        Right paren
         else if( ch .eq. ')' ) then
            if( intoken ) then
               ev_token_end(itok) = i-1
               call ev_eval_token(itok,ierr,ind)
               if( ierr .ne. EVAL_SUCCESS ) goto 10000
               intoken = .FALSE.
            endif
            itok = itok + 1
            if( itok .gt. EVAL_MAXTOKENS ) goto 10001
            ev_token_start(itok) = i
            ev_token_end(itok) = i
            ev_token_precedence(itok) = EVAL_PREC_NUL
            ev_token_type(itok) = EVAL_TOKEN_OPERATOR
c        * (multiply)
         else if( ch .eq. '*' ) then
            if( intoken ) then
               ev_token_end(itok) = i-1
               call ev_eval_token(itok,ierr,ind)
               if( ierr .ne. EVAL_SUCCESS ) goto 10000
               intoken = .FALSE.
            endif
            itok = itok + 1
            if( itok .gt. EVAL_MAXTOKENS ) goto 10001
            ev_token_start(itok) = i
c           Special case for '**'
            if( next .eq. '*' ) then
               i = i + 1
               ev_token_precedence(itok) = EVAL_PREC_POW
            else
               ev_token_precedence(itok) = EVAL_PREC_MULDIV
            endif
            ev_token_end(itok) = i
            ev_token_type(itok) = EVAL_TOKEN_OPERATOR
c        / (divide)
         else if( ch .eq. '/' ) then
            if( intoken ) then
               ev_token_end(itok) = i-1
               call ev_eval_token(itok,ierr,ind)
               if( ierr .ne. EVAL_SUCCESS ) goto 10000
               intoken = .FALSE.
            endif
            itok = itok + 1
            if( itok .gt. EVAL_MAXTOKENS ) goto 10001
            ev_token_start(itok) = i
            ev_token_end(itok) = i
            ev_token_precedence(itok) = EVAL_PREC_MULDIV
            ev_token_type(itok) = EVAL_TOKEN_OPERATOR
c        + (add)
         else if( ch .eq. '+' ) then
c           Special check for exponential notation
            if( intoken .and. (prev .eq. 'E' .or. prev .eq. 'D') .and.
     1      .NOT. ev_is_var(ev_expr(ev_token_start(itok):i-1)) ) then
               continue
            else
               if( intoken ) then
                  ev_token_end(itok) = i-1
                  call ev_eval_token(itok,ierr,ind)
                  if( ierr .ne. EVAL_SUCCESS ) goto 10000
                  intoken = .FALSE.
               endif
               itok = itok + 1
               if( itok .gt. EVAL_MAXTOKENS ) goto 10001
               ev_token_start(itok) = i
               ev_token_end(itok) = i
               ev_token_precedence(itok) = EVAL_PREC_ADDSUB
               ev_token_type(itok) = EVAL_TOKEN_OPERATOR
            endif
c        - (subtract)
         else if( ch .eq. '-' ) then
c           Special check for exponential notation
            if( intoken .and. (prev .eq. 'E' .or. prev .eq. 'D') .and.
     1      .NOT. ev_is_var(ev_expr(ev_token_start(itok):i-1)) ) then
               continue
            else
               if( intoken ) then
                  ev_token_end(itok) = i-1
                  call ev_eval_token(itok,ierr,ind)
                  if( ierr .ne. EVAL_SUCCESS ) goto 10000
                  intoken = .FALSE.
               endif
               itok = itok + 1
               if( itok .gt. EVAL_MAXTOKENS ) goto 10001
               ev_token_start(itok) = i
               ev_token_end(itok) = i
               ev_token_precedence(itok) = EVAL_PREC_ADDSUB
               ev_token_type(itok) = EVAL_TOKEN_OPERATOR
            endif
c        Comma, treated as an 'operator' for now
         else if( ch .eq. ',' ) then
            if( intoken ) then
               ev_token_end(itok) = i-1
               call ev_eval_token(itok,ierr,ind)
               if( ierr .ne. EVAL_SUCCESS ) goto 10000
               intoken = .FALSE.
            endif
            itok = itok + 1
            if( itok .gt. EVAL_MAXTOKENS ) goto 10001
            ev_token_start(itok) = i
            ev_token_end(itok) = i
            ev_token_precedence(itok) = EVAL_PREC_NUL
            ev_token_type(itok) = EVAL_TOKEN_OPERATOR
c        > (greater than), < (less than)
         else if( ch .eq. '>' .or. ch .eq. '<' ) then
            if( intoken ) then
               ev_token_end(itok) = i-1
               call ev_eval_token(itok,ierr,ind)
               if( ierr .ne. EVAL_SUCCESS ) goto 10000
               intoken = .FALSE.
            endif
            itok = itok + 1
            if( itok .gt. EVAL_MAXTOKENS ) goto 10001
            ev_token_start(itok) = i
c           Special case for '>=', '<='
            if( next .eq. '=' ) then
               i = i + 1
            endif
            ev_token_end(itok) = i
            ev_token_precedence(itok) = EVAL_PREC_COMPARE
            ev_token_type(itok) = EVAL_TOKEN_OPERATOR
c        && (and), || (or)
         else if( ch .eq. '&' .or. ch .eq. '|' ) then
            if( intoken ) then
               ev_token_end(itok) = i-1
               call ev_eval_token(itok,ierr,ind)
               if( ierr .ne. EVAL_SUCCESS ) goto 10000
               intoken = .FALSE.
            endif
            itok = itok + 1
            if( itok .gt. EVAL_MAXTOKENS ) goto 10001
            ev_token_start(itok) = i
c           Look for next character for '&&', '||'
            if( next .eq. '&' .or. next .eq. '|' ) then
               i = i + 1
               ev_token_precedence(itok) = EVAL_PREC_ANDOR
            else
               ev_token_precedence(itok) = EVAL_PREC_MULDIV
            endif
            ev_token_end(itok) = i
            ev_token_type(itok) = EVAL_TOKEN_OPERATOR
         else
c           What is left are: variables, functions and numbers
            if( .NOT. intoken ) then
               itok = itok + 1
               if( itok .gt. EVAL_MAXTOKENS ) goto 10001
               ev_token_start(itok) = i
               intoken = .TRUE.
            endif
         endif
      i = i + 1
      goto 1000
c
c     END DO
c
2000  continue
c
c     Flush any remaining tokens
      if( intoken ) then
         ev_token_end(itok) = i-1
         call ev_eval_token(itok,ierr,ind)
         if( ierr .ne. EVAL_SUCCESS ) goto 10000
         intoken = .FALSE.
      endif
c
      ev_num_tokens = itok
c
c     Look for unary + and -, and set precedence
      call ev_look_for_unary(ierr)
c
c     Perform some rudimentary syntax checks
c
c     Ensure any function is followed by a parenthesis
      do 3000 i=1,ev_num_tokens
         if( ev_token_type(i) .eq. EVAL_TOKEN_FUNCTION ) then
            if( ev_num_tokens .gt. 1 ) then
c              Look at next token
               ch = ev_expr(ev_token_start(i+1):ev_token_start(i+1))
               if( ch .ne. '(' ) goto 10003
            else
c              Error: function alone on line
               goto 10003
            endif
         endif
3000  continue
c
c     Successful tokenization
c
      ind = 0
c
      goto 10000
c
c     Error trap: too many tokens
10001 ierr = EVAL_TOKEN
      goto 10000
c
c     Error trap: garbage character in string
10002 ierr = EVAL_BADCHAR
      ind  = i
      goto 10000
c
c     Error trap: syntax error
10003 ierr = EVAL_BADSYNTAX
      ind  = ev_token_start(i)
      goto 10000
c
10000 continue
      return
      end
c
      subroutine ev_rpn(ierr)
      implicit none
c
c     Convert tokenized expression to reverse polish notation
c
      integer ierr
c
      include 'evaluate.cmn'
c
      integer EV_NULL
      parameter (EV_NULL = 0)
c
      integer i, itok, iprec, iop, token_type, token_prec,
     1        curfun, funlist(EVAL_MAXTOKENS)
      character*1 token_ch
      logical in_function
c
c     Initialize expression and operator stacks
      ev_expr_list_ptr = 0
      ev_expr_list(1)  = EV_NULL
      ev_op_list_ptr = 0
      ev_op_list(1)  = EV_NULL
c
c     Initialize function argument count machinery
      curfun = 0
      in_function = .FALSE.
      do i=1,EVAL_MAXTOKENS
         funlist(i) = 0
      end do
c
c     Do while for all tokens
      do 1000 itok=1,ev_num_tokens
         token_type = ev_token_type(itok)
         token_prec = ev_token_precedence(itok)
         token_ch   = ev_expr(ev_token_start(itok):ev_token_start(itok))
c        If token is operand, place on expression list
         if( token_type .eq. EVAL_TOKEN_OPERAND .or.
     1       token_type .eq. EVAL_TOKEN_VARIABLE ) then
            ev_expr_list_ptr = ev_expr_list_ptr + 1
            if( ev_expr_list_ptr .gt. EVAL_STACKSIZE ) goto 10002
            ev_expr_list(ev_expr_list_ptr) = itok
         elseif( token_ch .eq. '(' ) then
c           Place '(' on operator stack (precedence 0)
            ev_op_list_ptr = ev_op_list_ptr + 1
            if( ev_op_list_ptr .gt. EVAL_STACKSIZE ) goto 10002
            ev_op_list(ev_op_list_ptr) = itok
         elseif( token_ch .eq. ')' ) then
c           Pop operators off stack and add to expression list,
c           until a matching '(' is found
100         if( ev_op_list_ptr .gt. 0 ) then
               iop = ev_op_list(ev_op_list_ptr)
               if( ev_expr(ev_token_start(iop):ev_token_start(iop))
     1             .ne. '(' ) then
c                 Pop operators off op stack and onto expression list
                  ev_expr_list_ptr = ev_expr_list_ptr + 1
                  if( ev_expr_list_ptr .gt. EVAL_STACKSIZE ) goto 10002
                  ev_expr_list(ev_expr_list_ptr) = iop
c                 Check for pop of function, handle argument list
                  if( ev_token_type(iop) .eq. EVAL_TOKEN_FUNCTION ) then
                     curfun = curfun - 1
                     if( curfun .eq. 0 ) in_function = .FALSE.
                  endif
                  ev_op_list_ptr = ev_op_list_ptr - 1
               else
c                 Matching '(' found, decrement operator stack ptr
c                 and exit loop.
                  ev_op_list_ptr = ev_op_list_ptr - 1
                  goto 101
               endif
            else
c              Error, no matching '(' found
               goto 10001
            endif
            goto 100
c           Found matching '(' exit
101         continue
         else if ( token_ch .eq. ',' ) then
c           If last token on operator list is a function or operator,
c           (but not left parenthesis) put on expression list
            if(ev_op_list_ptr .gt. 0) then
               iop = ev_op_list(ev_op_list_ptr)
               if( ev_expr(ev_token_start(iop):ev_token_start(iop))
     1                  .ne. '(' ) then
                  ev_expr_list_ptr = ev_expr_list_ptr + 1
                  if( ev_expr_list_ptr .gt. EVAL_STACKSIZE ) goto 10002
                  ev_expr_list(ev_expr_list_ptr) = iop
                  ev_op_list_ptr = ev_op_list_ptr - 1
                  if( ev_token_type(iop) .eq. EVAL_TOKEN_FUNCTION ) then
                     curfun = curfun - 1
                  endif
               endif
            endif
            if( curfun .eq. 0 ) goto 10003
            ev_arg_count(funlist(curfun)) =
     1      ev_arg_count(funlist(curfun)) + 1
         else
            if( token_type .eq. EVAL_TOKEN_FUNCTION ) then
               if( in_function ) then
                  curfun = curfun + 1
                  if ( curfun .gt. EVAL_MAXTOKENS ) goto 10002
               else
                  curfun = 1
               endif
               funlist(curfun) = itok
               ev_arg_count(funlist(curfun)) = 1
               in_function = .TRUE.
            endif
c           Check operator precedence against precedence of top of stack
            if( ev_op_list_ptr .eq. 0 ) then
               iprec = -1
            else
               iprec = ev_token_precedence(ev_op_list(ev_op_list_ptr))
            endif
            if( token_prec .gt. iprec ) then
c              Place on operator stack
               ev_op_list_ptr = ev_op_list_ptr + 1
               if( ev_op_list_ptr .gt. EVAL_STACKSIZE ) goto 10002
               ev_op_list(ev_op_list_ptr) = itok
            else
c              Pop operators onto expression list until the token's
c              precedence is < the precedence of the operator on the
c              top of the stack (or '(' found)
200            ev_expr_list_ptr = ev_expr_list_ptr + 1
               if( ev_expr_list_ptr .gt. EVAL_STACKSIZE ) goto 10002
               ev_expr_list(ev_expr_list_ptr) =
     1                            ev_op_list(ev_op_list_ptr)
               ev_op_list_ptr = ev_op_list_ptr - 1
c              Check for pop of function, handle argument list
               if( ev_token_type(ev_expr_list(ev_expr_list_ptr)) .eq.
     1             EVAL_TOKEN_FUNCTION ) then
                  curfun = curfun - 1
                  if( curfun .eq. 0 ) in_function = .FALSE.
               endif
c              Check for empty list, exit if so
               if( ev_op_list_ptr .eq. 0 ) goto 201
c              Check precedence of operator on top of op list <
c              precedence current operator, exit if so
               iop = ev_op_list(ev_op_list_ptr)
               if( token_prec .gt. ev_token_precedence(iop) ) goto 201
c              Check for '(' on operator list
               if( ev_expr(ev_token_start(iop):ev_token_start(iop))
     1             .eq. '(' ) goto 201
               goto 200
201            continue
c              Push operator onto operator stack
               ev_op_list_ptr = ev_op_list_ptr + 1
               if( ev_op_list_ptr .gt. EVAL_STACKSIZE ) goto 10002
               ev_op_list(ev_op_list_ptr) = itok
            endif
         endif
1000  continue
c
c     Check for any remaining operators on op list, put on expression
c     list
300   if( ev_op_list_ptr .le. 0 ) goto 301
c        Check for '(' on operator list (mismatched parens)
         i = ev_token_start(ev_op_list(ev_op_list_ptr))
         if( ev_expr(i:i) .eq. '(' ) goto 10001
         ev_expr_list_ptr = ev_expr_list_ptr + 1
         if( ev_expr_list_ptr .gt. EVAL_STACKSIZE ) goto 10002
         ev_expr_list(ev_expr_list_ptr) = ev_op_list(ev_op_list_ptr)
         ev_op_list_ptr = ev_op_list_ptr - 1
      goto 300
301   continue
c
      goto 10000
c
c     Error condition: mismatched parenthesis
10001 ierr = EVAL_BADPAREN
      goto 10000
c
c     Error condition: overflow of internal stack
10002 ierr = EVAL_STACK
      goto 10000
c
c     Error condition: syntax error
10003 ierr = EVAL_BADSYNTAX
      goto 10000
c
10000 continue
      return
      end
c
      subroutine ev_compute_expr(value,result,ltest,ierr)
      implicit none
c
      logical ltest
      double precision value(*), result
      integer ierr
c
c     Perform expression evaluation. If logical 'ltest' is .TRUE.
c     then go through motions of expression evaluation, but do not
c     perform any computations (this is for catching syntax errors).
c
      include 'evaluate.cmn'
c
      integer i, j, itok, icode, token_type, ivar
      double precision op1, op2, tmp
      character*2 token_ch
      logical lop1, lop2
      integer iop1, iop2
c
      ev_stack_ptr = 0
      ev_stack(1)  = 0.0
      ev_integrate_id     = 0
      ev_differentiate_id = 0
      tmp = 0.0
c
      do 1000 i=1,ev_expr_list_ptr
         itok = ev_expr_list(i)
         token_type = ev_token_type(itok)
         if( token_type .eq. EVAL_TOKEN_OPERAND ) then
c           Push number onto stack
            call ev_push(ev_token_cvt_real(itok),icode)
            if( icode .ne. EVAL_SUCCESS ) goto 10005
         elseif( token_type .eq. EVAL_TOKEN_VARIABLE ) then
c           Match variable, push value onto variable stack
            ivar = 0
            do 1001 j=1,ev_num_vars
               if( ev_expr(ev_token_start(itok):ev_token_end(itok))
     1             .eq. ev_var(j) ) then
                  ivar = j
                  goto 1002
               endif
1001        continue
1002        continue
c           If no match, then exit with internal error
            if( ivar .eq. 0 ) goto 10008
            tmp = DBLE(value(ivar))
            call ev_push(tmp,icode)
            if( icode .ne. EVAL_SUCCESS ) goto 10005
         elseif( token_type .eq. EVAL_TOKEN_FUNCTION ) then
c           Pop off values from stack onto argument list
c           and pass into function evaluation routine.
            do j=1,ev_arg_count(itok)
               call ev_pop(ev_arg_list(j),icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
            end do
            call ev_func_eval(itok,ev_arg_list,ev_arg_count(itok),
     1                        ltest,icode)
            if( icode .ne. EVAL_SUCCESS ) goto 10002
c           Put return value back on stack
            call ev_push(ev_arg_list(1),icode)
            if( icode .ne. EVAL_SUCCESS ) goto 10005
         elseif( token_type .eq. EVAL_TOKEN_UNARY ) then
c           Decode operator
            token_ch = ev_expr(ev_token_start(itok):ev_token_end(itok))
            if( token_ch .eq. '+' ) then
c                 Unary + is really a NOP
            elseif( token_ch .eq. '-' ) then
c                 Unary -, change sign
                  if( ev_stack_ptr .le. 0 ) goto 10003
                  if( .NOT. ltest) ev_stack(1) = - ev_stack(1)
            else
c              This is not the most reasonable error message,
c              but it should "never" happen.
               goto 10006
            endif
         elseif( token_type .eq. EVAL_TOKEN_OPERATOR ) then
c           Decode operator
            token_ch = ev_expr(ev_token_start(itok):ev_token_end(itok))
            if( token_ch .eq. '+' ) then
               call ev_pop(op2,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               call ev_pop(op1,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               if( .NOT. ltest) tmp = op1 + op2
               call ev_push(tmp,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10005
            else if( token_ch .eq. '-' ) then
               call ev_pop(op2,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               call ev_pop(op1,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               if( .NOT. ltest ) tmp = op1 - op2
               call ev_push(tmp,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10005
            else if( token_ch .eq. '**') then
               call ev_pop(op2,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               call ev_pop(op1,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
c              Check for 'undefined exponentiation'
               if( .NOT. ltest ) then
                  if ( op1 .eq. 0.0D0 .and. op2 .le. 0.0D0 ) goto 10007
c                 Check for integer exponent
                  if( (op2 - INT(op2)) .eq. 0.0D0 ) then
                     call ev_pow_int(op1,op2,tmp)
                  else
                     if( op1 .lt. 0.0D0 ) goto 10007
                     tmp = op1 ** op2
                  endif
               endif
               call ev_push(tmp,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10005
            else if( token_ch .eq. '*' ) then
               call ev_pop(op2,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               call ev_pop(op1,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               if( .NOT. ltest ) tmp = op1 * op2
               call ev_push(tmp,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10005
            else if( token_ch .eq. '/' ) then
               call ev_pop(op2,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               call ev_pop(op1,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
c              Check for possible zero divide
               if( .NOT. ltest ) then
                  if( op2 .eq. 0.0D0 ) goto 10001
                  tmp = op1 / op2
               endif
               call ev_push(tmp,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10005
            else if( token_ch .eq. ',' ) then
c              Comma's should have been stripped by ev_rpn,
c              but check anyway
               goto 10003
            else if( token_ch .eq. '>' ) then
               call ev_pop(op2,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               call ev_pop(op1,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               if( .NOT. ltest ) tmp = op1 .gt. op2
               call ev_push(tmp,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10005
            else if( token_ch .eq. '>=' ) then
               call ev_pop(op2,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               call ev_pop(op1,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               if( .NOT. ltest ) tmp = op1 .ge. op2
               call ev_push(tmp,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10005
            else if( token_ch .eq. '<' ) then
               call ev_pop(op2,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               call ev_pop(op1,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               if( .NOT. ltest ) tmp = op1 .lt. op2
               call ev_push(tmp,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10005
            else if( token_ch .eq. '<=' ) then
               call ev_pop(op2,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               call ev_pop(op1,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               if( .NOT. ltest ) tmp = op1 .le. op2
               call ev_push(tmp,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10005
            else if( token_ch .eq. '|' ) then
               call ev_pop(op2,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               call ev_pop(op1,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               iop1 = op1
               iop2 = op2
               if( .NOT. ltest ) tmp = IOR(iop1,iop2)
               call ev_push(tmp,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10005
            else if( token_ch .eq. '&' ) then
               call ev_pop(op2,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               call ev_pop(op1,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               iop1 = op1
               iop2 = op2
               if( .NOT. ltest ) tmp = IAND(iop1,iop2)
               call ev_push(tmp,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10005
            else if( token_ch .eq. '||' ) then
               call ev_pop(op2,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               call ev_pop(op1,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               if( op1 .eq. 0 ) then
                  lop1 = .FALSE.
               else
                  lop1 = .TRUE.
               endif
               if( op2 .eq. 0 ) then
                  lop2 = .FALSE.
               else
                  lop2 = .TRUE.
               endif
               if( .NOT. ltest ) tmp = lop1 .or. lop2
               call ev_push(tmp,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10005
            else if( token_ch .eq. '&&' ) then
               call ev_pop(op2,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               call ev_pop(op1,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10003
               if( op1 .eq. 0 ) then
                  lop1 = .FALSE.
               else
                  lop1 = .TRUE.
               endif
               if( op2 .eq. 0 ) then
                  lop2 = .FALSE.
               else
                  lop2 = .TRUE.
               endif
               if( .NOT. ltest ) tmp = lop1 .and. lop2
               call ev_push(tmp,icode)
               if( icode .ne. EVAL_SUCCESS ) goto 10005
            else
c              Operator type not found. This shouldn't happen.
               goto 10006
            endif
         else
            goto 10004
         endif
1000  continue
c
c     Now, check to make sure we have only one value remaining on the stack
c     and there are no dangling arguments on arg list
      if( ev_stack_ptr .eq. 1 ) then
         result = ev_stack(1)
      else
         goto 10003
      endif
      goto 10000
c
c     Error condition: zero divide
10001 ierr = EVAL_OVERFLOW
      goto 10000
c
c     Error condition: function evaluation went wrong
10002 ierr = icode
      goto 10000
c
c     Error condition: syntax error
10003 ierr = EVAL_BADSYNTAX
      goto 10000
c
c     Error condition: unknown token type error
10004 ierr = EVAL_BADTOKEN
      goto 10000
c
c     Error condition: stack overflow
10005 ierr = EVAL_STACK
      goto 10000
c
c     Error condition: unknown operator error
10006 ierr = EVAL_BADOPERATOR
      goto 10000
c
c     Error condition: bad exponentiation
10007 ierr = EVAL_BADEXP
      goto 10000
c
c     Error condition: cannot match token with variable name
10008 ierr = EVAL_BADVAR
      goto 10000
c
10000 continue
      return
      end
c
      subroutine ev_eval_token(itok,ierr,ind)
      implicit none
c
      integer itok, ierr, ind
c
      include 'evaluate.cmn'
c
c     Subroutine tries to evaluate a token from the class:
c
c     variable
c     function
c     real number
c
      integer i, iend, ioerr
      logical found
      double precision tmp
      character*80 tmpstr   ! Temporary hack for DEC FORTRAN... jvf
c
c     Set up error indicator
      ind = ev_token_start(itok)
      iend = ev_token_end(itok)
c
c     If token starts with a letter, try a function or variable name
      if( INDEX('ABCDEFGHIJKLMNOPQRSTUVWXYZ',
     1          ev_expr(ind:ind) ) .ne. 0 ) then
c        Check functions first
         found = .FALSE.
         do 100 i=1,EVAL_MAXFUN
            if( ev_expr(ind:iend) .eq. ev_function_names(i) ) then
               found = .TRUE.
               ev_token_precedence(itok) = EVAL_PREC_FUNC
               ev_token_type(itok)       = EVAL_TOKEN_FUNCTION
               goto 200
            endif
100      continue
200      continue
         if( .NOT. found ) then
c           Next, check for variable name
            do 300 i=1,ev_num_vars
               if( ev_expr(ind:iend) .eq. ev_var(i) ) then
                  ev_token_precedence(itok) = EVAL_PREC_NUL
                  ev_token_type(itok)       = EVAL_TOKEN_VARIABLE
                  found = .TRUE.
                  goto 301
               endif
300         continue
301         continue
c           Check for unknown function/variable name
            if( .NOT. found ) ierr = EVAL_BADFUN
         endif
      else
c        Next, try for a real number
c        Temporary hack for DEC FORTRAN compiler
c        Program incorrectly reads
c        from ev_expr(ind:iend), and reports an input conversion error.
c        Bug fixed in V3.4 on OSF/1,Ultrix, 6.2 VMS
         tmpstr = ev_expr(ind:iend)
         read(tmpstr,*,end=1000,err=1000,iostat=ioerr) tmp
         ev_token_precedence(itok) = EVAL_PREC_NUL
         ev_token_type(itok)       = EVAL_TOKEN_OPERAND
         ev_token_cvt_real(itok)   = tmp
      endif
c
      return
c
c     Error trap, bad real number
1000  continue
      ierr = EVAL_BADCON
      return
      end
c
      subroutine ev_look_for_unary(ierr)
      implicit none
c
      integer ierr
c
      include 'evaluate.cmn'
c
      integer i
      character*1 prev, cur
c
      do 100 i=1,ev_num_tokens
         cur = ev_expr(ev_token_start(i):ev_token_end(i))
         if( cur .ne. '+' .and. cur .ne. '-' ) goto 100
         if( i .eq. 1 ) then
c           If + or - operator is first token then it is unary
            ev_token_type(i) = EVAL_TOKEN_UNARY
            ev_token_precedence(i) = EVAL_PREC_UNARY
c           Check for next token an operand, if so, increase prec
            if( i+1 .le. ev_num_tokens) then
               if( ev_token_type(i+1) .eq. EVAL_TOKEN_OPERAND ) then
                  ev_token_precedence(i) = EVAL_PREC_FUNC
               endif
            endif
         else
            prev = ev_expr(ev_token_start(i-1):ev_token_start(i-1))
            if( prev .eq. '(' .or. prev .eq. '+' .or. prev .eq. '-' .or.
     1          prev .eq. '*' .or. prev .eq. '/' .or. prev .eq. ',' .or.
     1          prev .eq. '&' .or. prev .eq. '|' .or. prev .eq. '<' .or.
     1          prev .eq. '>' )
     1          then
               ev_token_type(i) = EVAL_TOKEN_UNARY
               ev_token_precedence(i) = EVAL_PREC_UNARY
c              Check for next token an operand, if so, increase prec
               if( i+1 .le. ev_num_tokens) then
                  if( ev_token_type(i+1) .eq. EVAL_TOKEN_OPERAND ) then
                     ev_token_precedence(i) = EVAL_PREC_FUNC
                  endif
               endif
            endif
         endif
100   continue
c
      ierr = EVAL_SUCCESS
      return
      end
c
      subroutine ev_init_common
      implicit none
c
      include 'evaluate.cmn'
c
      logical lfirst
      data    lfirst /.TRUE./
      save    lfirst
c
      if( .NOT. lfirst ) return
      lfirst = .FALSE.
c
c     Names *MUST* be in upper case
c     If two names start with similar letters, then the longer one must
c     come first.
c     If you want to add a function name, change EVAL_MAXFUN in evaluate.cmn
c     modify ev_func_eval to recognize this function, and recompile code.
c     Function names *MUST* be larger than one character long
c     otherwise they would conflict with variable names
c     (NOTE: this is probably no longer true, now that we have
c     multi-character variables. jvf)
      ev_function_names(1) = 'EXP'
      ev_function_names(2) = 'SQRT'
      ev_function_names(3) = 'LOG10'
      ev_function_names(4) = 'LOG'
      ev_function_names(5) = 'ABS'
      ev_function_names(6) = 'SIN'
      ev_function_names(7) = 'COS'
      ev_function_names(8) = 'TAN'
      ev_function_names(9) = 'ASIN'
      ev_function_names(10) = 'ACOS'
      ev_function_names(11) = 'ATAN2'
      ev_function_names(12) = 'ATAN'
      ev_function_names(13) = 'MAX'
      ev_function_names(14) = 'MIN'
      ev_function_names(15) = 'SIND'
      ev_function_names(16) = 'COSD'
      ev_function_names(17) = 'TAND'
      ev_function_names(18) = 'RAN'
      ev_function_names(19) = 'INTEGRATE'
      ev_function_names(20) = 'DIFFERENTIATE'
      ev_function_names(21) = 'DERIV'
      ev_function_names(22) = 'MERGE'
      ev_function_names(23) = 'STEP'
c
      return
      end
c
      subroutine ev_initialize
      implicit none
c
      include 'evaluate.cmn'
c
c     Some quantities (mostly arrays) really don't need to be initialized
c     but we do in order to make a clean start (helps during debug, too).
c
      integer i
c
      do 100 i=1,EVAL_MAXTOKENS
         ev_token_start(i)    = 0
         ev_token_end(i)      = 0
         ev_token_precedence(i) = 0
         ev_token_type(i)     = 0
         ev_arg_count(i)      = 0
         ev_token_cvt_real(i) = 0.0
100   continue
c
      do 200 i=1,EVAL_STACKSIZE
         ev_expr_list(i) = 0
         ev_op_list(i)   = 0
         ev_stack(i)     = 0.0
         ev_arg_list(i)  = 0.0
         ev_integrate_count(i)= 0
         ev_integrate_sum(i)  = 0.0
         ev_integrate_lastx(i)= 0.0
         ev_integrate_lasty(i)= 0.0
         ev_differentiate_count(i)= 0
         ev_differentiate_lastx(i)= 0.0
         ev_differentiate_lasty(i)= 0.0
200   continue
c
      ev_expr_list_ptr   = 0
      ev_op_list_ptr     = 0
      ev_count           = 0
c
      ev_first_pass          = .TRUE.
      ev_differentiate_id    = 0
      ev_integrate_id        = 0
      ev_differentiate_npts  = 0
c
      return
      end
c
      subroutine ev_lookup_error_message(ierr,buffer)
      implicit none
c
      integer ierr
      character*(*) buffer
c
      include 'evaluate.cmn'
c
c     Returns an appropriate error message for each EVAL* error
c
c     If new error codes are 
      if( ierr .eq. EVAL_SUCCESS ) then
         buffer = 'successful completion'
      elseif( ierr .eq. EVAL_BADPAREN ) then
         buffer = 'unbalanced parenthesis'
      elseif( ierr .eq. EVAL_BADFUN ) then
         buffer = 'unknown function or variable name'
      elseif( ierr .eq. EVAL_BADVAR ) then
         buffer = 'illegal variable name'
      elseif( ierr .eq. EVAL_BADCON ) then
         buffer = 'illegal constant'
      elseif( ierr .eq. EVAL_BADSYNTAX ) then
         buffer = 'syntax error'
      elseif( ierr .eq. EVAL_OVERFLOW ) then
         buffer = 'expression caused floating overflow'
      elseif( ierr .eq. EVAL_BADARG ) then
         buffer = 'argument to function was invalid'
      elseif( ierr .eq. EVAL_BADCHAR ) then
         buffer = 'illegal character in expression'
      elseif( ierr .eq. EVAL_BADEXP ) then
         buffer = 'invalid exponentiation'
      elseif( ierr .eq. EVAL_BADVARNAME ) then
         buffer = 'invalid variable name'
      elseif( ierr .eq. EVAL_WRONG_NUM_ARGS ) then
         buffer = 'wrong number of arguments to function'
      elseif( ierr .eq. EVAL_STACK ) then
         buffer = 'internal: stack overflow'
      elseif( ierr .eq. EVAL_TOKEN ) then
         buffer = 'internal: too many tokens'
      elseif( ierr .eq. EVAL_LARGE ) then
         buffer = 'internal: expression too long'
      elseif( ierr .eq. EVAL_BADTOKEN ) then
         buffer = 'internal: invalid token type found ' //
     1            'during expression evaluation'
      elseif( ierr .eq. EVAL_BADOPERATOR ) then
         buffer = 'internal: invalid operator found ' //
     1            'during expression evaluation'
      elseif( ierr .eq. EVAL_NOVARMATCH ) then
         buffer = 'internal: could not match variable with token ' //
     1            'during expression evaluation'
      elseif( ierr .eq. EVAL_ALLOCATE ) then
         buffer = 'internal: error during ALLOCATE of ' //
     1            'an internal array'
      elseif( ierr .eq. EVAL_TOOMANYDIFFS ) then
         buffer = 'internal: too many deriv() functions used'
      elseif( ierr .eq. EVAL_TOOMANYINTS ) then
         buffer = 'internal: too many integrate() functions used'
      else
         buffer = 'unknown or invalid error code'
      endif
c
      return
      end
c
      subroutine ev_push(op,icode)
      implicit none
c
      double precision op
      integer icode
c
      include 'evaluate.cmn'
c
      integer i
c
c     Push op onto eval_stack, checking for stack overflow in process
      if( ev_stack_ptr+1 .gt. EVAL_STACKSIZE ) then
         icode = EVAL_STACK
         return
      endif
c
      do 100 i=ev_stack_ptr,1,-1
        ev_stack(i+1) = ev_stack(i)
100   continue
      ev_stack_ptr = ev_stack_ptr + 1
      ev_stack(1) = op
      icode = EVAL_SUCCESS
c
      return
      end
c
      subroutine ev_pop(op,icode)
      implicit none
c
      double precision op
      integer icode
c
      include 'evaluate.cmn'
c
      integer i
c
c     Pop op off eval_stack, checking for stack underflow in process
      if( ev_stack_ptr .le. 0 ) then
         icode = EVAL_STACK
         return
      endif
c
      op = ev_stack(1)
      do 100 i=2,ev_stack_ptr
        ev_stack(i-1) = ev_stack(i)
100   continue
      ev_stack_ptr = ev_stack_ptr - 1
      icode = EVAL_SUCCESS
c
      return
      end
c
      subroutine ev_func_eval(itok,op,nop,ltest,icode)
      implicit none
c
      integer itok, icode, nop
      logical ltest
      double precision op(nop)
c
      include 'evaluate.cmn'
c
      integer i, is, ie, itype, iargs, iseed
      logical varargs, lfirst_ran
      double precision delx, dely, dydx, temp,
     1                 delx1, delx2, dely1, dely2, dydx1, dydx2
      data lfirst_ran /.TRUE./
c
c     Decode function name and check for right number of arguments
      is = ev_token_start(itok)
      ie = ev_token_end(itok)
      varargs = .FALSE.
      iargs = 1
      if( ev_expr(is:ie) .eq. ev_function_names(1) ) then
c        EXP
         itype = 1
      else if( ev_expr(is:ie) .eq. ev_function_names(2) ) then
c        SQRT
         itype = 2
      else if( ev_expr(is:ie) .eq. ev_function_names(3) ) then
c        LOG10
         itype = 3
      else if( ev_expr(is:ie) .eq. ev_function_names(4) ) then
c        LOG
         itype = 4
      else if( ev_expr(is:ie) .eq. ev_function_names(5) ) then
c        ABS
         itype = 5
      else if( ev_expr(is:ie) .eq. ev_function_names(6) ) then
c        SIN
         itype = 6
      else if( ev_expr(is:ie) .eq. ev_function_names(7) ) then
c        COS
         itype = 7
      else if( ev_expr(is:ie) .eq. ev_function_names(8) ) then
c        TAN
         itype = 8
      else if( ev_expr(is:ie) .eq. ev_function_names(9) ) then
c        ASIN
         itype = 9
      else if( ev_expr(is:ie) .eq. ev_function_names(10) ) then
c        ACOS
         itype = 10
      else if( ev_expr(is:ie) .eq. ev_function_names(11) ) then
c        ATAN2
         itype = 11
         iargs = 2
      else if( ev_expr(is:ie) .eq. ev_function_names(12) ) then
c        ATAN
         itype = 12
      else if( ev_expr(is:ie) .eq. ev_function_names(13) ) then
c        MAX
         itype = 13
         varargs = .TRUE.
      else if( ev_expr(is:ie) .eq. ev_function_names(14) ) then
c        MIN
         itype = 14
         varargs = .TRUE.
      else if( ev_expr(is:ie) .eq. ev_function_names(15) ) then
c        SIND
         itype = 15
      else if( ev_expr(is:ie) .eq. ev_function_names(16) ) then
c        COSD
         itype = 16
      else if( ev_expr(is:ie) .eq. ev_function_names(17) ) then
c        TAND
         itype = 17
      else if( ev_expr(is:ie) .eq. ev_function_names(18) ) then
c        RAN
         if( lfirst_ran ) then
            lfirst_ran = .FALSE.
            iseed  = 2147483647  ! 2**31 - 1 
         endif
         itype = 18
      else if( ev_expr(is:ie) .eq. ev_function_names(19) ) then
c        INTEGRATE
         itype = 19
         iargs = 2
      else if( ev_expr(is:ie) .eq. ev_function_names(20) .or.
     1         ev_expr(is:ie) .eq. ev_function_names(21) ) then
c        DERIV (DIFFERENTIATE)
         itype = 20
         iargs = 2
      else if( ev_expr(is:ie) .eq. ev_function_names(22) ) then
c        MERGE
         itype = 21
         iargs = 3
      else if( ev_expr(is:ie) .eq. ev_function_names(23) ) then
c        STEP
         itype = 22
         iargs = 2
      else
c        Should never happen, but check for it anyway
         icode = EVAL_BADFUN
         return
      endif
c
      if( iargs .ne. nop .and. .NOT. varargs ) then
         icode = EVAL_WRONG_NUM_ARGS
         return
      endif
c
c     If this is the integrate function, increment the integrate function count
      if( itype .eq. 19 ) then
         ev_integrate_id = ev_integrate_id + 1
      endif
c
c     If this is the differentiate function, increment the deriv function count
      if( itype .eq. 20 ) then
         ev_differentiate_id = ev_differentiate_id + 1
      endif
c
c     if this is a "test" run, do not perform the actual evaluation
      if( ltest ) return
c
c     Decode function name
      if( itype .eq. 1 ) then
c        EXP function
         if( op(1) .gt. EVAL_MAXEXP ) then
            op(1) = EVAL_MAXFLOAT
         else
            op(1) = EXP(op(1))
         endif
         icode = EVAL_SUCCESS
      else if( itype .eq. 2) then
c        SQRT function
         if( op(1) .lt. 0.0 ) then
            op(1) = 0
            icode = EVAL_BADARG
         else
            op(1) = SQRT(op(1))
            icode = EVAL_SUCCESS
         endif
      else if( itype .eq. 3) then
c        LOG10 function
         if( op(1) .le. 0.0 ) then
            op(1) = 0
            icode = EVAL_BADARG
         else
            op(1) = LOG10(op(1))
            icode = EVAL_SUCCESS
         endif
      else if( itype .eq. 4) then
c        LOG function
         if( op(1) .le. 0.0 ) then
            op(1) = 0
            icode = EVAL_BADARG
         else
            op(1) = LOG(op(1))
            icode = EVAL_SUCCESS
         endif
      else if( itype .eq. 5) then
c        ABS function
         op(1) = ABS(op(1))
         icode = EVAL_SUCCESS
      else if( itype .eq. 6) then
c        SIN function
         op(1) = SIN(op(1))
         icode = EVAL_SUCCESS
      else if( itype .eq. 7) then
c        COS function
         op(1) = COS(op(1))
         icode = EVAL_SUCCESS
      else if( itype .eq. 8) then
c        TAN function
         op(1) = TAN(op(1))
         icode = EVAL_SUCCESS
      else if( itype .eq. 9) then
c        ASIN function
         if( op(1) .lt. -1.0D0 .or. op(1) .gt. 1.0D0 ) then
            op(1) = 0
            icode = EVAL_BADARG
         else
            op(1) = ASIN(op(1))
            icode = EVAL_SUCCESS
         endif
      else if( itype .eq. 10) then
c        ACOS function
         if( op(1) .lt. -1.0D0 .or. op(1) .gt. 1.0D0 ) then
            op(1) = 0
            icode = EVAL_BADARG
         else
            op(1) = ACOS(op(1))
            icode = EVAL_SUCCESS
         endif
      else if( itype .eq. 11) then
c        ATAN2 function
         op(1) = ATAN2(op(2),op(1))
         icode = EVAL_SUCCESS
      else if( itype .eq. 12) then
c        ATAN function
         op(1) = ATAN(op(1))
         icode = EVAL_SUCCESS
      else if( itype .eq. 13) then
c        MAX function
         do i=2,nop
            op(1) = MAX(op(1),op(i))
         end do
         icode = EVAL_SUCCESS
      else if( itype .eq. 14) then
c        MIN function
         do i=2,nop
            op(1) = MIN(op(1),op(i))
         end do
         icode = EVAL_SUCCESS
      else if( itype .eq. 15) then
c        SIND function
         op(1) = SIND(op(1))
         icode = EVAL_SUCCESS
      else if( itype .eq. 16) then
c        COSD function
         op(1) = COSD(op(1))
         icode = EVAL_SUCCESS
      else if( itype .eq. 17) then
c        TAND function
c        Some simple firewalling against inifinites produced by tan function
c        Should (of course) reduce argument of tan to [-90,90]
         if( op(1) .eq. 90.0D0 ) then
            op(1) = EVAL_MAXFLOAT
         else if( op(1) .eq. -90.0D0 ) then
            op(1) = -EVAL_MAXFLOAT
         else
            op(1) = TAND(op(1))
         endif
         icode = EVAL_SUCCESS
      else if( itype .eq. 18) then
c        RAN function
         op(1) = RAN(iseed)
         icode = EVAL_SUCCESS
      else if( itype .eq. 19) then
c        INTEGRATE function: INTEGRATE(y,x)
c        op(2) == y == f(x); op(1) == x;
c        (due to RPN evaluation, args pop off right to left)
         i = ev_integrate_id
         ev_integrate_count(i) = ev_integrate_count(i) + 1
         if( ev_integrate_count(i) .gt. 1 ) then
            ev_integrate_sum(i) = ev_integrate_sum(i) +
     1         0.5D0*(op(2)+ev_integrate_lasty(i))*
     1               (op(1) - ev_integrate_lastx(i))
         endif
         ev_integrate_lasty(i) = op(2)
         ev_integrate_lastx(i) = op(1)
         op(1) = ev_integrate_sum(i)
         icode = EVAL_SUCCESS
      else if( itype .eq. 20) then
c        DIFFERENTIATE function: DIFFERENTIATE(y,x)
c        op(2) == y == f(x); op(1) == x;
c        NOTE: we now support a two pass differentiator, so we can
c              do a better job on derivatives.
c        (due to RPN evaluation, args pop off right to left)
         ev_differentiate_count(1) = ev_differentiate_count(1) + 1
         if( ev_first_pass ) then
            ev_differentiate_x((ev_differentiate_count(1))) = op(1)
            ev_differentiate_y((ev_differentiate_count(1))) = op(2)
            if( ev_differentiate_count(1) .gt. 1 ) then
               delx = (op(1) - ev_differentiate_lastx(1))
               dely = (op(2) - ev_differentiate_lasty(1))
               if( delx .eq. 0.0 ) then
                 dydx = SIGN(EVAL_MAXFLOAT,dely)
               else
                 dydx = dely/delx
               endif
            else
               dydx = 0.0
            endif
            ev_differentiate_lasty(1) = op(2)
            ev_differentiate_lastx(1) = op(1)
         else
            i = ev_differentiate_count(1)
            if( ev_differentiate_npts .eq. 1 ) then
               dydx = 0.0D0
            else if( i .eq. 1 ) then
               call ev_deriv(ev_differentiate_x(i),
     1                       ev_differentiate_y(i), 2, dydx)
            else if( i .eq. ev_differentiate_npts ) then
               call ev_deriv(ev_differentiate_x(i-1),
     1                       ev_differentiate_y(i-1), 2, dydx)
            else if( i .eq. 2 .or. i .eq. ev_differentiate_npts-1 ) then
               call ev_deriv(ev_differentiate_x(i-1),
     1                       ev_differentiate_y(i-1), 3, dydx)
            else
               call ev_deriv(ev_differentiate_x(i-2),
     1                       ev_differentiate_y(i-2), 5, dydx)
            endif
         endif
         op(1) = dydx
         icode = EVAL_SUCCESS
      else if( itype .eq. 21) then
c        MERGE(a,b,c)
c        if (c) then return(a) else return(b)
c        (due to RPN evaluation, args pop off right to left)
         if( op(1) .ne. 0 ) then
           temp = op(3)
         else
           temp = op(2)
         endif
         op(1) = temp
         icode = EVAL_SUCCESS
      else if( itype .eq. 22) then
c        STEP(a,b)
c        Stepping (iota) function = a + b*(count-1)
         temp = op(2)
         op(1) = op(2) + op(1)*DBLE(ev_count - 1)
         icode = EVAL_SUCCESS
      endif
c
      return
      end
c
      logical function ev_is_var(string)
      implicit none
c
      character*(*) string
c
c     Returns .TRUE. if string is a user defined variable name
c     Else returns .FALSE.
c
      include 'evaluate.cmn'
c
      integer i
c
      ev_is_var = .FALSE.
      do 100 i=1,ev_num_vars
         if( string .eq. ev_var(i) ) then
            ev_is_var = .TRUE.
            goto 101
         endif
100   continue
101   continue
c
      return
      end
c
      subroutine ev_check_varname(name,ivar,ierr)
      implicit none
c
c     Check variable name for validity
c
      integer ivar, ierr
      character*(*) name
c
      include 'evaluate.cmn'
c
      integer i, ilen, string_len
      character*(EVAL_MAXVARLEN) buff
c
      ilen = string_len(name)
      if( ilen .le. 0 .or. ilen .gt. EVAL_MAXVARLEN) then
         ierr = EVAL_BADVARNAME
         goto 10000
      endif
c
      buff = name
      call upcase(buff)
c
c     Check for invalid characters, must start with a letter
      if(INDEX('ABCDEFGHIJKLMNOPQRSTUVWXYZ_',buff(1:1))
     1      .eq. 0 ) then
         ierr = EVAL_BADVARNAME
         goto 10000
      endif
c     Rest must be alphanumeric
      do 100 i=2,ilen
         if(INDEX('0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ_',buff(i:i))
     1      .eq. 0 ) then
            ierr = EVAL_BADVARNAME
            goto 10000
         endif
100   continue
c
c     Has this name been used before?
      do 200 i=1,ivar
         if (ev_var(i) .eq. buff) then
            ierr = EVAL_BADVARNAME
            goto 10000
         endif
200   continue
c
10000 continue
c
      return
      end
c
      subroutine ev_pow_int(num,exponent,result)
      implicit none
c
      double precision num, exponent, result
c
      integer i, iabs, j
c
      i = INT(exponent)
      iabs = ABS(i)
      result = 1.0D0
      if( i .gt. 0 ) then
         do 100 j=1,iabs
            result = result*num
100      continue
      else
         do 101 j=1,iabs
            result = result/num
101      continue
      endif
c
      return
      end
c
      subroutine ev_deriv(x,y,n,a)
c
      implicit none
c
      integer n
      double precision x(n), y(n), a
c
      include 'evaluate.cmn'
c
c     locals
      integer k
      double precision sumx, sumxy, xmean, ymean
c
      xmean=0.0D0
      do k=1,n
        xmean=xmean+x(k)
      enddo
      xmean=xmean/n
c
      ymean=0.0D0
      do k=1,n
        ymean=ymean+y(k)
      enddo
      ymean=ymean/n
c
      sumx=0.0D0
      do k=1,n
        sumx=sumx+(x(k)-xmean)*(x(k)-xmean)
      enddo
c
      sumxy=0.0D0
      do k=1,n
        sumxy=sumxy+(x(k)-xmean)*(y(k)-ymean)
      enddo
c
      if( sumx .eq. 0.0D0 ) then
         a = SIGN(EVAL_MAXFLOAT,sumxy)
      else
         a = sumxy/sumx
      endif
c
      return
      end
